RNA m6A methylation participates in regulation of postnatal development of the mouse cerebellum

RNA m6A methylation participates in regulation of postnatal development of the mouse cerebellum 6 6 Background: N -methyladenosine (m A) is an important epitranscriptomic mark with high abundance in the brain. Recently, it has been found to be involved in the regulation of memory formation and mammalian cortical neurogenesis. However, while it is now established that m A methylation occurs in a spatially restricted manner, its functions in specific brain regions still await elucidation. Results: We identify widespread and dynamic RNA m A methylation in the developing mouse cerebellum and further uncover distinct features of continuous and temporal-specific m A methylation across the four postnatal developmental processes. Temporal-specific m A peaks from P7 to P60 exhibit remarkable changes in their distribution patterns along the mRNA transcripts. We also show spatiotemporal-specific expression of m A writers METTL3, METTL14, and WTAP and erasers ALKBH5 and FTO in the mouse cerebellum. Ectopic expression of METTL3 mediated by lentivirus infection leads to disorganized structure of both Purkinje and glial cells. In addition, under hypobaric hypoxia exposure, Alkbh5-deletion causes abnormal cell proliferation and differentiation in the cerebellum through disturbing the balance of RNA m A methylation in different cell fate determination genes. Notably, nuclear export of the hypermethylated RNAs is enhanced in the cerebellum of Alkbh5-deficient mice exposed to hypobaric hypoxia. Conclusions: Together, our findings provide strong evidence that RNA m A methylation is controlled in a precise spatiotemporal manner and participates in the regulation of postnatal development of the mouse cerebellum. Keywords: N -methyladenosine, RNA methylation, ALKBH5, METTL3, Cerebellar development Background epitranscriptomic mark found in mRNAs [2, 3]. Similar to Epigenetic regulation, including histone modifications, DNA and protein modifications, m Ain mRNA is DNA modifications, and non-coding RNAs, plays crucial reversible and thus open to dynamic regulation [4, 5]. roles in both embryonic and adult neurogenesis [1]. The Levels of RNA methylation are finely balanced through an 6 6 most recently discovered regulatory modification, N - interplay among m A methyltransferases (writers), methyladenosine (m A), is a highly abundant, so-called demethylases (erasers), and binding proteins (readers) [6]. m A has been shown to regulate RNA processing, including RNA splicing [7, 8], nuclear export [9], RNA * Correspondence: niuym@ibms.pumc.edu.cn; songshh@big.ac.cn; wmtong@ibms.pumc.edu.cn degradation [10, 11], and translation [12–14]. In Chunhui Ma, Mengqi Chang and Hongyi Lv contributed equally to this agreement with its functions in RNA metabolism, m Ais work. 1 an important regulatory factor in different cellular Department of Pathology, Institute of Basic Medical Sciences Chinese Academy of Medical Science, School of Basic Medicine Peking Union Medical processes, including heat shock response [13], DNA College; Neuroscience Center, Chinese Academy of Medical Sciences, Beijing damage response [15], cell fate determination [16–19], 100005, China 2 and innate immunity [20]. In addition to its functions in BIG Data Center, Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing 100101, China these cellular processes, in vivo studies in model Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Ma et al. Genome Biology (2018) 19:68 Page 2 of 18 6 6 organisms have revealed crucial roles of m Ain m A peaks (Additional file 1: Figure S1d). Distribution spermatogenesis [9, 21–23], embryogenesis [24, 25], and analysis of all the m A peaks in mRNAs revealed that T-cell homeostasis [26]. their methylation patterns were similar from P7 to P60 Although m A is abundant in the brain, its dynamic (Additional file 1: Figure S1e). Intriguingly, we observed regulation and biological significance for brain development significant enrichment of m A peaks in start codon remain largely unknown [3]. Among the m A-related genes regions in addition to the coding sequence (CDS) and identified, the RNA demethylase FTO is proven to be stop codon regions. A statistical analysis of the m A peaks crucial for various brain functions, including D2R-D3R- located in different transcript segments showed that the GIRK-mediated signaling [27], adult neurogenesis [28], and proportion of m A peaks in start codon regions increased memory formation [29, 30]. A recent report confirms the gradually from P7 to P60 (Additional file 1: Figure S1f and role of m A methylation in temporal progression of Additional file 1: Table S3). mammalian cortical neurogenesis and in regulating axon To explore how m A changes in postnatal mouse regeneration in the adult mammalian nervous system [31, cerebellum, m A peaks were subjected to pairwise 32]. To better understand the underlying mechanisms comparison between each two adjacent stages. In total, responsible for the temporal and spatial complexity of the we identified 12,452 “ON” switches (emerging m A brain, it is important to characterize the physiological peaks in the later stage) and 11,192 “OFF” switches 6 6 functions of m A in a context-dependent manner. However, (resolving m A peaks from the former to the later stage) the biological functions of other m A-related genes in (P < 0.05, enrichment score > 1.5) during the postnatal different brain regions remain unclear. development (Fig. 1a, b and Additional file 1: Figure The cerebellum consists of a multi-layered structure, S1g). The m A ON and OFF switches also exhibited rendering it an important model for the study of neur- apparent differences in their distribution patterns. The onal cell proliferation, migration, and differentiation m A ON switches were more likely to be enriched in [33]. Considerable progress has been made in identifying start codon regions rather than stop codon regions. In the epigenetic mechanisms responsible for cerebellar de- contrast, the numbers of m A OFF switches velopment, such as DNA 5-hydroxymethylation [34] and surrounding stop codon regions were higher than that of chromatin accessibility [35]. Interestingly, the cerebel- ON switches (Fig. 1c and Additional file 1: Figure S1h). lum of Fto-deficient mouse was found to be smaller than Together, these data imply that m A marks occurring at that of wild-type mouse [28]; however, mechanisms re- earlier or later stages of cerebellar development have sponsible for this growth defect remain unclear. In different preferences with regard to their distribution agreement with the spatial and temporal regulation in along with the whole mRNA transcripts. the brain, cerebellum-specific methylation has been To evaluate the biological significance of genes with dy- identified [36]. Despite these findings, the precise role of namic RNA m A modification, we next performed Gene RNA methylation in the cerebellum and its temporal Ontology (GO) analysis for those genes with ON/OFF control during postnatal mouse cerebellar development switches (Additional file 2: Table S4). The genes with still await elucidation. different types of switches appeared to possess different Here we investigated the dynamic regulation of m Ain functions. For example, the genes containing P7–P14 OFF the developing mouse cerebellum and determined its switches were mostly annotated to biological processes biological significance by disturbing the expression of the such as cell cycle, cell division, and DNA repair (Fig. 1d). 6 6 m AwriterMETTL3orthe m AeraserALKBH5in mice In contrast, the newly occurring methylation at P14, P21, exposed to hypobaric hypoxia. This study uncovered that the and P60 represented by the ON switches was detected in temporal regulation of RNA m A methylation is essential for many genes involved in signal transduction, cell adhesion, precise control of postnatal cerebellar development. learning, and synaptic plasticity (Fig. 1e). Together, these data confirm that extensive RNA methylation and Results demethylation occur over the course of neuronal Temporal regulation of m A methylation in postnatal differentiation in vivo, in both proliferating and fully mouse cerebellum differentiated neural cells. Moreover, the distinct functions 6 6 To understand the biological impact of m A in mouse of genes containing m A OFF or ON switches suggest 6 6 cerebellum, we first performed m Aanalysisusing that m A is a prerequisite for those genes to exert their postnatal cerebella at day 7 (P7), P14, P21, and P60 functions at each developmental stage. (Additional file 1: Table S1). We identified 10,449, 10,266, 10,419, and 10,783 methylated poly(A) RNAs at the four Temporal-specific m A methylation acts in concert with stages, respectively (Additional file 1: Table S2, Additional cerebellar developmental control file 1: Figure S1a–c). Data analysis showed that UGGACU As m A is developmentally regulated in mouse cerebella, was the most conserved consensus sequence among all we next analyzed the methylation profiles over the four Ma et al. Genome Biology (2018) 19:68 Page 3 of 18 Fig. 1 Developmentally regulated RNA methylation during mouse cerebellar development. a Numbers of ON or OFF m A switches between P7, P14, P21, and P60 cerebella pairwise comparisons. The numbers of the poly(A) RNAs containing the m A switches are included in parentheses. b IGV plots 6 6 showing examples of ON and OFF m A switches at each transition period. Grey reads originate from input libraries and red reads originate from m A- IP libraries. Y-axis represents normalized numbers of reads count. Arrows indicate the direction of gene transcription. Dashed boxes indicate the position of ON/OFF switches. c Distribution of each type of ON and OFF m A switch along the whole mRNA transcripts. d, e Most impacted GO biological process terms of the methylated RNAs containing OFF (d)or ON(e)m A switches across the four developmental stages postnatal stages. We found 8367 continuously methylated direction (Fig. 2c and Additional file 1: Figure S2b–c). RNAs (CMRs) throughout cerebellar development, Given the different distribution between ON and OFF together with 634, 260, 315, and 512 specifically switches (Fig. 1c and Additional file 1: Figure S1h), we methylated RNAs (SMRs) at P7, P14, P21, and P60, further analyzed the distribution of m A peaks in respectively (Fig. 2a, b and Additional file 1: Figure S2a). continuously and temporal-specifically methylated We then investigated whether the SMRs and CMRs mRNAs. In the mouse cerebellum of P7, the temporal- possessed different features in their methylation during specific m A peaks were significantly enriched at stop postnatal cerebellar development. Compared with the codon regions, while the P60-specific m A peaks were SMRs, the CMRs exhibited higher levels of both mostly enriched in start codon regions. The P14- or P21- methylation and expression throughout the specific m A peaks displayed an apparent transition from developmental process. Moreover, the methylation levels stop codon to start codon during progression of cerebellar of SMRs displayed a gradual reduction from P7 to P60, development (Fig. 2d). Furthermore, we quantified the while their expression levels changed in the opposite numbers of temporal-specific m A peaks located in the Ma et al. Genome Biology (2018) 19:68 Page 4 of 18 Fig. 2 Distinct features of temporal-specific m A methylation during postnatal development of the mouse cerebellum. a Venn diagram showing the numbers and relationship of methylated poly(A) RNAs in the mouse cerebellum from P7 to P60. b IGV plots showing examples of continuous and temporally specific methylation. Grey reads originate from input libraries and red reads originate from m A-IP libraries. Y-axis represents normalized numbers of read counts. Arrows indicate the direction of gene transcription. Dashed boxes indicate the position of temporal-specific 6 6 m A peaks. c Box plots showing the relative methylation levels of SMRs and CMRs as evaluated by the enrichment scores of m A peaks. SMR specifically methylated RNA, CMR continuously methylated RNA. ***P < 0.001 by Wilcoxon test. d Normalized distribution of temporal-specific m A peaks in each developmental stage along the whole mRNA transcripts. e, f Most impacted GO biological process terms (BP)(e) and cellular component terms (CC)(f) of SMRs in the four developmental stages five regions of mRNA transcripts. From P7 to P60, we the cells in the cerebellum at P60 are fully differentiated. observed a sharp increase in the numbers of temporal- To confirm that such a dynamic distribution pattern was specific m A peaks in start codon regions, as well as a related to the cell states, we compared the RNA methyla- decrease in the CDS and stop codon regions (Additional tion profiles in mouse cerebellum at P7 and P60 directly file 1: Figure S2d). To ensure that these changes were real using the peaks identified by exomePeak. Although the rather than artifacts, we randomly extracted the same proportion of m A peaks located in CDS were similar numbers of CMR peaks as those of temporal-specific between the two stages, we observed a higher percentage peaks and repeated the distribution analysis. Unlike the of P7-specific m A peaks near stop codon regions, but a temporal-specific peaks shown in Fig. 2d, these CMR higher percentage of P60-specific m A peaks in start peaks displayed similar distribution patterns from P7 to codon regions (Additional file 1: Figure S2f). These data P60 (Additional file 1: Figure S2e). The majority of cells in show that for those temporal-specifically methylated the cerebellum of P7 are proliferating and immature, while RNAs, the patterns of m A deposition onto RNA exhibit Ma et al. Genome Biology (2018) 19:68 Page 5 of 18 significant changes along with the progression of (Additional file 1: Figure S4c), among which 839 RNAs cerebellar development. showed either positive (674) or negative (165) correlation To explore the functional significance of temporal- between methylation and expression levels across the four specific methylation, we next performed GO analysis for stages (Additional file 5: Table S7). Results of the cluster those genes with continuous or temporal-specific methy- and subsequent GO analyses suggest that for these RNAs, lation. In line with the progression of cerebellar develop- m A participates in the developmental control of the ment, the four groups of genes containing temporal- mouse cerebellum through modulating their gene specific m A marks were annotated to different biological expression (Additional file 1: Figure S4d, e). processes (Fig. 2e and Additional file 3: Table S5). Genes with P7-specific methylation (such as Aspm in Fig. 2b) Altered expression of either METTL3 or ALKBH5 causes were enriched in processes such as cell cycle, cell division, defective mouse cerebellar development and DNA damage response, which are prerequisites in To gain insights into how the dynamic m A methyl marks proliferating neural cells. In contrast, many genes display- are regulated during cerebellar development, we analyzed ing P60-specific methylation (such as Glra1 in Fig. 2b) the protein expression profiles of three m A writers were annotated to biological processes including trans- (METTL3, METTL14, and WTAP) and two m A erasers port, oxidation-reduction process, and metabolic pro- (ALKBH5 and FTO). All five proteins were highly cesses, which are required for mature neuronal activities. expressed at the early stage of cerebellar development (P7) In agreement with their diverging biological functions, the and showed a gradual reduction towards the maturation genes with temporal-specific methylation were also char- of cerebellar neurons (P60) (Fig. 3a). Given the acteristic of their distinct cellular localization (Fig. 2f). In spatiotemporal specificity of gene expression in the brain, contrast, genes encoded by the P7 and P60 CMRs exhib- we performed immunohistochemical staining to detect ited high similarity in their annotation of enriched bio- their expression in situ. In the cerebellum at P7, all five logical processes and cellular components (Additional file proteins were detected in the external granule cell layer 1: Figure S2g and Additional file 3: Table S5). Together, (EGL), Purkinje cell layer (PCL), and inner granule cell these results suggest that the majority of genes with con- layer (IGL); however, their expression levels varied among tinuous methylation play a fundamental role throughout different types of cells (Fig. 3b). Furthermore, along with cerebellar development, while temporal-specific methyla- the progression of cerebellar development, we observed a tion only occurs at specified times for those genes, thus reduction of their expression levels in internal granular ensuring the proper progression of cerebellar development layers but an elevation in Purkinje cells (Fig. 3c). in a temporal-dependent manner. In agreement with the dynamic expression of the 6 6 6 In parallel to the methylation analysis using the m A m A-related proteins, the global m A levels of poly(A) peaks identified using exomePeak, we re-analyzed the RNA also decreased from P7 to P60 (Fig. 3d). The 6 6 6 m A-seq data at P7 and P60 using another canonical highest m A level at P7 suggests that m A methylation peak-calling algorithm, MACS2. The P7 SMRs, P60 plays an important role at the earlier developmental SMRs, and CMRs exhibited distinct features in their rela- stage. To confirm this, we reduced local m A levels in tive methylation levels, distribution patterns of the m A the developing cerebellum (P7) by knocking down peaks, and their functional annotations (Additional file 1: Mettl3 using lentivirus system and examined the Figure S3 and Additional file 4: Table S6). The very high consequent morphological changes (Fig. 4a–d). Purkinje similar results from the two sets of independent analysis cells in the control group displayed an orderly alignment further strengthened our awareness of the importance of along the outer surface of IGLs with dendritic temporal specific m A methylation during mouse arborization. However, knockdown of Mettl3 resulted in cerebellar development. We then evaluated how the RNA a severe alteration in Purkinje cell numbers and laminal methylation impacted on gene expression. When each structure and and stunted dendrites (Fig. 4e). Moreover, sample was analyzed individually, RNA methylation and GFAP immunostaining revealed that glial cell fibers expression levels at the global scale were negatively were severely disorganized in Mettl3-knocked down correlated with each other (Additional file 1: Figure S4a). regions (Fig. 4f). We also examined the effect of Mettl3 Next, pairwise comparison of RNA expression between overexpression on cerebellar development (Additional each two adjacent stages identified a total of 2451 (P value file 1: Figure S5a, b). Although the global m A levels < 0.05) differentially expressed RNAs throughout increased modestly (Additional file 1: Figure S5c), development (Additional file 1: Figure S4b), among which overexpression of Mettl3 led to apparent morphology approximately 90% were validated by another parallel changes, as was observed in the Mettl3-knocked down differential expression analysis using a STAR/edgeR cerebellum (Additional file 1: Figure S5d-S5e), suggest- package. Cluster analysis of all expressed RNAs from P7 ing that appropriate m A levels are crucial for cerebellar to P60 produced six clusters of gene expression patterns development. Ma et al. Genome Biology (2018) 19:68 Page 6 of 18 Fig. 3 Dynamic expression of m A-related genes and RNA methylation levels during mouse cerebellar development. a Western blot analysis showing the expression levels of m A methyltransferases and demethylases in the mouse cerebellum from P7 to P60. n = 12 biological replicates. GAPDH is used as an internal control. b, c Immunostaining of mouse cerebellum at P7 and P60 using antibodies against the five m A-related proteins. n = 4 biological replicates. Enlarged images in the dashed boxes are shown in c. Scale bar in b represents 50 μm, and scale bar in c represents 10 μm. EGL external granule cell layer, IGL internal granule cell layer, PCL Purkinje cell layer, ML molecular layer. The dashed circles indicate the representative Purkinje cells. d Dynamic poly(A) RNA m A methylation levels during mouse cerebellar development. Mouse cerebella at P7, P14, P21, and P60 were included in this assay. Experiments were repeated three times using 15 mice in total. Representative data are shown here. ***P < 0.001 Given the spatiotemporal expression of ALKBH5 in the speculated that this phenomenon was likely due to genetic mouse cerebellum, we next tested its role in cerebellar de- redundancy and a compensatory response [37]. To over- velopment. Compared to that of wild-type (WT) mice, the come these compensatory mechanisms, we made use of the cerebellum of Alkbh5-knockout (KO) mice lacked any de- fact that the mammalian brain is sensitive to neuronal dam- tectable changes in weight and morphology (Additional file age by hypoxia [38], while ALKBH5 is unique among the 1: Figure S6). Considering the potential functions of FTO ALKB families in its response to hypoxia stimulation [39]. and other undefined demethylases in the cerebellum, we Therefore, we challenged the mice to extreme physiological Ma et al. Genome Biology (2018) 19:68 Page 7 of 18 Fig. 4 Altered morphology in the mouse cerebellum infected with lentivirus for Mettl3 knockdown. a, b GFP fluorescence image (a) and hematoxylin-eosin (HE) staining (b) of the mouse cerebellum dissected 10-days post-lentivirus infection. Arrows indicate the position of injection. Scale bar in a represents 1 mm. Scale bar in b represents 200 μm. c Western blot analysis to confirm the knockdown efficiency of Mettl3 in the neuro2a cell line. β-ACTIN was used as an internal control. sh-1 and sh-2 represent two kinds of Mettl3 knockdown lentiviral vectors. sh-1 was used for subsequent in vivo analysis. d Decreased m A levels resulting from Mettl3 knocking down as analyzed using UHPLC-MS/MS. n =4. *P < 0.05. e, f Immunofluorescence staining with antibodies against Calbindin-D28K (e) and GFAP (f) was performed to detect Purkinje cells (e) and astrocytes (f) upon lentivirus infection (n = 3 biological replicates). GFP fluorescence indicates the area with Mettl3 knockdown. Sections were counterstained with DAPI to visualize the nuclei. Scale bar represents 100 μm conditions in order to exhaust the potential compensatory Alkbh5-deficient cerebellum was significantly reduced as mechanisms. We propose that ALKBH5 might be involved reflected by NeuN immunostaining (Fig. 5f). These data in the protective mechanism of the brain against damage suggest that Alkbh5 deficiency disrupted the normal caused by hypoxia. To confirm this, WT and KO mice were progression of granular neurons from proliferation to exposed to hypobaric hypoxia for 48 h. We found that in differentiation. In addition, we found that Alkbh5 deficiency KO mice the sizes of the whole brain and the cerebellum considerably reduced dendritic arborization of Purkinje were significantly reduced compared to their littermate cells (Fig. 5g), concomitant with an increase in controls (Fig. 5a, b). Immunostaining analysis revealed a disorganization of the radial fibers in glial cells (Fig. 5h). significant increase in the numbers of Ki67 proliferating Following hypoxic treatment and return to normoxic cells and phospho-histone 3 (PH3 ) mitotic cells in the conditions, we examined the cerebellar sections again at EGL compared to the WT counterparts (Fig. 5c, d). another stage (P14). The internal granule neurons, Purkinje Consistently, we also detected an increase in the number of cells, and glial cells exhibited no visible difference between cells in the S phase of the cell cycle (positive BrdU WT and KO mice (Additional file 1:FigureS7a–c). immunoactivity) in the cerebellum of KO mice (Fig. 5e). However, we could still observe a slight increase in the Furthermore, the number of mature neurons in the IGL of numbers of proliferating granule cells as indicated by Ki67 Ma et al. Genome Biology (2018) 19:68 Page 8 of 18 Fig. 5 Defective cerebellar development in Alkbh5 knockout mice after being exposed to hypobaric hypoxia for 48 h. a Brain weight of individual wild-type (WT) and Alkbh5 knockout (KO) mice. Numbers of biological replicates are included in parentheses. Black lines in the graph indicate the mean. ***P < 0.001. b Representative HE staining images of wild-type (WT) and Alkbh5 knockout (KO) mouse cerebellum (P7). n = 7 for WT and KO, respectively. Scale bar, 100 μm. c–f Representative images from immunostaining analysis with antibodies against Ki67 (c), phospho-H3 (PH3) (d), BrdU (e), and NeuN (f) in the cerebellum of WT and KO mice. Quantification of immuno-reactive cells is shown in the right panels. Numbers of biological replicates are included in parentheses. Black lines in the graph indicate the mean. *P < 0.05, **P < 0.01, ns not significant. Sections in c and e were stained with DAPI and sections in d and f were counterstained with eosin to visualize nuclei. Scale bar, 100 μm. g, h Representative images from immunohistochemical analysis with antibodies against Calbindin-D28K (g) and GFAP (h) to detect Purkinje cells (g) and astrocytes (h). n = 7 for WT and KO, respectively. Scale bar, 100 μm. *P < 0.05, **P < 0.01, ***P < 0.001, ns not significant and PH3 immunoreactivity (Additional file 1:FigureS7d, deficiency demonstrates that proper orchestration of the e). In addition, the cerebellum in Alkbh5 KO mice spatiotemporal expression of m A writers and erasers is remained smaller than that of the littermate controls necessary for mouse cerebellar development. (Additional file 1: Figure S7f). These data indicate that in spite of the undefined compensatory response to partially Imbalance of RNA m A methylation in the Alkbh5-deficient recover the morphology changes, Alkbh5 deficiency had a mouse cerebellum exposed to hypobaric hypoxia profound and deleterious effect on cerebellar development To examine how Alkbh5 deficiency affects m A under hypoxic conditions. methylation of the cerebellar poly(A) RNA, we first Together, the abnormal cerebellar development resulting performed UHPLC-MS/MS analysis on the poly(A) from either ectopic expression of Mettl3 or Alkbh5 RNA from cerebella at P7 and found a slight increase in Ma et al. Genome Biology (2018) 19:68 Page 9 of 18 6 6 the global m A levels in the KO mice exposed to increase of m A at the global level, Alkbh5 deficiency hypobaric hypoxia (Additional file 1: Figure S8a). To led to disordered m A levels of thousands of RNAs. gain more detailed insight into the epitranscriptome- We next performed GO analysis to evaluate whether wide changes of all methylated RNAs, we performed the defective cerebellar development in KO mice re- 6 6 m A-seq using the poly(A) RNAs from both WT and sulted from imbalanced m A methylation (Fig. 6c, d and KO mice exposed to hypobaric hypoxia (Additional file Additional file 1: Figure S8e). Partial or complete loss of 1: Table S1). A slight increase in the number of m A methylation mainly occurred to the genes participating peaks and methylated RNAs was observed in the in the control of cell division (such as Cenpe in cerebella of KO mice (Additional file 1: Table S2), as Additional file 1: Figure S8c), cell cycle (such as Cdca2 well as higher enrichment scores of m A peaks (Fig. 6a). in Additional file 1: Figure S8d), and cell projection In addition, we identified 1348 poly(A) RNAs with gain organization (such as Erbb4 in Additional file 1: Figure of methylation and 711 poly(A) RNAs with loss of S8d). In contrast, the hypermethylated RNAs were methylation in KO mice cerebellum (P7) (Fig. 6b; found to be related to metabolic processes (such as Ccr5 Additional file 1: Figure S8b, c). By comparing those in Additional file 1: Figure S8c), ion transport (such as abnormally methylated RNAs to the SMRs during Camk2g in Additional file 1: Figure S8d), and axon cerebellar development of the WT mice, we found that guidance (such as Rora in Additional file 1: Figure S8d). upon Alkbh5 deficiency, 53 original P7 SMRs lost their Previous studies showed that ALKBH5 regulates RNA nu- methylation, while 42 original P14 SMRs, 53 P21 SMRs, clear export via an m A-dependent manner. We therefore and 104 P60 SMRs became methylated at P7. In asked whether the changes in m A methylation resulting addition to the gain- and loss-of-methylation RNAs, we from Alkbh5 deficiency also affect the RNA nuclear export also identified a considerable number of RNAs exhibit- in vivo. Based on functional enrichment analysis, we chose ing a significant increase (514) or decrease (81) in their several genes in different functional pathways, such as cell methylation levels (Fig. 6b). Therefore, in spite of a mild division (Mphosph9)[40], membrane potential (Opa1)[41], Fig. 6 Imbalanced m A RNA methylation in the cerebellum of Alkbh5-deficient mouse exposed to hypobaric hypoxia. a Boxplots showing the relative methylation levels of all m A peaks between wild-type (WT)and Alkbh5-deficient (KO) mouse cerebellum. ***P < 0.001 by Wilcoxon test. b Venn dia- gram showing the relationship among the RNAs containing gain-of-methylation (Gain), lost-of-methylation (Loss), and differentially methylated RNAs (DMRs). c, d GO analysis of genes encoded by RNAs containing loss- (c)or gain-of-methylation (d) in KO mouse cerebellum. BP biological process, CC cellular component. e Gene-specific m A-IP-qPCR results showing the relative methylation levels of four RNAs in the cerebellum of KO mice compared to those in WT mice. n =7. f Subcellular localization of the four RNAs as shown by the ratio of RNA abundance between the cytoplasmic and nuclear fractions isolated from the cerebellum of WT and KO mice. n = 10 for WT and 8 or 10 for KO, respectively. *P < 0.05, **P <0.01, ns not significant Ma et al. Genome Biology (2018) 19:68 Page 10 of 18 microtubule cytoskeleton organization (Wdpcp)[42], and improved our understanding of the mechanism related to ion transport (Letm1)[43], for further investigation. We first mouse cerebellar development [34, 35, 44]. Given the role 6 6 conducted gene-specific m A qPCR assays to validate the of m A in fate determination of embryonic stem cells and increased m A levels of these four mRNAs in the mammalian cortical neurogenesis [16, 17, 31, 45], we cerebellum of KO mice (Fig. 6e). Subsequently, subcellular hypothesized that m A also participates in neuronal cell fractionation analysis of cerebellar cells was performed proliferation and differentiation in the mouse cerebellum. (Additional file 1: Figure S8f, g). Compared to WT mice In this study, we show that a considerable number of cerebellum, all the four RNAs invariably exhibited increased RNAs in the mouse cerebellum exhibit developmentally abundance in the cytoplasm in the Alkbh5 deficient regulated methylation. GO analysis implied that cerebellum (Fig. 6f). Together, these data indicate that numerous genes regulating cell proliferation, Alkbh5 deficiency led to dysregulated RNA nuclear export differentiation, or metabolism acquire m A marks in a in the cerebellum. temporal-specific manner, among which the expression of Given the role of m A in regulating RNA stability, we 839 RNAs was under tight regulation by their methylation also analyzed the gene expression changes resulting from status. Notably, we found that the number of differentially Alkbh5 deficiency. Among the 497 differentially expressed methylated RNAs is much higher than that of differen- RNAs, 81 exhibited changed methylation and were tially expressed RNAs as observed during cerebellar devel- enriched in functional pathways, including transport, cell opment and in Alkbh5 KO mice. We speculate that this is division, and cell cycle control (Additional file 6:Table probably due to the additional functions of RNA methyla- S8). For instance, the genes involved in cell cycle control, tion in RNA splicing, nuclear export, and protein transla- such as Ddx11, Ccnb1,and Cbx1, displayed a reduction in tion [5]. For example, although the hypermethylated their RNA expression together with their reduced RNAs presented in Fig. 6e did not change in their expres- methylation. For genes associated with transport or sion levels, their nuclear export was enhanced by Alkbh5 metabolism, Alkbh5 deficiency-induced RNA hypermethy- deficiency. Since it is the transcripts in the cytoplasmic lation led to either a decrease (such as Cacna2d3, Notch3, pool that are utilized for RNA translation to protein and and Jam3) or increase (such as Slc6a20a, Dlk1, Slc2a4, RNA degradation, appropriate m A levels are necessary and Ascl2) in their RNA expression in the cerebellum for maintaining RNA nuclear export in a balanced state. (Additional file 1: Figure S8h). These results imply that the selective m A marks, through Taken together, our data confirm that, by altering the their various modulatory effects on RNA processing, original m A level, Alkbh5 deficiency disturbed the RNA function as epitranscriptomic switches that either activate metabolism of a subset of cell fate determination genes or suppress a series of physiological events during mouse and thus led to defective cerebellar development in the cerebellar development. In line with our findings in the mice exposed to hypobaric hypoxia. mouse cerebellum, dynamic changes in m ARNA methylation were also observed during postnatal liver Discussion development in pig [46] and mouse spermatogenesis [22]. Brain development is known to be regulated in a precise The dynamic RNA methylation levels in the developing spatiotemporal manner. We have recently identified cerebellum reveal how important it is to maintain the bal- region-specific m A with distinct methylation patterns ance of RNA methylation in the developing mouse cerebel- and biological functions for the mouse cerebellum and lum, which is accomplished by both the m A cerebral cortex [36]. In the present study, we further methyltransferases and demethylases. The morphological demonstrate that developmentally regulated m A changes observed upon transient knockdown or methylation is one of the important regulatory factors overexpression confirm that METTL3 is involved in for postnatal mouse cerebellar development. Temporal- regulating cerebellar development. In addition to specific methylation differs from continuous methylation methyltransferase, RNA m A demethylase is another in its methylation level, biological functions, and important factor to maintain RNA methylation at dynamic distribution of m A marks. Furthermore, we appropriate levels. The important role of FTO in the uncovered that abnormal expression of m A mouse cerebellum was shown in a previous report using a methyltransferases (METTL3) and demethylases Fto-knockout mouse model, which is characterized by a (ALKBH5) caused imbalanced RNA methylation, which cerebellum that is significantly smaller than that in WT further led to defective cerebellar development. We also litter mates [28]. ALKBH5, another demethylase identified discovered that Alkbh5 deficiency led to facilitated earlier, regulates spermatogenesis [20] and maintains nuclear export of those hypermethylated RNAs. tumorigenicity of glioblastoma stem-like cells [47]. How- Previous identified changes in genome-wide gene ex- ever, its function in the brain remains unclear. Here we pression and epigenetic modifications throughout the only observed defects in cerebellar development in Alkbh5- whole cerebellar developmental processes have greatly deficienct mice upon exposure to hypobaric hypoxia. It is Ma et al. Genome Biology (2018) 19:68 Page 11 of 18 likely that ALKBH5, FTO, and probably other undefined developmental process. Our results show that the positive demethylases overlap in their demethylation activity and or negative regulation of m A marks on these cell can partially substitute for each other’s functions under fate-determining genes is important to drive the proper normal conditions [37]. On the other side, our results also progression of postnatal cerebellar development. reveal the distinct functions between ALKBH5 and FTO. Different from previous in vitro analysis showing that In line with a previous report [9], we provide in vivo evi- m A marks are predominantly located near stop codons dence that ALKBH5 regulates the nuclear export of the and 3′ UTRs, previous in vivo methylation analysis of hypermethylated RNAs, while FTO regulates RNA splicing mouse cerebellum and cerebral cortex identified the in a methylation-dependent manner [7]. In this study, we m A marks in start codon regions as well as stop codons observed increased cell proliferation in the Alkbh5-defi- and 3′ UTRs [36]. Here, we identified similar patterns of cienct mouse cerebellum exposed to hypobaric hypoxia. In m A distribution throughout the process of postnatal contrast, Fto deficiency in adult mouse hippocampus led to cerebellar development. Importantly, we observed a reduced cell proliferation in adult neural stem cells [28]. significant change in the distribution of temporal- The discrepancy in these results confirmed that RNA specific m A marks from P7 to P60. We deduce that the methylation constitutes a complex process and suggests distribution patterns of temporal-specific m A peaks that the regulation and function of m A methylation in the possibly reflect the difference in methylation sites mouse brain should be investigated in a more precise between proliferating, immature and fully differentiated, spatiotemporal-dependent manner. mature neural cells. In the mouse cerebellum of P7, In the WT mouse cerebellum at P7, the m A during which the majority of neural cells are actively modification selectively marked a group of RNAs proliferating, the temporal-specific m A peaks were responsible for cell cycle control and cell division. In enriched in CDS, stop codons, and 3′ UTRs, which was Alkbh5 KO mice exposed to hypobaric hypoxia, we similar to the results obtained from the m A-seq detected many of these types of RNAs with aberrant analysis performed using cell lines [2, 3, 48]. In contrast, methylation, which may further disturb RNA metabolism, distribution of temporal-specific m A peaks in start such as RNA splicing, nuclear export, RNA decay, or codon regions was most significant in mouse cerebellum translation. In agreement with these findings, we observed at P60, which mainly consists of well-differentiated and an increase in the numbers of mitotic cells in the mature neural cells. In support of our findings, m A cerebellum of Alkbh5 KO mice upon exposure to marks in start codon regions were also observed in the hypobaric hypoxia. Mechanistically, we showed the fully differentiated tissues from mouse cerebral cortex hypermethylation and enhanced nuclear export of [36] (unpublished data), A. thaliana [49], rice [50], liver Mphosph9, which encodes a centrosomal protein related to [46], and muscles [51]. Collectively, these data imply cell division and cell cycle control [40]. Therefore, that the methylation in start codon regions occurs not abnormal RNA methylation is likely one of the reasons only in mature neural cells of the mouse cerebellum, but resulting in dysregulated cell proliferation and cell cycle also in many other types of fully differentiated cells in progression observed in the Alkbh5-deficient cerebellum. diverse organisms. We thus propose that the position of 6 6 Consistent with our findings, deficiency in the m A reader m A marks included in the transcript may be an Ythdf2 delayed DNA synthesis and mitotic processes in a important factor that determines the functions mediated zebrafish model [24]. In another study, deletion of Mettl14 by m A. Further investigation will be necessary to clarify in embryonic mouse brain also led to prolonged cell cycle the precise roles of m A in different physiological of cortical neural progenitors [31]. Together, these processes. However, the m A-IP-seq used in this study 6 6 phenotypes provide strong evidence that balanced m A was unable to identify the exact position of the m A levels are essential for cell cycle control. Moreover, we marks. We cannot rule out the possibility that the m A found that a higher number of RNAs were marks located in 5′ UTR may be misclassified into start hypermethylated in Alkbh5-deficient cerebella of P7, which codon regions for those poly(A) RNAs with a very short were originally methylated only in fully differentiated, 5′ UTR. Thereby, single nucleotide-resolution m A mature neural cells in WT cerebellum. Such untimely and analysis [48, 52] should be utilized next to identify the aberrant methylation finally appears to result in premature exact positions of m A marks, which will further enable differentiation and aberrant neurogenesis in the exploration of the functions and regulatory mechanism cerebellum. This may explain the phenotypes of lower of m A deposited in different regions of the transcripts. numbers of mature granule neurons in IGLs, together with In addition, the different patterns of m A observed disorganized Purkinje cells, and astrocytes. Generally, cell between cell lines and tissues indicates the complexity of fate determination requires a burst of expression of a in vivo temporal regulation of RNA methylation, distinct group of genes at each specified stage, which is suggesting the necessity to characterize RNA spatially and temporally regulated throughout the whole methylation in a context-dependent manner. Ma et al. Genome Biology (2018) 19:68 Page 12 of 18 Gene regulation in the brain is characterized by its con- Antibodies siderable cellular heterogeneity, including gene expression Detailed information on the primary and secondary anti- [53], DNA hydroxymethylation [34], and m ARNA bodies used in this study is listed in Additional file 1: methylation [36]. Accordingly, we found that during Table S9. cerebellar development, the expression of m Awriter and eraser genes in granule neurons and Purkinje neurons RNA isolation, poly(A) RNA purification, and m A-IP changed in opposite directions. Although the poly(A) RNAs Mice were sacrificed by cervical dislocation and cerebel- analyzed in this study were isolated from a mixture of all lum was removed and immediately snap-frozen in liquid neural cell types in the mouse cerebellum, the methylation nitrogen. Total RNA was extracted using TriReagent profiles depicted here mainly reflect the nature of m A (Sigma, T9424) according to the manufacturer’s instruc- modification in granule neurons as they constitute about tions. In general, four mice of P7, four mice of P14, and 90% of the cerebellar neural cells [35]. For further insights, three mice of P21 or P60 were used for one m A-IP it will be important to elucidate the methylation profiles reaction according to the established protocol [54]. and related biological functions in Purkinje cells and glial Briefly, poly(A) RNA was purified from total RNA using cells, respectively. A more detailed, cell type-specific methy- Poly(A)Purist™ MAG Kit (Thermo Fisher, AM1922). lation analysis will be necessary to gain insights into the bio- Next, poly(A) RNA was fragmented to 100-150 nucleotides 6 6 logical functions and developmental regulation of m A (for m A-seq) using RNA Fragmentation Reagents modification in individual types of neural cells. (Thermo Fisher, AM8740) as instructed by the manufacture. Fragmented poly(A) RNA (1 μg) was Conclusions incubated with 2.5 μgofanti-N -methyl-adenosine Here, we describe the first m A RNA methylation antibody at 4 °C for 4 h, followed by addition of landscape during postnatal development of the mouse Dynabeads™ Protein A for Immunoprecipitation (Thermo cerebellum and reveal the unique features of continuously Fisher, 10002D). After incubation at 4 °C for 2 h, beads and temporally specific methylation in the mouse were extensively washed five times with IPP buffer cerebellum from P7 to P60. Our results indicate that (150 mM NaCl, 0.1% NP-40, 10 mM Tris-HCl, pH 7.4); temporal regulation of RNA m A methylation, orchestrated immunoprecipitated RNA was recovered through elution 6 6 by spatiotemporal expression of m A writers and erasers, is with m A nucleotide followed by ethanol precipitation. essential in regulation of postnatal cerebellar development. For m A-seq, a cDNA library was constructed using a We also provide the first in vivo evidence that ALKBH5- KAPA Stranded mRNA-seq kit (Kapa Biosystems, induced abnormal RNA methylation affects RNA nuclear KR0960) and subjected to next-generation sequencing on export. Further studies will be necessary to uncover the de- an Illumina Hiseq X Ten system. Input RNA before im- tailed molecular mechanism of the regulation and bio- munoprecipitation was applied to RNA-seq in parallel. 6 6 6 logical functions of m A during mouse cerebellar For gene-specific m A-IP-qPCR, m A-IP was development. In addition, characterization of the role of performed using 3 μg of fragmented poly(A) RNA (300– m A in the cerebellum will offer a new viewpoint to 500 nucleotides). Reverse transcription was carried out elucidate the mechanism of related neurological diseases. with an equal ratio of RNA from input and IP product by using ReverTra Ace® qPCR RT Master Mix Methods (TOYOBO, FSQ-301). Quantitative real-time PCR was Animal care performed by using THUNDERBIRD™ SYBR® qPCR Mix Wild-type C57BL/6 mice were purchased from Vital River (TOYOBO, QPS-201). Gapdh was used as a negative company (Beijing, China). The Alkbh5-deficient mice were control. The primers used in this study are listed in used and genotyped as described previously [9]. For all ex- Additional file 1: Table S10. periments performed here, equal numbers of male and fe- male mice were included for analysis. All animal Protein extraction and western blot analysis experiments and euthanasia were approved and performed Cerebella were dissected from mice of different ages and in accordance with the guidelines of Animal Care and Use triturated in RIPA lysis buffer with protease inhibitor Committee of IBMS/PUMC. The IRB (Institutional Review cocktail (Roche, 04693159001) and phosphatase inhibi- Board) approval number is ACUC-A02–2014-001. tors. Tissue lysate (60–80 μg) was subjected to 10% SDS-PAGE. After electrotransfer, the blots were blocked Cell lines in 5% milk in Tris-buffered saline with Tween 20 (TBST) Neuro2a cells and HEK293T cells were purchased from for 1 h at room temperature, then incubated at 4 °C National Infrastructure of Cell Line Resource (Beijing, overnight with the primary antibodies. The transferred China). The cell lines have been authenticated and tested membranes were incubated with appropriate HRP- for mycoplasma contamination by the provider. conjugated secondary antibodies, exposed to enhanced Ma et al. Genome Biology (2018) 19:68 Page 13 of 18 chemiluminescence (ECL) solution, and visualized by gel VSV-G (kindly provided by Prof. Qi Xu). shRNA of image analysis system (Tanon, 5800). mouse Mettl3 was cloned into pLL3.7 vector between HpaI and XhoI restriction enzyme sites. Two sets of Histology and immunostaining analysis shRNAs were used for each gene and were evaluated for After anesthesia with saturated tribromoethanol solu- their knockdown efficiency in Neuro2a cells via western tion, mice were perfused with 4% paraformaldehyde; the blot analysis. The shRNA with higher efficiency was whole brains were dissected and post-fixed in 4% para- chosen for subsequent in vivo virus infection of mouse formaldehyde for 24–48 h at 4 °C. For paraffin sections, cerebellum. The oligos of Mettl3 shRNA were: 1, 5′-aacc the cerebellum was dehydrated with ethanol, clarified GCACACTGATGAATCTTTAGGTtcaagag ACCTAAA- with xylene, and then embedded in paraffin. We used 4- GATTCATCAGTGTGCtttttc-3′;2,5′-aaccCTGCAAA- μm-thick sections for hematoxylin and eosin (H&E) TATGTTCACTATGAtcaagag TCATAGTGAACATATT staining or immunohistochemistry analysis. For cryosec- TGCAGtttttc-3′. The negative control containing tions, the cerebella were dehydrated in 30% sucrose for scrambled shRNA was 5′-aaccTTCTCCGAACGTGT- more than 24 h, embedded in OCT compound and fro- CACGTtcaagagACGTGACACGTTCGGAGAAtttttc-3′. zen at − 40 °C. We used 12-μm-thick sections for im- The plasmids generated here were written in short as munofluorescence analysis. Immunohistochemistry and OE, sh-1, sh-2. immunofluorescence staining were performed according to standard protocols [55]. Lentivirus infection of mouse cerebellum For BrdU incorporation analysis, the mice were intra- All lentivirus vectors were purified and amplified using peritoneally injected with BrdU (Sigma, B5002; final QIAGEN Plasmid maxi kit (Qiagen, number 12163). 293 T concentration 50 μg/g) and analyzed 2 hours after BrdU cells were co-transfected with a mixture of lentiviral plas- injection. Immunofluorescence analysis was performed mids (including 17.5 μg of pLV-derived plasmid, 13.125 μg as described above by using anti-BrdU antibody. of pH 1, 4.375 μgofpH2foroverexpressionsystem; 14 μg Images of H&E and immunohistochemical staining were of pLL3.7-derived plasmid, 7 μg of pMDLg/pRRE, 7 μgof recorded using a Panoramic MIDI II digital slide scanner pCMV-VSV-G, 7 μg of pRSV-Rev for knockdown system) (3D Histech). Results of immunofluorescence analysis using DNA transfection reagent (Neofect, TF20121201). were imaged with a confocal microscope (Zeiss, LSM780). Ssupernatants were collected 48 h after transfection and concentrated by ultracentrifugation. Titration of lentivirus Quantitative m A level measurement using UHPLC-MS/MS was carried out by transducing 293 T cells in tenfold serial The absolute amount of m A in the poly(A) RNA was dilution. Generally the lentivirususedinthisstudy hada 8 10 measured as described before [36]. Briefly, 100 ng of titer of 1 × 10 to 1 × 10 TU/ml. poly(A) RNA was digested with nuclease P1 (Sigma) at In vivo lentivirus infection was performed as previ- 37 °C for 2 h, followed by treatment with calf intestine ously described with minor modification [56]. Briefly, alkaline phosphatase (CIAP, Promega) at 37 °C for 2 h. wild-type C57BL/6 mice (P7) were anesthetized on ice The solution was filtered and analyzed using ultra- and then injected with a mixture of lentivirus (2 μl), performance liquid chromatography and a triple-4 quad- polybrene (800 ng/μl), and fastgreen (0.1%) into the rupole mass spectrometer (AB SCIEX QTRAP 5500). interlobular space between lobules V and VI. Pups were m A level in each sample was calculated by comparing revived at 37 °C and then returned to original cages. with the standard curve obtained from pure nucleoside After 7 (for overexpression) or 10 (for knockdown) days standards loaded simultaneously. The ratio of m AtoA post-lentivirus infection, immunofluorescence analysis was calculated to obtain the global methylation level. was performed as described above to detect the pheno- types of cerebellum upon lentivirus infection. Plasmid construction The lentiviral overexpression vectors used in this study Hypobaric hypoxia treatment of mice include pLV-EF1a-EGFP(2A) Puro, pH 1, and pH 2 vec- WT or Alkbh5 KO mice of P5, together with their mothers, tors (kindly provided by Prof. Kaili Ma). cDNA of mouse were placed into a hypobaric hypoxic chamber and exposed Mettl3 was cloned into pLV-EF1a-EGFP(2A) Puro vector to a simulated atmospheric pressure of 10.6 kPa digested with EcoR I and BamH I restriction enzymes. (377 mmHg), which was equivalent to an altitude of The overexpression efficiency of lentivirus was verified 5500 m. All animals were kept at constant temperature in Neuro2a cells by western blotting. The primers used (25–30 °C) on a daily light schedule of 12 h of light vs. dark were: forward, 5′-cggaattc ATGTCGGACACGTGGAG- with normal activity. The hypobaric hypoxic condition was 3′; reverse, 5′- cgggatccCGCTATAAATTCTTAGG-3′. maintained and monitored continuously with the sensor in- −/− The lentiviral knockdown system used in this study in- side the chamber. Wild-type and Alkbh5 pups were cludes pLL3.7, pRSV-rev, pMDLg/pRRE, and pCMV- dissected for phenotype analysis after incubation in the Ma et al. Genome Biology (2018) 19:68 Page 14 of 18 hypobaric chamber for 48 h. To examine the phenotypes of two mismatched bases; 2) bases with low quality score hypoxia-treated mice at later developmental stages, the WT (< 20) were trimmed off from the 3′ end; 3) reads with and KO mice were returned to normoxic conditions to- length longer than 70 nucleotides and more than 70% gether with their mothers and raised for another 7 days be- bases with quality score > 25 were retained. These high- fore dissection for phenotype analysis. quality reads were mapped against the mouse genome (mm10) allowing up to two mismatches using TopHat Magnetic resonance imaging analysis software (version 2.0.13) [57, 58]. Only uniquely mapped Mice were anaesthetized with 0.75% amobarbital before reads were kept for the subsequent analysis. analysis. Mouse brain MRI was performed on a 7.0 T BioClinscan Animal MRI System (Bruker, Germany) with RNA expression analysis Siemens manipulating software. T2-weighted images of The FPKM (fragments per kilobase of transcript per million the whole brain (including the sagittal, coronal, and trans- mapped reads) values of RNAs in each sample were calcu- verse views) were acquired with a 2D-TSE (2D-turbo lated using Cufflinks toolkit (version 2.0.2) [59]. Only the spin-echo) sequence. The scanning time for each mouse transcripts with FPKM value > 0.2 were considered as was about 10 min. The scanning parameters were as fol- expressed transcripts [60]. The differentially expressed lows: sagittal view, TR = 2360 ms, TE = 41 ms, number of RNAs (DERs, P < 0.05) between any two samples were averages = 1, 15 axial slices with a slice thickness (ST) of 0. identified using Cuffdiff. The normalized FPKM value for 7 mm, a field of view (FOV) of 3.52 cm × 3.52 cm, and a each sample was used for clustering using the KMC matrix of 320 × 320 were positioned over the whole brain method of MeV software with the following parameter set- with a pixel spacing of 0.11 mm; coronal view, TR = tings: 1) Pearson correlation clustering; 2) K-means cluster- 3510 ms, TE = 53 ms, number of averages = 1, axial slices ing for six clusters; 3) other parameters set as defaults [61]. = 19, ST = 0.5 mm, image matrix size = 320 × 384, pixel To verify the reliability of the DERs, we performed a par- spacing = 0.078 mm, FOV 2.50 cm × 3.00 cm; transverse allel differential expression analysis by STAR/edgeR pack- view, TR = 3380 ms, TE = 41 ms, number of averages = 1, ages. High quality reads in each sample were mapped to axial slices = 21, ST = 0.7 mm, image matrix size = 320 × the mouse genome (mm10) using STAR software and the 320, pixel spacing = 0.094 mm, FOV 2.82 cm × 2.82 cm. uniquely mapped reads were kept for the subsequent ana- Then the signals of different brain regions were measured lysis [62, 63]. HTSeq software [64] was utilized to calculate and compared using a RadiAnt DICOM Viewer. the read count of RNAs and edgeR [65] was applied to identify the DERs (P < 0.05). All software applications were Subcellular fractionation analysis of RNA run with default parameters. Fresh cerebella of WT and KO mice at P7 were minced and incubated in Digestion Solution (30 U/mL papain, 240 μg/ m A peak calling and motif analysis 6 6 mL cystein, 400 μg/mL DNAseItype IV)at37°Cfor 1h. RNA m A-modified regions, also called m A peaks, were The reaction was stopped by adding Ovomuccoid inhibitor identified using exomePeak software (version 2.7.0) with solution (1125 μg/mL Ovomuccoid trypsine inhibitor, FDR (false discovery rate) < 0.05 [66, 67]. The consensus 525 μg/mL BSA, 400 μg/mL DNase I type IV) and incuba- sequence motifs enriched in m A peaks were identified by tion at 37 °C for 4 min. After purification, cerebellar cells HOMER [68]. For further comparison of m A modification were subjected to cellular fractionation using a PARIS™ Kit between samples, we used coverageBed of BEDToods (Thermo Fisher, AM1921) according to the manufacturer’s (version 2.26.0) [69]withthe “-s”, “-splited”, “-counts”,and protocol. Nuclear and cytoplasmic RNAs (120 ng) were re- “-F 0.50” parameters to calculate the read count of each verse transcribed using ReverTra Ace® qPCR RT Master peak. The “IP FPKM”, “input FPKM”,and “Enrichment Mix (TOYOBO, FSQ-301). Quantitative real-time PCR was score” of peaks were calculated as following methods: performed using THUNDERBIRD™ SYBR® qPCR Mix a ¼ A x10 = B xC ð1Þ i; j i; j j i; j (TOYOBO, QPS-201). Nd1 was used as a cytoplasmic marker while U1 was used as a nuclear marker. The primers b ¼ D x10 = E xC ð2Þ i; j i; j j i; j used in this study are listed in Additional file 1:Table S10. c ¼ a =b ð3Þ i; j i; j i; j RNA-seq data processing and reads mapped For each sample, single-end reads were used for bio- where a denotes the peak “IP FPKM” value of the ith i,j informatics analysis. Quality control of raw data was methylation peak in the IP sample from the jth bio- done using FastQC software (version 0.10.1). Sequencing logical sample; A denotes the total reads mapped to the i,j data were preprocessed with in-house Perl scripts using ith methylation peak in the IP sample from the jth bio- the following criteria: 1) adaptor sequence was removed logical sample; B denotes the total unique reads mapped by finding the sequence AGATCGGAAG with at most to the mouse reference (mm10) in the IP sample from Ma et al. Genome Biology (2018) 19:68 Page 15 of 18 the jth biological sample; C denotes the length (base) bin was calculated to represent the occupancy of m A i,j of the ith methylation peak in the IP sample from the peaks along the whole transcripts. jth biological sample; b denotes the peak “input FPKM” We further counted the number of m A peaks in the i,j value of the ith methylation peak in the input sample 5′ UTR, start codon region, CDS, stop codon region and from the jth biological sample; D denotes the total 3′ UTR, respectively. Specifically, a 300-nucleotide i,j reads mapped to the ith methylation peak in the input region centered on start codons or stop codons was de- sample from the jth biological sample; E denotes the fined as start codon regions or stop codon regions [36]. total unique reads mapped to the mouse reference (mm10) in the input sample from the jth biological sam- 6 Correlation analysis between RNA m A methylation level ple; c denotes the peak “enrichment score” value of the i,j and RNA expression level ith methylation peak in the IP/input sample from the jth The global clustering analysis of RNA methylation levels biological sample. (enrichment score) and expression levels (FPKM) among M A peaks that satisfied 1) exomePeak FDR value < 0. the four stages was performed using Cluster 3.0 soft- 6 6 05, 2) m A peak “IP FPKM” value > 1, and 3) m A peak ware. The enrichment scores and FPKM of each RNA in “enrichment score” value > 1.5 were used for further each sample were adjusted with “mean” parameters in comparative analysis. The continuously methylated “Normalize genes” and “Center genes”. The hierarchical RNAs (CMRs) were defined as RNAs containing at least clustering method was employed using “Complete link- one m A peak in all the four samples, while the age”. The clustered heatmap figure was generated using temporal-specific methylated RNAs (SMRs) were the TreeView-1.1.6r4-win software. RNAs with m A modification only in one sample, and In order to obtain a dynamic view of the regulatory effect 6 6 m A peaks of SMRs were defined as specific m A peaks. 6 of RNA m A methylation on RNA abundance, we By merging all the m A peaks of four developmental calculated the correlation coefficients and the associated P stages using mergeBed of BEDTools with “-s”, “-c”, and 6 valuebetween theenrichment scoresofm A peaks and “-o” parameters, we identified “ON” or “OFF” RNA m A FPKM of methylated RNAs at the four stages using switches during cerebellar development. All m A peaks Pearson test in R package. High-confidence regulatory gene absent in the former stage but present in the later stage pairs were defined by the correlation coefficients (|r| > 0.95) 6 6 were called “ON” m A switches, while all m A peaks with P value < 0.05. Cluster analysis was performed for displaying changes in the opposite direction were called those positively (r > 0.95) and negatively (r < − 0.95) corre- “OFF” m A switches. In addition, the differentially lated RNAs using MeV software, respectively. The MeV methylated RNAs (DMRs) between wild-type and cluster results were integrated and plotted using R pro- Alkbh5 knockout mice were further identified using exo- gramming language. mePeak software with default parameters, and high- confidence DMR m A peaks matching the criteria of GO analysis FDR < 0.05 and log2 (fold change) > 1 or < − 1 were The DAVID tool [74] was used for GO analysis by apply- selected for further analysis. ing default parameters, except that only those transcripts For validation analysis of the distribution patterns of 6 fulfilling the condition of FPKM > 0.2 in the input samples stage-specific m A peaks which were obtained by using were set as background. Either bubble plots or column the exomePeak software, peak calling analysis was plots were generated based on the enriched GO terms performed in parallel using the MACS2 software for the using the ggplot2 in R package [75]. Color intensity indi- P7 and P60 samples with the default parameters [54, 70]. cates the value of −log (P value), while the length of col- The integrative genomics viewer (IGV) tool was used 6 umns or the size of the circles indicate the gene counts. A for visualization of m A peaks along the whole transcript 6 full list of all selected terms of the biological process, [71, 72]. The heatmaps of the enrichment score for m A cellular components, and molecular functions category peaks were plotted using the pheatmap package in R [73]. are provided in Supplemental tables (Additional file 2-4: TableS4-S6; Additionalfile 6: Table S8). Characterization of m A peak distribution patterns Distribution of m A peaks along mRNAs were Kyoto Encyclopedia of Genes and Genomes analysis determined as previously described with minor revision All modified genes were annotated to Kyoto Encyclopedia [3, 18]. To characterize the distribution patterns of m A of Genes and Genomes (KEGG) pathways using KAAS peaks, a reference mouse transcriptome was built using (KEGG Automatic Annotation Server) [76], and the the longest transcript of each gene. Next, each of the 5′ enriched KEGG pathways for the genes encoded by P7 or UTR, CDS, and 3′ UTR regions were split into 100 bins P60 temporal-specific SMRs or CMRs were identified with equal length. The percentage of m A peaks in each using Fisher’s exact test (Additional file 4: Table S6) [77]. Ma et al. Genome Biology (2018) 19:68 Page 16 of 18 6 6 Quantification and statistical analysis Genomes; KO: Alkbh5-knockout; m A: N -methyladenosine; P7: Postnatal day 7; PCL: Purkinje cell layer; PH3: Phospho-histone 3; SMRs: Specifically The differences in distribution for m A enrichment score methylated RNAs; UHPLC-MS/MS: Ultra high performance liquid and log FPKM between samples were detected by chromatography-mass spectrum/mass spectrum analysis; TBST: Tris Buffered Wilcoxon test. Pearson test was used to perform Saline with Tween 20; WT: Wild-type correlation analysis. All statistical analysis and graphs of Acknowledgements results in Figs. 3d, 4d, 5a, 5c–5f and 6f, and Additional file We thank Prof. Xu Q and Dr. Ma K-L for their generosity in offering lentivirus 1: Figures S5c, S6a, S7d, S7e, S8a and S8g were assessed vectors. using two-tailed unpaired Student’s t-test and performed Funding using GraphPad Prism 6.0 software. Results are presented This work was supported by the National Natural Science Foundation of China as mean ± s.e.m. Immunoreactive cells in 3–5 randomly se- (31471288 to NY, 31471343 to TWM), CAMS Initiative for Innovative Medicine (2016-I2M-1-004 to NY, 2016-I2M-2-001 to TWM), and The Youth Innovation Pro- lected lobules of one cerebellum were counted using Stata- motion Association of Chinese Academy of Science (2017141 to SS.) Quest 5.0 software. For each subject, data were collected from five to seven mice. For each mouse, data were ob- Availability of data and materials The datasets have been deposited in the Sequenced Read Archive (SRA) tained from the immunostaining results of two to three under accession number PRJNA400297 (SRX3136161-SRX3136172) [78], and near-midline slices. also in the Genome Sequence Archive [79] in BIG Data Center [80], Beijing Institute of Genomics (BIG), Chinese Academy of Science under accession number CRA000472 [81]. Additional files Authors’ contributions NY and TWM conceived the project and, with SS and KA, designed the data Additional file 1: Figure S1. Dynamic RNA methylation in the analysis. MC performed phenotype analysis. CM performed m A-IP and developing mouse cerebellum. Figure S2. Comparison between phenotype analysis. LH, ZW and TX performed bioinformatic analyses. ZZW continuous and temporal-specific methylation during mouse cerebellar performed lentivirus infection experiments. HX performed mouse brain development. Figure S3. Comparison of RNA methylation between P7 dissection experiments. WG performed plasmid construction experiments. ZS and P60 based on m A peaks identified using MACS2. Figure S4. and ZY assisted with immunostaining analyses. WD assisted with western Correlation between RNA methylation and gene expression during blot analyses. LC and LQ performed MRI analysis. NY, MC and LH wrote the cerebellar development. Figure S5. Morphology analysis of mouse manuscript with input from all authors. All authors read and approved the cerebellum upon lentivirus infection for Mettl3 overexpression. Figure S6. final manuscript. Phenotype analysis of the cerebellum in the WT and KO mice under normoxic condition. Figure S7. Morphology analysis of the cerebellum in Ethics approval WT and KO mice exposed to hypobaric hypoxia and normoxia All animal experiments and euthanasia were approved and performed in successively. Figure S8. Dysregulated RNA methylation resulting from accordance with the guidelines of Animal Care and Use Committee of IBMS/ Alkbh5 deficiency in mouse cerebellum exposed to hypobaric hypoxia. PUMC. The IRB (Institutional Review Board) approval number is ACUC-A02– Table S1. Data quality and processing information of m A-seq of poly(A) 2014-001. RNA from wild-type mouse cerebellum (P7, P14, P21, and P60), the cere- bellum of wild-type (WT) and Alkbh5 knockout (KO) mice exposed to Competing interests hypobaric hypoxia (P7). Table S2. Statistics of m A peaks and expressed The authors declare that they have no competing interests. RNAs in wild-type mouse cerebellum (P7, P14, P21, and P60), the cerebel- lum of wild-type (WT) and Alkbh5 knockout (KO) mice exposed to hypo- baric hypoxia (P7). Table S3. Numbers of m A peaks located in different Publisher’sNote regions of mRNA transcripts in wild-type mouse cerebellum at P7, P14, Springer Nature remains neutral with regard to jurisdictional claims in P21, and P60. Table S9. List of antibodies and their applications used in published maps and institutional affiliations. this study. Table S10. List of primers for RT-qPCR used in this study. (PDF 14643 kb) Author details Additional file 2: Table S4. GO analysis of genes containing m AON Department of Pathology, Institute of Basic Medical Sciences Chinese and OFF switches during mouse cerebellar development. (XLSX 564 kb) Academy of Medical Science, School of Basic Medicine Peking Union Medical College; Neuroscience Center, Chinese Academy of Medical Sciences, Beijing Additional file 3: Table S5. GO analysis of genes encoded by the CMRs 100005, China. BIG Data Center, Beijing Institute of Genomics, Chinese at P7 and P60, and SMRs at the four developmental stages. (XLSX 128 kb) Academy of Sciences, Beijing 100101, China. University of Chinese Academy Additional file 4: Table S6. GO and KEGG pathway enrichment of Sciences, People’s Republic of China, Beijing 100049, China. State Key analysis of the genes encoded by SMRs and CMRs compared between Laboratory of Molecular Oncology, National Cancer Center/Cancer Institute the m A peaks at P7 and P60 that were identified by using MACS2 and Hospital, Chinese Academy of Medical Sciences, Peking Union Medical software. (XLSX 303 kb) College, Beijing 100021, China. Department of Microbiology, Oslo University Additional file 5: Table S7 List of the 839 RNAs with strong positive or Hospital, Rikshospitalet, Oslo NO-0027, Norway. Department of Molecular negative correlation between their methylation levels and expression Medicine, Institute of Basic Medical Sciences, University of Oslo, NO-0317 levels. (XLSX 82 kb) Oslo, Norway. Additional file 6: Table S8 GO analysis of the genes with altered m A Received: 3 December 2017 Accepted: 26 April 2018 levels between wild-type and Alkbh5-deficient mice exposed to hypobaric hypoxia. (XLSX 83 kb) References Abbreviations 1. Yao B, Christian KM, He C, Jin P, Ming GL, Song H. Epigenetic mechanisms BP: Biological processes; CC: Cellular components; CDS: Coding sequence; in neurogenesis. Nat Rev Neurosci. 2016;17:537–49. CMRs: Continuously methylated RNAs; DMRs: Differentially methylated RNAs; 2. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar ECL: Enhanced Chemiluminescence; EGL: External granule cell layer; L, Osenberg S, Cesarkas K, Jacob-Hirsch J, Amariglio N, Kupiec M, et al. FDR: False discovery rate; GO: Gene Ontology; HE: Hematoxylin-eosin Topology of the human and mouse m6A RNA methylomes revealed by staining; IGL: Inner granule cell layer; KEGG: Kyoto Encyclopedia of Genes and m6A-seq. Nature. 2012;485:201–6. Ma et al. Genome Biology (2018) 19:68 Page 17 of 18 3. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. 26. Li HB, Tong J, Zhu S, Batista PJ, Duffy EE, Zhao J, Bailis W, Cao G, Comprehensive analysis of mRNA methylation reveals enrichment in 3' Kroehling L, Chen Y, et al. m6A mRNA methylation controls T cell UTRs and near stop codons. Cell. 2012;149:1635–46. homeostasis by targeting the IL-7/STAT5/SOCS pathways. Nature. 2017; 4. Niu Y, Zhao X, Wu YS, Li MM, Wang XJ, Yang YG. N6-methyl-adenosine 548:338–42. (m6A) in RNA: an old modification with a novel epigenetic function. 27. Hess ME, Hess S, Meyer KD, Verhagen LA, Koch L, Bronneke HS, Dietrich MO, Genomics Proteomics Bioinformatics. 2013;11:8–17. Jordan SD, Saletore Y, Elemento O, et al. The fat mass and obesity associated gene (Fto) regulates activity of the dopaminergic midbrain 5. Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2017;18:31–42. circuitry. Nat Neurosci. 2013;16:1042–8. 6. Cao G, Li HB, Yin Z, Flavell RA. Recent advances in dynamic m6A RNA 28. Li L, Zang L, Zhang F, Chen J, Shen H, Shu L, Liang F, Feng C, Chen D, Tao modification. Open Biol. 2016;6:160003 H, et al. Fat mass and obesity-associated (FTO) protein regulates adult 7. Zhao X, Yang Y, Sun BF, Shi Y, Yang X, Xiao W, Hao YJ, Ping XL, Chen YS, neurogenesis. Hum Mol Genet. 2017;26:2398-2411. Wang WJ, et al. FTO-dependent demethylation of N6-methyladenosine 29. Walters BJ, Mercaldo V, Gillon CJ, Yip M, Neve RL, Boyce FM, Frankland PW, regulates mRNA splicing and is required for adipogenesis. Cell Res. 2014; Josselyn SA. The role of the RNA demethylase FTO (Fat mass and obesity- 24:1403–19. associated) and mRNA methylation in hippocampal memory formation. 8. Xiao W, Adhikari S, Dahal U, Chen YS, Hao YJ, Sun BF, Sun HY, Li A, Ping XL, Neuropsychopharmacology. 2017;42:1502–10. Lai WY, et al. Nuclear m(6)A Reader YTHDC1 Regulates mRNA Splicing. Mol 30. Widagdo J, Zhao QY, Kempen MJ, Tan MC, Ratnu VS, Wei W, Leighton L, Cell. 2016;61:507–19. Spadaro PA, Edson J, Anggono V, Bredy TW. Experience-dependent accumulation of N6-methyladenosine in the prefrontal cortex is associated 9. Zheng G, Dahl JA, Niu Y, Fedorcsak P, Huang CM, Li CJ, Vagbo CB, Shi Y, with memory processes in mice. J Neurosci. 2016;36:6771–7. Wang WL, Song SH, et al. ALKBH5 is a mammalian RNA demethylase that 31. Yoon K-J, Ringeling FR, Vissers C, Jacob F, Pokrass M, Jimenez-Cyrus D, Su Y, impacts RNA metabolism and mouse fertility. Mol Cell. 2013;49:18–29. Kim N-S, Zhu Y, Zheng L, et al. Temporal control of mammalian cortical 10. Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, Fu Y, Parisien M, Dai Q, Jia neurogenesis by m6A methylation. Cell. 2017;171:877-889.e17. G, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature. 2014;505:117–20. 32. Weng YL, Wang X, An R, Cassin J, Vissers C, Liu Y, Liu Y, Xu T, Wang X, 11. Ke S, Pandya-Jones A, Saito Y, Fak JJ, Vagbo CB, Geula S, Hanna JH, Black Wong SZH, et al. Epitranscriptomic m(6)A regulation of axon regeneration DL, Darnell JE Jr, Darnell RB. m6A mRNA modifications are deposited in in the adult mammalian nervous system. Neuron. 2018;97:313–25. e316 nascent pre-mRNA and are not required for splicing but do specify 33. Goldowitz D, Hamre K. The cells and molecules that make a cerebellum. cytoplasmic turnover. Genes Dev. 2017;31:990–1006. Trends Neurosci. 1998;21:375–82. 12. Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, Weng X, Chen K, Shi H, 34. Szulwach KE, Li X, Li Y, Song CX, Wu H, Dai Q, Irier H, Upadhyay AK, Gearing He C. N(6)-methyladenosine modulates messenger RNA translation M, Levey AI, et al. 5-hmC-mediated epigenetic dynamics during postnatal efficiency. Cell. 2015;161:1388–99. neurodevelopment and aging. Nat Neurosci. 2011;14:1607–16. 13. Zhou J, Wan J, Gao X, Zhang X, Jaffrey SR, Qian SB. Dynamic m(6)A mRNA 35. Frank CL, Liu F, Wijayatunge R, Song L, Biegler MT, Yang MG, Vockley CM, methylation directs translational control of heat shock response. Nature. Safi A, Gersbach CA, Crawford GE, West AE. Regulation of chromatin 2015;526:591–4. accessibility and Zic binding at enhancers in the developing cerebellum. 14. Coots RA, Liu X-M, Mao Y, Dong L, Zhou J, Wan J, Zhang X, Qian S-B. Nat Neurosci. 2015;18:647–56. m6A facilitates eIF4F-independent mRNA translation. Mol Cell. 2017;68: 36. Chang M, Lv H, Zhang W, Ma C, He X, Zhao S, Zhang Z-W, Zeng Y-X, Song 504–14. e507 S, Niu Y, Tong W-M. Region-specific RNA m A methylation represents a 15. Xiang Y, Laurent B, Hsu CH, Nachtergaele S, Lu Z, Sheng W, Xu C, Chen H, new layer of control in the gene regulatory network in the mouse brain. Ouyang J, Wang S, et al. RNA m6A methylation regulates the ultraviolet- Open Biol. 2017;7:170166 induced DNA damage response. Nature. 2017;543:573–6. 37. Barbaric I, Miller G, Dear TN. Appearances can be deceiving: phenotypes of 16. Batista PJ, Molinie B, Wang J, Qu K, Zhang J, Li L, Bouley DM, Lujan E, Haddad knockout mice. Brief Funct Genomic Proteomic. 2007;6:91–103. B, Daneshvar K, et al. m(6)A RNA modification controls cell fate transition in 38. Bunn HF, Poyton RO. Oxygen sensing and molecular adaptation to hypoxia. mammalian embryonic stem cells. Cell Stem Cell. 2014;15:707–19. Physiol Rev. 1996;76:839–85. 17. Aguilo F, Zhang F, Sancho A, Fidalgo M, Di Cecilia S, Vashisht A, Lee DF, 39. Zhang C, Samanta D, Lu H, Bullen JW, Zhang H, Chen I, He X, Semenza GL. Chen CH, Rengasamy M, Andino B, et al. Coordination of m(6)A mRNA Hypoxia induces the breast cancer stem cell phenotype by HIF-dependent methylation and gene transcription by ZFP217 regulates pluripotency and and ALKBH5-mediated m6A-demethylation of NANOG mRNA. Proc Natl reprogramming. Cell Stem Cell. 2015;17:689–704. Acad Sci U S A. 2016;113:E2047–56. 18. Chen T, Hao YJ, Zhang Y, Li MM, Wang M, Han W, Wu Y, Lv Y, Hao J, Wang 40. Jakobsen L, Vanselow K, Skogs M, Toyoda Y, Lundberg E, Poser I, Falkenby L, et al. m(6)A RNA methylation is regulated by microRNAs and promotes LG, Bennetzen M, Westendorf J, Nigg EA, et al. Novel asymmetrically reprogramming to pluripotency. Cell Stem Cell. 2015;16:289–301. localizing components of human centrosomes identified by complementary 19. Zhang C, Chen Y, Sun B, Wang L, Yang Y, Ma D, Lv J, Heng J, Ding Y, Xue Y, proteomics methods. EMBO J. 2011;30:1520–35. et al. m6A modulates haematopoietic stem and progenitor cell 41. Sarzi E, Angebault C, Seveno M, Gueguen N, Chaix B, Bielicki G, Boddaert N, specification. Nature. 2017;549:273–6. Mausset-Bonnefont AL, Cazevieille C, Rigau V, et al. The human OPA1delTTAG mutation induces premature age-related systemic 20. Zheng Q, Hou J, Zhou Y, Li Z, Cao X. The RNA helicase DDX46 inhibits neurodegeneration in mouse. Brain. 2012;135:3599–613. innate immunity by entrapping m6A-demethylated antiviral transcripts in the nucleus. Nat Immunol. 2017;18:1094–103. 42. Cui C, Chatterjee B, Lozito TP, Zhang Z, Francis RJ, Yagi H, Swanhart LM, 21. Hsu PJ, Zhu Y, Ma H, Guo Y, Shi X, Liu Y, Qi M, Lu Z, Shi H, Wang J, et al. Sanker S, Francis D, Yu Q, et al. Wdpcp, a PCP protein required for Ythdc2 is an N6-methyladenosine binding protein that regulates ciliogenesis, regulates directional cell migration and cell polarity by direct mammalian spermatogenesis. Cell Res. 2017;27:1115–27. modulation of the actin cytoskeleton. PLoS Biol. 2013;11:e1001720. 22. Lin Z, Hsu PJ, Xing X, Fang J, Lu Z, Zou Q, Zhang KJ, Zhang X, Zhou Y, 43. Jiang D, Zhao L, Clish CB, Clapham DE. Letm1, the mitochondrial Ca2+/H+ Zhang T, et al. Mettl3−/Mettl14-mediated mRNA N6-methyladenosine antiporter, is essential for normal glucose metabolism and alters brain modulates murine spermatogenesis. Cell Res. 2017;27:1216-1230. function in Wolf-Hirschhorn syndrome. Proc Natl Acad Sci U S A. 2013;110: 23. Xu K, Yang Y, Feng GH, Sun BF, Chen JQ, Li YF, Chen YS, Zhang XX, Wang E2249–54. CX, Jiang LY, et al. Mettl3-mediated m6A regulates spermatogonial 44. Sato A, Sekine Y, Saruta C, Nishibe H, Morita N, Sato Y, Sadakata T, differentiation and meiosis initiation. Cell Res. 2017;27:1100–14. Shinoda Y, Kojima T, Furuichi T. Cerebellar development transcriptome database (CDT-DB): profiling of spatio-temporal gene expression during 24. Zhao BS, Wang X, Beadell AV, Lu Z, Shi H, Kuuspalu A, Ho RK, He C. m6A- the postnatal development of mouse cerebellum. Neural Netw. 2008;21: dependent maternal mRNA clearance facilitates zebrafish maternal-to- 1056–69. zygotic transition. Nature. 2017;542:475–8. 25. Ivanova I, Much C, Di Giacomo M, Azzi C, Morgan M, Moreira PN, Monahan 45. Geula S, Moshitch-Moshkovitz S, Dominissini D, Mansour AA, Kol N, Salmon- J, Carrieri C, Enright AJ, O'Carroll D. The RNA m6A reader YTHDF2 is Divon M, Hershkovitz V, Peer E, Mor N, Manor YS, et al. Stem cells. m6A essential for the post-transcriptional regulation of the maternal mRNA methylation facilitates resolution of naive pluripotency toward transcriptome and oocyte competence. Mol Cell. 2017;67:1059–67. e1054 differentiation. Science. 2015;347:1002–6. Ma et al. Genome Biology (2018) 19:68 Page 18 of 18 46. He S, Wang H, Liu R, He M, Che T, Jin L, Deng L, Tian S, Li Y, Lu H, et al. 68. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, mRNA N6-methyladenosine methylation of postnatal liver development in Singh H, Glass CK. Simple combinations of lineage-determining pig. PLoS One. 2017;12:e0173421. transcription factors prime cis-regulatory elements required for macrophage 47. Zhang S, Zhao BS, Zhou A, Lin K, Zheng S, Lu Z, Chen Y, Sulman EP, Xie K, and B cell identities. Mol Cell. 2010;38:576–89. Bogler O, et al. m(6)A Demethylase ALKBH5 maintains tumorigenicity of 69. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing glioblastoma stem-like cells by sustaining FOXM1 expression and cell genomic features. Bioinformatics. 2010;26:841–2. proliferation program. Cancer Cell. 2017;31:591–606. e596 70. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, Liu XS. Model-based analysis of ChIP-Seq 48. Ke S, Alemu EA, Mertens C, Gantman EC, Fak JJ, Mele A, Haripal B, Zucker-Scharff (MACS). Genome Biol. 2008;9:R137. I, Moore MJ, Park CY, et al. A majority of m6A residues are in the last exons, 71. Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, allowing the potential for 3' UTR regulation. Genes Dev. 2015;29:2037–53. Mesirov JP. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6. 49. Luo GZ, MacQueen A, Zheng G, Duan H, Dore LC, Lu Z, Liu J, Chen K, Jia G, 72. Thorvaldsdottir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer Bergelson J, He C. Unique features of the m6A methylome in Arabidopsis (IGV): high-performance genomics data visualization and exploration. Brief thaliana. Nat Commun. 2014;5:5630. Bioinform. 2013;14:178–92. 50. Li Y, Wang X, Li C, Hu S, Yu J, Song S. Transcriptome-wide N(6)- 73. Kolde R: pheatmap. 2015. https://CRAN.R-project.org/package=pheatmap. methyladenosine profiling of rice callus and leaf reveals the presence of Published 11 Dec 2015. tissue-specific competitors involved in selective mRNA modification. RNA 74. Huang d W, Sherman BT, Lempicki RA. Systematic and integrative analysis of Biol. 2014;11:1180–8. large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57. 51. Tao X, Chen J, Jiang Y, Wei Y, Chen Y, Xu H, Zhu L, Tang G, Li M, Jiang A, 75. Wickham H: ggplot2: elegant graphics for data analysis. New York: Springer; et al. Transcriptome-wide N 6 -methyladenosine methylome profiling of porcine muscle and adipose tissues reveals a potential mechanism for 76. Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic transcriptional regulation and differential methylation pattern. BMC genome annotation and pathway reconstruction server. Nucleic Acids Res. Genomics. 2017;18:336. 2007;35:W182–5. 52. Linder B, Grozhik AV, Olarerin-George AO, Meydan C, Mason CE, Jaffrey SR. 77. Huang d W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: Single-nucleotide-resolution mapping of m6A and m6Am throughout the paths toward the comprehensive functional analysis of large gene lists. transcriptome. Nat Methods. 2015;12:767–72. Nucleic Acids Res. 2009;37:1–13. 53. LoVerso PR, Cui F. Cell type-specific transcriptome profiling in mammalian 78. Chunhui Ma, Mengqi Chang, Hongyi Lv, Yamei Niu, Shuhui Song, Wei-Min brains. Front Biosci (Landmark Ed). 2016;21:973–85. Tong, et al: RNA m6A methylation participates in regulation of postnatal 54. Dominissini D, Moshitch-Moshkovitz S, Salmon-Divon M, Amariglio N, development of the mouse cerebellum. SRA PRJNA400297. 2017. https:// Rechavi G. Transcriptome-wide mapping of N(6)-methyladenosine by www.ncbi.nlm.nih.gov/sra/?term=PRJNA400297. Accessed 25 Apr 2018. m(6)A-seq based on immunocapturing and massively parallel sequencing. 79. Wang Y, Song F, Zhu J, Zhang S, Yang Y, Chen T, Tang B, Dong L, Ding N, Nat Protoc. 2013;8:176–89. Zhang Q, et al. GSA: Genome Sequence Archive. Genomics Proteomics 55. Sanchez-Ortiz E, Cho W, Nazarenko I, Mo W, Chen J, Parada LF. NF1 Bioinformatics. 2017;15:14–8. regulation of RAS/ERK signaling is required for appropriate granule neuron 80. Members BIGDC. The BIG Data Center: from deposition to integration to progenitor expansion and migration in cerebellar development. Genes Dev. translation. Nucleic Acids Res. 2017;45:D18–24. 2014;28:2407–20. 81. Chunhui Ma, Mengqi Chang, Hongyi Lv, Yamei Niu, Shuhui Song, Wei-Min 56. Oomoto I, Suzuki-Hirano A, Umeshima H, Han YW, Yanagisawa H, Carlton P, Tong, et al: RNA m6A methylation participates in regulation of postnatal Harada Y, Kengaku M, Okamoto A, Shimogori T, Wang DO. ECHO-liveFISH: development of the mouse cerebellum. GSA CRA000472. 2017. http://bigd. in vivo RNA labeling reveals dynamic regulation of nuclear RNA foci in big.ac.cn/search?dbId=gsa&q=CRA000472. Accessed 25 Apr 2018. living tissues. Nucleic Acids Res. 2015;43:e126. 57. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R36. 58. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25:1105–11. 59. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7:562–78. 60. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5. 61. Saeed AI, Bhagabati NK, Braisted JC, Liang W, Sharov V, Howe EA, Li J, Thiagarajan M, White JA, Quackenbush J. TM4 microarray software suite. Methods Enzymol. 2006;411:134–93. 62. Dobin A, Gingeras TR. Optimizing RNA-Seq mapping with STAR. Methods Mol Biol. 2016;1415:245–62. 63. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21. 64. Anders S, Pyl PT, Huber W. HTSeq-a Python framework to work with high- throughput sequencing data. Bioinformatics. 2015;31:166–9. 65. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40. 66. Meng J, Cui X, Rao MK, Chen Y, Huang Y. Exome-based analysis for RNA epigenome sequencing data. Bioinformatics. 2013;29:1565–7. 67. Liu L, Zhang SW, Zhang YC, Liu H, Zhang L, Chen R, Huang Y, Meng J. Decomposition of RNA methylome reveals co-methylation patterns induced by latent enzymatic regulators of the epitranscriptome. Mol BioSyst. 2015;11: 262–74. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology Springer Journals

RNA m6A methylation participates in regulation of postnatal development of the mouse cerebellum

Free
18 pages

Loading next page...
 
/lp/springer_journal/rna-m6a-methylation-participates-in-regulation-of-postnatal-McbWCl1s49
Publisher
Springer Journals
Copyright
Copyright © 2018 by The Author(s).
Subject
Life Sciences; Animal Genetics and Genomics; Human Genetics; Plant Genetics and Genomics; Microbial Genetics and Genomics; Bioinformatics; Evolutionary Biology
eISSN
1474-760X
D.O.I.
10.1186/s13059-018-1435-z
Publisher site
See Article on Publisher Site

Abstract

6 6 Background: N -methyladenosine (m A) is an important epitranscriptomic mark with high abundance in the brain. Recently, it has been found to be involved in the regulation of memory formation and mammalian cortical neurogenesis. However, while it is now established that m A methylation occurs in a spatially restricted manner, its functions in specific brain regions still await elucidation. Results: We identify widespread and dynamic RNA m A methylation in the developing mouse cerebellum and further uncover distinct features of continuous and temporal-specific m A methylation across the four postnatal developmental processes. Temporal-specific m A peaks from P7 to P60 exhibit remarkable changes in their distribution patterns along the mRNA transcripts. We also show spatiotemporal-specific expression of m A writers METTL3, METTL14, and WTAP and erasers ALKBH5 and FTO in the mouse cerebellum. Ectopic expression of METTL3 mediated by lentivirus infection leads to disorganized structure of both Purkinje and glial cells. In addition, under hypobaric hypoxia exposure, Alkbh5-deletion causes abnormal cell proliferation and differentiation in the cerebellum through disturbing the balance of RNA m A methylation in different cell fate determination genes. Notably, nuclear export of the hypermethylated RNAs is enhanced in the cerebellum of Alkbh5-deficient mice exposed to hypobaric hypoxia. Conclusions: Together, our findings provide strong evidence that RNA m A methylation is controlled in a precise spatiotemporal manner and participates in the regulation of postnatal development of the mouse cerebellum. Keywords: N -methyladenosine, RNA methylation, ALKBH5, METTL3, Cerebellar development Background epitranscriptomic mark found in mRNAs [2, 3]. Similar to Epigenetic regulation, including histone modifications, DNA and protein modifications, m Ain mRNA is DNA modifications, and non-coding RNAs, plays crucial reversible and thus open to dynamic regulation [4, 5]. roles in both embryonic and adult neurogenesis [1]. The Levels of RNA methylation are finely balanced through an 6 6 most recently discovered regulatory modification, N - interplay among m A methyltransferases (writers), methyladenosine (m A), is a highly abundant, so-called demethylases (erasers), and binding proteins (readers) [6]. m A has been shown to regulate RNA processing, including RNA splicing [7, 8], nuclear export [9], RNA * Correspondence: niuym@ibms.pumc.edu.cn; songshh@big.ac.cn; wmtong@ibms.pumc.edu.cn degradation [10, 11], and translation [12–14]. In Chunhui Ma, Mengqi Chang and Hongyi Lv contributed equally to this agreement with its functions in RNA metabolism, m Ais work. 1 an important regulatory factor in different cellular Department of Pathology, Institute of Basic Medical Sciences Chinese Academy of Medical Science, School of Basic Medicine Peking Union Medical processes, including heat shock response [13], DNA College; Neuroscience Center, Chinese Academy of Medical Sciences, Beijing damage response [15], cell fate determination [16–19], 100005, China 2 and innate immunity [20]. In addition to its functions in BIG Data Center, Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing 100101, China these cellular processes, in vivo studies in model Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Ma et al. Genome Biology (2018) 19:68 Page 2 of 18 6 6 organisms have revealed crucial roles of m Ain m A peaks (Additional file 1: Figure S1d). Distribution spermatogenesis [9, 21–23], embryogenesis [24, 25], and analysis of all the m A peaks in mRNAs revealed that T-cell homeostasis [26]. their methylation patterns were similar from P7 to P60 Although m A is abundant in the brain, its dynamic (Additional file 1: Figure S1e). Intriguingly, we observed regulation and biological significance for brain development significant enrichment of m A peaks in start codon remain largely unknown [3]. Among the m A-related genes regions in addition to the coding sequence (CDS) and identified, the RNA demethylase FTO is proven to be stop codon regions. A statistical analysis of the m A peaks crucial for various brain functions, including D2R-D3R- located in different transcript segments showed that the GIRK-mediated signaling [27], adult neurogenesis [28], and proportion of m A peaks in start codon regions increased memory formation [29, 30]. A recent report confirms the gradually from P7 to P60 (Additional file 1: Figure S1f and role of m A methylation in temporal progression of Additional file 1: Table S3). mammalian cortical neurogenesis and in regulating axon To explore how m A changes in postnatal mouse regeneration in the adult mammalian nervous system [31, cerebellum, m A peaks were subjected to pairwise 32]. To better understand the underlying mechanisms comparison between each two adjacent stages. In total, responsible for the temporal and spatial complexity of the we identified 12,452 “ON” switches (emerging m A brain, it is important to characterize the physiological peaks in the later stage) and 11,192 “OFF” switches 6 6 functions of m A in a context-dependent manner. However, (resolving m A peaks from the former to the later stage) the biological functions of other m A-related genes in (P < 0.05, enrichment score > 1.5) during the postnatal different brain regions remain unclear. development (Fig. 1a, b and Additional file 1: Figure The cerebellum consists of a multi-layered structure, S1g). The m A ON and OFF switches also exhibited rendering it an important model for the study of neur- apparent differences in their distribution patterns. The onal cell proliferation, migration, and differentiation m A ON switches were more likely to be enriched in [33]. Considerable progress has been made in identifying start codon regions rather than stop codon regions. In the epigenetic mechanisms responsible for cerebellar de- contrast, the numbers of m A OFF switches velopment, such as DNA 5-hydroxymethylation [34] and surrounding stop codon regions were higher than that of chromatin accessibility [35]. Interestingly, the cerebel- ON switches (Fig. 1c and Additional file 1: Figure S1h). lum of Fto-deficient mouse was found to be smaller than Together, these data imply that m A marks occurring at that of wild-type mouse [28]; however, mechanisms re- earlier or later stages of cerebellar development have sponsible for this growth defect remain unclear. In different preferences with regard to their distribution agreement with the spatial and temporal regulation in along with the whole mRNA transcripts. the brain, cerebellum-specific methylation has been To evaluate the biological significance of genes with dy- identified [36]. Despite these findings, the precise role of namic RNA m A modification, we next performed Gene RNA methylation in the cerebellum and its temporal Ontology (GO) analysis for those genes with ON/OFF control during postnatal mouse cerebellar development switches (Additional file 2: Table S4). The genes with still await elucidation. different types of switches appeared to possess different Here we investigated the dynamic regulation of m Ain functions. For example, the genes containing P7–P14 OFF the developing mouse cerebellum and determined its switches were mostly annotated to biological processes biological significance by disturbing the expression of the such as cell cycle, cell division, and DNA repair (Fig. 1d). 6 6 m AwriterMETTL3orthe m AeraserALKBH5in mice In contrast, the newly occurring methylation at P14, P21, exposed to hypobaric hypoxia. This study uncovered that the and P60 represented by the ON switches was detected in temporal regulation of RNA m A methylation is essential for many genes involved in signal transduction, cell adhesion, precise control of postnatal cerebellar development. learning, and synaptic plasticity (Fig. 1e). Together, these data confirm that extensive RNA methylation and Results demethylation occur over the course of neuronal Temporal regulation of m A methylation in postnatal differentiation in vivo, in both proliferating and fully mouse cerebellum differentiated neural cells. Moreover, the distinct functions 6 6 To understand the biological impact of m A in mouse of genes containing m A OFF or ON switches suggest 6 6 cerebellum, we first performed m Aanalysisusing that m A is a prerequisite for those genes to exert their postnatal cerebella at day 7 (P7), P14, P21, and P60 functions at each developmental stage. (Additional file 1: Table S1). We identified 10,449, 10,266, 10,419, and 10,783 methylated poly(A) RNAs at the four Temporal-specific m A methylation acts in concert with stages, respectively (Additional file 1: Table S2, Additional cerebellar developmental control file 1: Figure S1a–c). Data analysis showed that UGGACU As m A is developmentally regulated in mouse cerebella, was the most conserved consensus sequence among all we next analyzed the methylation profiles over the four Ma et al. Genome Biology (2018) 19:68 Page 3 of 18 Fig. 1 Developmentally regulated RNA methylation during mouse cerebellar development. a Numbers of ON or OFF m A switches between P7, P14, P21, and P60 cerebella pairwise comparisons. The numbers of the poly(A) RNAs containing the m A switches are included in parentheses. b IGV plots 6 6 showing examples of ON and OFF m A switches at each transition period. Grey reads originate from input libraries and red reads originate from m A- IP libraries. Y-axis represents normalized numbers of reads count. Arrows indicate the direction of gene transcription. Dashed boxes indicate the position of ON/OFF switches. c Distribution of each type of ON and OFF m A switch along the whole mRNA transcripts. d, e Most impacted GO biological process terms of the methylated RNAs containing OFF (d)or ON(e)m A switches across the four developmental stages postnatal stages. We found 8367 continuously methylated direction (Fig. 2c and Additional file 1: Figure S2b–c). RNAs (CMRs) throughout cerebellar development, Given the different distribution between ON and OFF together with 634, 260, 315, and 512 specifically switches (Fig. 1c and Additional file 1: Figure S1h), we methylated RNAs (SMRs) at P7, P14, P21, and P60, further analyzed the distribution of m A peaks in respectively (Fig. 2a, b and Additional file 1: Figure S2a). continuously and temporal-specifically methylated We then investigated whether the SMRs and CMRs mRNAs. In the mouse cerebellum of P7, the temporal- possessed different features in their methylation during specific m A peaks were significantly enriched at stop postnatal cerebellar development. Compared with the codon regions, while the P60-specific m A peaks were SMRs, the CMRs exhibited higher levels of both mostly enriched in start codon regions. The P14- or P21- methylation and expression throughout the specific m A peaks displayed an apparent transition from developmental process. Moreover, the methylation levels stop codon to start codon during progression of cerebellar of SMRs displayed a gradual reduction from P7 to P60, development (Fig. 2d). Furthermore, we quantified the while their expression levels changed in the opposite numbers of temporal-specific m A peaks located in the Ma et al. Genome Biology (2018) 19:68 Page 4 of 18 Fig. 2 Distinct features of temporal-specific m A methylation during postnatal development of the mouse cerebellum. a Venn diagram showing the numbers and relationship of methylated poly(A) RNAs in the mouse cerebellum from P7 to P60. b IGV plots showing examples of continuous and temporally specific methylation. Grey reads originate from input libraries and red reads originate from m A-IP libraries. Y-axis represents normalized numbers of read counts. Arrows indicate the direction of gene transcription. Dashed boxes indicate the position of temporal-specific 6 6 m A peaks. c Box plots showing the relative methylation levels of SMRs and CMRs as evaluated by the enrichment scores of m A peaks. SMR specifically methylated RNA, CMR continuously methylated RNA. ***P < 0.001 by Wilcoxon test. d Normalized distribution of temporal-specific m A peaks in each developmental stage along the whole mRNA transcripts. e, f Most impacted GO biological process terms (BP)(e) and cellular component terms (CC)(f) of SMRs in the four developmental stages five regions of mRNA transcripts. From P7 to P60, we the cells in the cerebellum at P60 are fully differentiated. observed a sharp increase in the numbers of temporal- To confirm that such a dynamic distribution pattern was specific m A peaks in start codon regions, as well as a related to the cell states, we compared the RNA methyla- decrease in the CDS and stop codon regions (Additional tion profiles in mouse cerebellum at P7 and P60 directly file 1: Figure S2d). To ensure that these changes were real using the peaks identified by exomePeak. Although the rather than artifacts, we randomly extracted the same proportion of m A peaks located in CDS were similar numbers of CMR peaks as those of temporal-specific between the two stages, we observed a higher percentage peaks and repeated the distribution analysis. Unlike the of P7-specific m A peaks near stop codon regions, but a temporal-specific peaks shown in Fig. 2d, these CMR higher percentage of P60-specific m A peaks in start peaks displayed similar distribution patterns from P7 to codon regions (Additional file 1: Figure S2f). These data P60 (Additional file 1: Figure S2e). The majority of cells in show that for those temporal-specifically methylated the cerebellum of P7 are proliferating and immature, while RNAs, the patterns of m A deposition onto RNA exhibit Ma et al. Genome Biology (2018) 19:68 Page 5 of 18 significant changes along with the progression of (Additional file 1: Figure S4c), among which 839 RNAs cerebellar development. showed either positive (674) or negative (165) correlation To explore the functional significance of temporal- between methylation and expression levels across the four specific methylation, we next performed GO analysis for stages (Additional file 5: Table S7). Results of the cluster those genes with continuous or temporal-specific methy- and subsequent GO analyses suggest that for these RNAs, lation. In line with the progression of cerebellar develop- m A participates in the developmental control of the ment, the four groups of genes containing temporal- mouse cerebellum through modulating their gene specific m A marks were annotated to different biological expression (Additional file 1: Figure S4d, e). processes (Fig. 2e and Additional file 3: Table S5). Genes with P7-specific methylation (such as Aspm in Fig. 2b) Altered expression of either METTL3 or ALKBH5 causes were enriched in processes such as cell cycle, cell division, defective mouse cerebellar development and DNA damage response, which are prerequisites in To gain insights into how the dynamic m A methyl marks proliferating neural cells. In contrast, many genes display- are regulated during cerebellar development, we analyzed ing P60-specific methylation (such as Glra1 in Fig. 2b) the protein expression profiles of three m A writers were annotated to biological processes including trans- (METTL3, METTL14, and WTAP) and two m A erasers port, oxidation-reduction process, and metabolic pro- (ALKBH5 and FTO). All five proteins were highly cesses, which are required for mature neuronal activities. expressed at the early stage of cerebellar development (P7) In agreement with their diverging biological functions, the and showed a gradual reduction towards the maturation genes with temporal-specific methylation were also char- of cerebellar neurons (P60) (Fig. 3a). Given the acteristic of their distinct cellular localization (Fig. 2f). In spatiotemporal specificity of gene expression in the brain, contrast, genes encoded by the P7 and P60 CMRs exhib- we performed immunohistochemical staining to detect ited high similarity in their annotation of enriched bio- their expression in situ. In the cerebellum at P7, all five logical processes and cellular components (Additional file proteins were detected in the external granule cell layer 1: Figure S2g and Additional file 3: Table S5). Together, (EGL), Purkinje cell layer (PCL), and inner granule cell these results suggest that the majority of genes with con- layer (IGL); however, their expression levels varied among tinuous methylation play a fundamental role throughout different types of cells (Fig. 3b). Furthermore, along with cerebellar development, while temporal-specific methyla- the progression of cerebellar development, we observed a tion only occurs at specified times for those genes, thus reduction of their expression levels in internal granular ensuring the proper progression of cerebellar development layers but an elevation in Purkinje cells (Fig. 3c). in a temporal-dependent manner. In agreement with the dynamic expression of the 6 6 6 In parallel to the methylation analysis using the m A m A-related proteins, the global m A levels of poly(A) peaks identified using exomePeak, we re-analyzed the RNA also decreased from P7 to P60 (Fig. 3d). The 6 6 6 m A-seq data at P7 and P60 using another canonical highest m A level at P7 suggests that m A methylation peak-calling algorithm, MACS2. The P7 SMRs, P60 plays an important role at the earlier developmental SMRs, and CMRs exhibited distinct features in their rela- stage. To confirm this, we reduced local m A levels in tive methylation levels, distribution patterns of the m A the developing cerebellum (P7) by knocking down peaks, and their functional annotations (Additional file 1: Mettl3 using lentivirus system and examined the Figure S3 and Additional file 4: Table S6). The very high consequent morphological changes (Fig. 4a–d). Purkinje similar results from the two sets of independent analysis cells in the control group displayed an orderly alignment further strengthened our awareness of the importance of along the outer surface of IGLs with dendritic temporal specific m A methylation during mouse arborization. However, knockdown of Mettl3 resulted in cerebellar development. We then evaluated how the RNA a severe alteration in Purkinje cell numbers and laminal methylation impacted on gene expression. When each structure and and stunted dendrites (Fig. 4e). Moreover, sample was analyzed individually, RNA methylation and GFAP immunostaining revealed that glial cell fibers expression levels at the global scale were negatively were severely disorganized in Mettl3-knocked down correlated with each other (Additional file 1: Figure S4a). regions (Fig. 4f). We also examined the effect of Mettl3 Next, pairwise comparison of RNA expression between overexpression on cerebellar development (Additional each two adjacent stages identified a total of 2451 (P value file 1: Figure S5a, b). Although the global m A levels < 0.05) differentially expressed RNAs throughout increased modestly (Additional file 1: Figure S5c), development (Additional file 1: Figure S4b), among which overexpression of Mettl3 led to apparent morphology approximately 90% were validated by another parallel changes, as was observed in the Mettl3-knocked down differential expression analysis using a STAR/edgeR cerebellum (Additional file 1: Figure S5d-S5e), suggest- package. Cluster analysis of all expressed RNAs from P7 ing that appropriate m A levels are crucial for cerebellar to P60 produced six clusters of gene expression patterns development. Ma et al. Genome Biology (2018) 19:68 Page 6 of 18 Fig. 3 Dynamic expression of m A-related genes and RNA methylation levels during mouse cerebellar development. a Western blot analysis showing the expression levels of m A methyltransferases and demethylases in the mouse cerebellum from P7 to P60. n = 12 biological replicates. GAPDH is used as an internal control. b, c Immunostaining of mouse cerebellum at P7 and P60 using antibodies against the five m A-related proteins. n = 4 biological replicates. Enlarged images in the dashed boxes are shown in c. Scale bar in b represents 50 μm, and scale bar in c represents 10 μm. EGL external granule cell layer, IGL internal granule cell layer, PCL Purkinje cell layer, ML molecular layer. The dashed circles indicate the representative Purkinje cells. d Dynamic poly(A) RNA m A methylation levels during mouse cerebellar development. Mouse cerebella at P7, P14, P21, and P60 were included in this assay. Experiments were repeated three times using 15 mice in total. Representative data are shown here. ***P < 0.001 Given the spatiotemporal expression of ALKBH5 in the speculated that this phenomenon was likely due to genetic mouse cerebellum, we next tested its role in cerebellar de- redundancy and a compensatory response [37]. To over- velopment. Compared to that of wild-type (WT) mice, the come these compensatory mechanisms, we made use of the cerebellum of Alkbh5-knockout (KO) mice lacked any de- fact that the mammalian brain is sensitive to neuronal dam- tectable changes in weight and morphology (Additional file age by hypoxia [38], while ALKBH5 is unique among the 1: Figure S6). Considering the potential functions of FTO ALKB families in its response to hypoxia stimulation [39]. and other undefined demethylases in the cerebellum, we Therefore, we challenged the mice to extreme physiological Ma et al. Genome Biology (2018) 19:68 Page 7 of 18 Fig. 4 Altered morphology in the mouse cerebellum infected with lentivirus for Mettl3 knockdown. a, b GFP fluorescence image (a) and hematoxylin-eosin (HE) staining (b) of the mouse cerebellum dissected 10-days post-lentivirus infection. Arrows indicate the position of injection. Scale bar in a represents 1 mm. Scale bar in b represents 200 μm. c Western blot analysis to confirm the knockdown efficiency of Mettl3 in the neuro2a cell line. β-ACTIN was used as an internal control. sh-1 and sh-2 represent two kinds of Mettl3 knockdown lentiviral vectors. sh-1 was used for subsequent in vivo analysis. d Decreased m A levels resulting from Mettl3 knocking down as analyzed using UHPLC-MS/MS. n =4. *P < 0.05. e, f Immunofluorescence staining with antibodies against Calbindin-D28K (e) and GFAP (f) was performed to detect Purkinje cells (e) and astrocytes (f) upon lentivirus infection (n = 3 biological replicates). GFP fluorescence indicates the area with Mettl3 knockdown. Sections were counterstained with DAPI to visualize the nuclei. Scale bar represents 100 μm conditions in order to exhaust the potential compensatory Alkbh5-deficient cerebellum was significantly reduced as mechanisms. We propose that ALKBH5 might be involved reflected by NeuN immunostaining (Fig. 5f). These data in the protective mechanism of the brain against damage suggest that Alkbh5 deficiency disrupted the normal caused by hypoxia. To confirm this, WT and KO mice were progression of granular neurons from proliferation to exposed to hypobaric hypoxia for 48 h. We found that in differentiation. In addition, we found that Alkbh5 deficiency KO mice the sizes of the whole brain and the cerebellum considerably reduced dendritic arborization of Purkinje were significantly reduced compared to their littermate cells (Fig. 5g), concomitant with an increase in controls (Fig. 5a, b). Immunostaining analysis revealed a disorganization of the radial fibers in glial cells (Fig. 5h). significant increase in the numbers of Ki67 proliferating Following hypoxic treatment and return to normoxic cells and phospho-histone 3 (PH3 ) mitotic cells in the conditions, we examined the cerebellar sections again at EGL compared to the WT counterparts (Fig. 5c, d). another stage (P14). The internal granule neurons, Purkinje Consistently, we also detected an increase in the number of cells, and glial cells exhibited no visible difference between cells in the S phase of the cell cycle (positive BrdU WT and KO mice (Additional file 1:FigureS7a–c). immunoactivity) in the cerebellum of KO mice (Fig. 5e). However, we could still observe a slight increase in the Furthermore, the number of mature neurons in the IGL of numbers of proliferating granule cells as indicated by Ki67 Ma et al. Genome Biology (2018) 19:68 Page 8 of 18 Fig. 5 Defective cerebellar development in Alkbh5 knockout mice after being exposed to hypobaric hypoxia for 48 h. a Brain weight of individual wild-type (WT) and Alkbh5 knockout (KO) mice. Numbers of biological replicates are included in parentheses. Black lines in the graph indicate the mean. ***P < 0.001. b Representative HE staining images of wild-type (WT) and Alkbh5 knockout (KO) mouse cerebellum (P7). n = 7 for WT and KO, respectively. Scale bar, 100 μm. c–f Representative images from immunostaining analysis with antibodies against Ki67 (c), phospho-H3 (PH3) (d), BrdU (e), and NeuN (f) in the cerebellum of WT and KO mice. Quantification of immuno-reactive cells is shown in the right panels. Numbers of biological replicates are included in parentheses. Black lines in the graph indicate the mean. *P < 0.05, **P < 0.01, ns not significant. Sections in c and e were stained with DAPI and sections in d and f were counterstained with eosin to visualize nuclei. Scale bar, 100 μm. g, h Representative images from immunohistochemical analysis with antibodies against Calbindin-D28K (g) and GFAP (h) to detect Purkinje cells (g) and astrocytes (h). n = 7 for WT and KO, respectively. Scale bar, 100 μm. *P < 0.05, **P < 0.01, ***P < 0.001, ns not significant and PH3 immunoreactivity (Additional file 1:FigureS7d, deficiency demonstrates that proper orchestration of the e). In addition, the cerebellum in Alkbh5 KO mice spatiotemporal expression of m A writers and erasers is remained smaller than that of the littermate controls necessary for mouse cerebellar development. (Additional file 1: Figure S7f). These data indicate that in spite of the undefined compensatory response to partially Imbalance of RNA m A methylation in the Alkbh5-deficient recover the morphology changes, Alkbh5 deficiency had a mouse cerebellum exposed to hypobaric hypoxia profound and deleterious effect on cerebellar development To examine how Alkbh5 deficiency affects m A under hypoxic conditions. methylation of the cerebellar poly(A) RNA, we first Together, the abnormal cerebellar development resulting performed UHPLC-MS/MS analysis on the poly(A) from either ectopic expression of Mettl3 or Alkbh5 RNA from cerebella at P7 and found a slight increase in Ma et al. Genome Biology (2018) 19:68 Page 9 of 18 6 6 the global m A levels in the KO mice exposed to increase of m A at the global level, Alkbh5 deficiency hypobaric hypoxia (Additional file 1: Figure S8a). To led to disordered m A levels of thousands of RNAs. gain more detailed insight into the epitranscriptome- We next performed GO analysis to evaluate whether wide changes of all methylated RNAs, we performed the defective cerebellar development in KO mice re- 6 6 m A-seq using the poly(A) RNAs from both WT and sulted from imbalanced m A methylation (Fig. 6c, d and KO mice exposed to hypobaric hypoxia (Additional file Additional file 1: Figure S8e). Partial or complete loss of 1: Table S1). A slight increase in the number of m A methylation mainly occurred to the genes participating peaks and methylated RNAs was observed in the in the control of cell division (such as Cenpe in cerebella of KO mice (Additional file 1: Table S2), as Additional file 1: Figure S8c), cell cycle (such as Cdca2 well as higher enrichment scores of m A peaks (Fig. 6a). in Additional file 1: Figure S8d), and cell projection In addition, we identified 1348 poly(A) RNAs with gain organization (such as Erbb4 in Additional file 1: Figure of methylation and 711 poly(A) RNAs with loss of S8d). In contrast, the hypermethylated RNAs were methylation in KO mice cerebellum (P7) (Fig. 6b; found to be related to metabolic processes (such as Ccr5 Additional file 1: Figure S8b, c). By comparing those in Additional file 1: Figure S8c), ion transport (such as abnormally methylated RNAs to the SMRs during Camk2g in Additional file 1: Figure S8d), and axon cerebellar development of the WT mice, we found that guidance (such as Rora in Additional file 1: Figure S8d). upon Alkbh5 deficiency, 53 original P7 SMRs lost their Previous studies showed that ALKBH5 regulates RNA nu- methylation, while 42 original P14 SMRs, 53 P21 SMRs, clear export via an m A-dependent manner. We therefore and 104 P60 SMRs became methylated at P7. In asked whether the changes in m A methylation resulting addition to the gain- and loss-of-methylation RNAs, we from Alkbh5 deficiency also affect the RNA nuclear export also identified a considerable number of RNAs exhibit- in vivo. Based on functional enrichment analysis, we chose ing a significant increase (514) or decrease (81) in their several genes in different functional pathways, such as cell methylation levels (Fig. 6b). Therefore, in spite of a mild division (Mphosph9)[40], membrane potential (Opa1)[41], Fig. 6 Imbalanced m A RNA methylation in the cerebellum of Alkbh5-deficient mouse exposed to hypobaric hypoxia. a Boxplots showing the relative methylation levels of all m A peaks between wild-type (WT)and Alkbh5-deficient (KO) mouse cerebellum. ***P < 0.001 by Wilcoxon test. b Venn dia- gram showing the relationship among the RNAs containing gain-of-methylation (Gain), lost-of-methylation (Loss), and differentially methylated RNAs (DMRs). c, d GO analysis of genes encoded by RNAs containing loss- (c)or gain-of-methylation (d) in KO mouse cerebellum. BP biological process, CC cellular component. e Gene-specific m A-IP-qPCR results showing the relative methylation levels of four RNAs in the cerebellum of KO mice compared to those in WT mice. n =7. f Subcellular localization of the four RNAs as shown by the ratio of RNA abundance between the cytoplasmic and nuclear fractions isolated from the cerebellum of WT and KO mice. n = 10 for WT and 8 or 10 for KO, respectively. *P < 0.05, **P <0.01, ns not significant Ma et al. Genome Biology (2018) 19:68 Page 10 of 18 microtubule cytoskeleton organization (Wdpcp)[42], and improved our understanding of the mechanism related to ion transport (Letm1)[43], for further investigation. We first mouse cerebellar development [34, 35, 44]. Given the role 6 6 conducted gene-specific m A qPCR assays to validate the of m A in fate determination of embryonic stem cells and increased m A levels of these four mRNAs in the mammalian cortical neurogenesis [16, 17, 31, 45], we cerebellum of KO mice (Fig. 6e). Subsequently, subcellular hypothesized that m A also participates in neuronal cell fractionation analysis of cerebellar cells was performed proliferation and differentiation in the mouse cerebellum. (Additional file 1: Figure S8f, g). Compared to WT mice In this study, we show that a considerable number of cerebellum, all the four RNAs invariably exhibited increased RNAs in the mouse cerebellum exhibit developmentally abundance in the cytoplasm in the Alkbh5 deficient regulated methylation. GO analysis implied that cerebellum (Fig. 6f). Together, these data indicate that numerous genes regulating cell proliferation, Alkbh5 deficiency led to dysregulated RNA nuclear export differentiation, or metabolism acquire m A marks in a in the cerebellum. temporal-specific manner, among which the expression of Given the role of m A in regulating RNA stability, we 839 RNAs was under tight regulation by their methylation also analyzed the gene expression changes resulting from status. Notably, we found that the number of differentially Alkbh5 deficiency. Among the 497 differentially expressed methylated RNAs is much higher than that of differen- RNAs, 81 exhibited changed methylation and were tially expressed RNAs as observed during cerebellar devel- enriched in functional pathways, including transport, cell opment and in Alkbh5 KO mice. We speculate that this is division, and cell cycle control (Additional file 6:Table probably due to the additional functions of RNA methyla- S8). For instance, the genes involved in cell cycle control, tion in RNA splicing, nuclear export, and protein transla- such as Ddx11, Ccnb1,and Cbx1, displayed a reduction in tion [5]. For example, although the hypermethylated their RNA expression together with their reduced RNAs presented in Fig. 6e did not change in their expres- methylation. For genes associated with transport or sion levels, their nuclear export was enhanced by Alkbh5 metabolism, Alkbh5 deficiency-induced RNA hypermethy- deficiency. Since it is the transcripts in the cytoplasmic lation led to either a decrease (such as Cacna2d3, Notch3, pool that are utilized for RNA translation to protein and and Jam3) or increase (such as Slc6a20a, Dlk1, Slc2a4, RNA degradation, appropriate m A levels are necessary and Ascl2) in their RNA expression in the cerebellum for maintaining RNA nuclear export in a balanced state. (Additional file 1: Figure S8h). These results imply that the selective m A marks, through Taken together, our data confirm that, by altering the their various modulatory effects on RNA processing, original m A level, Alkbh5 deficiency disturbed the RNA function as epitranscriptomic switches that either activate metabolism of a subset of cell fate determination genes or suppress a series of physiological events during mouse and thus led to defective cerebellar development in the cerebellar development. In line with our findings in the mice exposed to hypobaric hypoxia. mouse cerebellum, dynamic changes in m ARNA methylation were also observed during postnatal liver Discussion development in pig [46] and mouse spermatogenesis [22]. Brain development is known to be regulated in a precise The dynamic RNA methylation levels in the developing spatiotemporal manner. We have recently identified cerebellum reveal how important it is to maintain the bal- region-specific m A with distinct methylation patterns ance of RNA methylation in the developing mouse cerebel- and biological functions for the mouse cerebellum and lum, which is accomplished by both the m A cerebral cortex [36]. In the present study, we further methyltransferases and demethylases. The morphological demonstrate that developmentally regulated m A changes observed upon transient knockdown or methylation is one of the important regulatory factors overexpression confirm that METTL3 is involved in for postnatal mouse cerebellar development. Temporal- regulating cerebellar development. In addition to specific methylation differs from continuous methylation methyltransferase, RNA m A demethylase is another in its methylation level, biological functions, and important factor to maintain RNA methylation at dynamic distribution of m A marks. Furthermore, we appropriate levels. The important role of FTO in the uncovered that abnormal expression of m A mouse cerebellum was shown in a previous report using a methyltransferases (METTL3) and demethylases Fto-knockout mouse model, which is characterized by a (ALKBH5) caused imbalanced RNA methylation, which cerebellum that is significantly smaller than that in WT further led to defective cerebellar development. We also litter mates [28]. ALKBH5, another demethylase identified discovered that Alkbh5 deficiency led to facilitated earlier, regulates spermatogenesis [20] and maintains nuclear export of those hypermethylated RNAs. tumorigenicity of glioblastoma stem-like cells [47]. How- Previous identified changes in genome-wide gene ex- ever, its function in the brain remains unclear. Here we pression and epigenetic modifications throughout the only observed defects in cerebellar development in Alkbh5- whole cerebellar developmental processes have greatly deficienct mice upon exposure to hypobaric hypoxia. It is Ma et al. Genome Biology (2018) 19:68 Page 11 of 18 likely that ALKBH5, FTO, and probably other undefined developmental process. Our results show that the positive demethylases overlap in their demethylation activity and or negative regulation of m A marks on these cell can partially substitute for each other’s functions under fate-determining genes is important to drive the proper normal conditions [37]. On the other side, our results also progression of postnatal cerebellar development. reveal the distinct functions between ALKBH5 and FTO. Different from previous in vitro analysis showing that In line with a previous report [9], we provide in vivo evi- m A marks are predominantly located near stop codons dence that ALKBH5 regulates the nuclear export of the and 3′ UTRs, previous in vivo methylation analysis of hypermethylated RNAs, while FTO regulates RNA splicing mouse cerebellum and cerebral cortex identified the in a methylation-dependent manner [7]. In this study, we m A marks in start codon regions as well as stop codons observed increased cell proliferation in the Alkbh5-defi- and 3′ UTRs [36]. Here, we identified similar patterns of cienct mouse cerebellum exposed to hypobaric hypoxia. In m A distribution throughout the process of postnatal contrast, Fto deficiency in adult mouse hippocampus led to cerebellar development. Importantly, we observed a reduced cell proliferation in adult neural stem cells [28]. significant change in the distribution of temporal- The discrepancy in these results confirmed that RNA specific m A marks from P7 to P60. We deduce that the methylation constitutes a complex process and suggests distribution patterns of temporal-specific m A peaks that the regulation and function of m A methylation in the possibly reflect the difference in methylation sites mouse brain should be investigated in a more precise between proliferating, immature and fully differentiated, spatiotemporal-dependent manner. mature neural cells. In the mouse cerebellum of P7, In the WT mouse cerebellum at P7, the m A during which the majority of neural cells are actively modification selectively marked a group of RNAs proliferating, the temporal-specific m A peaks were responsible for cell cycle control and cell division. In enriched in CDS, stop codons, and 3′ UTRs, which was Alkbh5 KO mice exposed to hypobaric hypoxia, we similar to the results obtained from the m A-seq detected many of these types of RNAs with aberrant analysis performed using cell lines [2, 3, 48]. In contrast, methylation, which may further disturb RNA metabolism, distribution of temporal-specific m A peaks in start such as RNA splicing, nuclear export, RNA decay, or codon regions was most significant in mouse cerebellum translation. In agreement with these findings, we observed at P60, which mainly consists of well-differentiated and an increase in the numbers of mitotic cells in the mature neural cells. In support of our findings, m A cerebellum of Alkbh5 KO mice upon exposure to marks in start codon regions were also observed in the hypobaric hypoxia. Mechanistically, we showed the fully differentiated tissues from mouse cerebral cortex hypermethylation and enhanced nuclear export of [36] (unpublished data), A. thaliana [49], rice [50], liver Mphosph9, which encodes a centrosomal protein related to [46], and muscles [51]. Collectively, these data imply cell division and cell cycle control [40]. Therefore, that the methylation in start codon regions occurs not abnormal RNA methylation is likely one of the reasons only in mature neural cells of the mouse cerebellum, but resulting in dysregulated cell proliferation and cell cycle also in many other types of fully differentiated cells in progression observed in the Alkbh5-deficient cerebellum. diverse organisms. We thus propose that the position of 6 6 Consistent with our findings, deficiency in the m A reader m A marks included in the transcript may be an Ythdf2 delayed DNA synthesis and mitotic processes in a important factor that determines the functions mediated zebrafish model [24]. In another study, deletion of Mettl14 by m A. Further investigation will be necessary to clarify in embryonic mouse brain also led to prolonged cell cycle the precise roles of m A in different physiological of cortical neural progenitors [31]. Together, these processes. However, the m A-IP-seq used in this study 6 6 phenotypes provide strong evidence that balanced m A was unable to identify the exact position of the m A levels are essential for cell cycle control. Moreover, we marks. We cannot rule out the possibility that the m A found that a higher number of RNAs were marks located in 5′ UTR may be misclassified into start hypermethylated in Alkbh5-deficient cerebella of P7, which codon regions for those poly(A) RNAs with a very short were originally methylated only in fully differentiated, 5′ UTR. Thereby, single nucleotide-resolution m A mature neural cells in WT cerebellum. Such untimely and analysis [48, 52] should be utilized next to identify the aberrant methylation finally appears to result in premature exact positions of m A marks, which will further enable differentiation and aberrant neurogenesis in the exploration of the functions and regulatory mechanism cerebellum. This may explain the phenotypes of lower of m A deposited in different regions of the transcripts. numbers of mature granule neurons in IGLs, together with In addition, the different patterns of m A observed disorganized Purkinje cells, and astrocytes. Generally, cell between cell lines and tissues indicates the complexity of fate determination requires a burst of expression of a in vivo temporal regulation of RNA methylation, distinct group of genes at each specified stage, which is suggesting the necessity to characterize RNA spatially and temporally regulated throughout the whole methylation in a context-dependent manner. Ma et al. Genome Biology (2018) 19:68 Page 12 of 18 Gene regulation in the brain is characterized by its con- Antibodies siderable cellular heterogeneity, including gene expression Detailed information on the primary and secondary anti- [53], DNA hydroxymethylation [34], and m ARNA bodies used in this study is listed in Additional file 1: methylation [36]. Accordingly, we found that during Table S9. cerebellar development, the expression of m Awriter and eraser genes in granule neurons and Purkinje neurons RNA isolation, poly(A) RNA purification, and m A-IP changed in opposite directions. Although the poly(A) RNAs Mice were sacrificed by cervical dislocation and cerebel- analyzed in this study were isolated from a mixture of all lum was removed and immediately snap-frozen in liquid neural cell types in the mouse cerebellum, the methylation nitrogen. Total RNA was extracted using TriReagent profiles depicted here mainly reflect the nature of m A (Sigma, T9424) according to the manufacturer’s instruc- modification in granule neurons as they constitute about tions. In general, four mice of P7, four mice of P14, and 90% of the cerebellar neural cells [35]. For further insights, three mice of P21 or P60 were used for one m A-IP it will be important to elucidate the methylation profiles reaction according to the established protocol [54]. and related biological functions in Purkinje cells and glial Briefly, poly(A) RNA was purified from total RNA using cells, respectively. A more detailed, cell type-specific methy- Poly(A)Purist™ MAG Kit (Thermo Fisher, AM1922). lation analysis will be necessary to gain insights into the bio- Next, poly(A) RNA was fragmented to 100-150 nucleotides 6 6 logical functions and developmental regulation of m A (for m A-seq) using RNA Fragmentation Reagents modification in individual types of neural cells. (Thermo Fisher, AM8740) as instructed by the manufacture. Fragmented poly(A) RNA (1 μg) was Conclusions incubated with 2.5 μgofanti-N -methyl-adenosine Here, we describe the first m A RNA methylation antibody at 4 °C for 4 h, followed by addition of landscape during postnatal development of the mouse Dynabeads™ Protein A for Immunoprecipitation (Thermo cerebellum and reveal the unique features of continuously Fisher, 10002D). After incubation at 4 °C for 2 h, beads and temporally specific methylation in the mouse were extensively washed five times with IPP buffer cerebellum from P7 to P60. Our results indicate that (150 mM NaCl, 0.1% NP-40, 10 mM Tris-HCl, pH 7.4); temporal regulation of RNA m A methylation, orchestrated immunoprecipitated RNA was recovered through elution 6 6 by spatiotemporal expression of m A writers and erasers, is with m A nucleotide followed by ethanol precipitation. essential in regulation of postnatal cerebellar development. For m A-seq, a cDNA library was constructed using a We also provide the first in vivo evidence that ALKBH5- KAPA Stranded mRNA-seq kit (Kapa Biosystems, induced abnormal RNA methylation affects RNA nuclear KR0960) and subjected to next-generation sequencing on export. Further studies will be necessary to uncover the de- an Illumina Hiseq X Ten system. Input RNA before im- tailed molecular mechanism of the regulation and bio- munoprecipitation was applied to RNA-seq in parallel. 6 6 6 logical functions of m A during mouse cerebellar For gene-specific m A-IP-qPCR, m A-IP was development. In addition, characterization of the role of performed using 3 μg of fragmented poly(A) RNA (300– m A in the cerebellum will offer a new viewpoint to 500 nucleotides). Reverse transcription was carried out elucidate the mechanism of related neurological diseases. with an equal ratio of RNA from input and IP product by using ReverTra Ace® qPCR RT Master Mix Methods (TOYOBO, FSQ-301). Quantitative real-time PCR was Animal care performed by using THUNDERBIRD™ SYBR® qPCR Mix Wild-type C57BL/6 mice were purchased from Vital River (TOYOBO, QPS-201). Gapdh was used as a negative company (Beijing, China). The Alkbh5-deficient mice were control. The primers used in this study are listed in used and genotyped as described previously [9]. For all ex- Additional file 1: Table S10. periments performed here, equal numbers of male and fe- male mice were included for analysis. All animal Protein extraction and western blot analysis experiments and euthanasia were approved and performed Cerebella were dissected from mice of different ages and in accordance with the guidelines of Animal Care and Use triturated in RIPA lysis buffer with protease inhibitor Committee of IBMS/PUMC. The IRB (Institutional Review cocktail (Roche, 04693159001) and phosphatase inhibi- Board) approval number is ACUC-A02–2014-001. tors. Tissue lysate (60–80 μg) was subjected to 10% SDS-PAGE. After electrotransfer, the blots were blocked Cell lines in 5% milk in Tris-buffered saline with Tween 20 (TBST) Neuro2a cells and HEK293T cells were purchased from for 1 h at room temperature, then incubated at 4 °C National Infrastructure of Cell Line Resource (Beijing, overnight with the primary antibodies. The transferred China). The cell lines have been authenticated and tested membranes were incubated with appropriate HRP- for mycoplasma contamination by the provider. conjugated secondary antibodies, exposed to enhanced Ma et al. Genome Biology (2018) 19:68 Page 13 of 18 chemiluminescence (ECL) solution, and visualized by gel VSV-G (kindly provided by Prof. Qi Xu). shRNA of image analysis system (Tanon, 5800). mouse Mettl3 was cloned into pLL3.7 vector between HpaI and XhoI restriction enzyme sites. Two sets of Histology and immunostaining analysis shRNAs were used for each gene and were evaluated for After anesthesia with saturated tribromoethanol solu- their knockdown efficiency in Neuro2a cells via western tion, mice were perfused with 4% paraformaldehyde; the blot analysis. The shRNA with higher efficiency was whole brains were dissected and post-fixed in 4% para- chosen for subsequent in vivo virus infection of mouse formaldehyde for 24–48 h at 4 °C. For paraffin sections, cerebellum. The oligos of Mettl3 shRNA were: 1, 5′-aacc the cerebellum was dehydrated with ethanol, clarified GCACACTGATGAATCTTTAGGTtcaagag ACCTAAA- with xylene, and then embedded in paraffin. We used 4- GATTCATCAGTGTGCtttttc-3′;2,5′-aaccCTGCAAA- μm-thick sections for hematoxylin and eosin (H&E) TATGTTCACTATGAtcaagag TCATAGTGAACATATT staining or immunohistochemistry analysis. For cryosec- TGCAGtttttc-3′. The negative control containing tions, the cerebella were dehydrated in 30% sucrose for scrambled shRNA was 5′-aaccTTCTCCGAACGTGT- more than 24 h, embedded in OCT compound and fro- CACGTtcaagagACGTGACACGTTCGGAGAAtttttc-3′. zen at − 40 °C. We used 12-μm-thick sections for im- The plasmids generated here were written in short as munofluorescence analysis. Immunohistochemistry and OE, sh-1, sh-2. immunofluorescence staining were performed according to standard protocols [55]. Lentivirus infection of mouse cerebellum For BrdU incorporation analysis, the mice were intra- All lentivirus vectors were purified and amplified using peritoneally injected with BrdU (Sigma, B5002; final QIAGEN Plasmid maxi kit (Qiagen, number 12163). 293 T concentration 50 μg/g) and analyzed 2 hours after BrdU cells were co-transfected with a mixture of lentiviral plas- injection. Immunofluorescence analysis was performed mids (including 17.5 μg of pLV-derived plasmid, 13.125 μg as described above by using anti-BrdU antibody. of pH 1, 4.375 μgofpH2foroverexpressionsystem; 14 μg Images of H&E and immunohistochemical staining were of pLL3.7-derived plasmid, 7 μg of pMDLg/pRRE, 7 μgof recorded using a Panoramic MIDI II digital slide scanner pCMV-VSV-G, 7 μg of pRSV-Rev for knockdown system) (3D Histech). Results of immunofluorescence analysis using DNA transfection reagent (Neofect, TF20121201). were imaged with a confocal microscope (Zeiss, LSM780). Ssupernatants were collected 48 h after transfection and concentrated by ultracentrifugation. Titration of lentivirus Quantitative m A level measurement using UHPLC-MS/MS was carried out by transducing 293 T cells in tenfold serial The absolute amount of m A in the poly(A) RNA was dilution. Generally the lentivirususedinthisstudy hada 8 10 measured as described before [36]. Briefly, 100 ng of titer of 1 × 10 to 1 × 10 TU/ml. poly(A) RNA was digested with nuclease P1 (Sigma) at In vivo lentivirus infection was performed as previ- 37 °C for 2 h, followed by treatment with calf intestine ously described with minor modification [56]. Briefly, alkaline phosphatase (CIAP, Promega) at 37 °C for 2 h. wild-type C57BL/6 mice (P7) were anesthetized on ice The solution was filtered and analyzed using ultra- and then injected with a mixture of lentivirus (2 μl), performance liquid chromatography and a triple-4 quad- polybrene (800 ng/μl), and fastgreen (0.1%) into the rupole mass spectrometer (AB SCIEX QTRAP 5500). interlobular space between lobules V and VI. Pups were m A level in each sample was calculated by comparing revived at 37 °C and then returned to original cages. with the standard curve obtained from pure nucleoside After 7 (for overexpression) or 10 (for knockdown) days standards loaded simultaneously. The ratio of m AtoA post-lentivirus infection, immunofluorescence analysis was calculated to obtain the global methylation level. was performed as described above to detect the pheno- types of cerebellum upon lentivirus infection. Plasmid construction The lentiviral overexpression vectors used in this study Hypobaric hypoxia treatment of mice include pLV-EF1a-EGFP(2A) Puro, pH 1, and pH 2 vec- WT or Alkbh5 KO mice of P5, together with their mothers, tors (kindly provided by Prof. Kaili Ma). cDNA of mouse were placed into a hypobaric hypoxic chamber and exposed Mettl3 was cloned into pLV-EF1a-EGFP(2A) Puro vector to a simulated atmospheric pressure of 10.6 kPa digested with EcoR I and BamH I restriction enzymes. (377 mmHg), which was equivalent to an altitude of The overexpression efficiency of lentivirus was verified 5500 m. All animals were kept at constant temperature in Neuro2a cells by western blotting. The primers used (25–30 °C) on a daily light schedule of 12 h of light vs. dark were: forward, 5′-cggaattc ATGTCGGACACGTGGAG- with normal activity. The hypobaric hypoxic condition was 3′; reverse, 5′- cgggatccCGCTATAAATTCTTAGG-3′. maintained and monitored continuously with the sensor in- −/− The lentiviral knockdown system used in this study in- side the chamber. Wild-type and Alkbh5 pups were cludes pLL3.7, pRSV-rev, pMDLg/pRRE, and pCMV- dissected for phenotype analysis after incubation in the Ma et al. Genome Biology (2018) 19:68 Page 14 of 18 hypobaric chamber for 48 h. To examine the phenotypes of two mismatched bases; 2) bases with low quality score hypoxia-treated mice at later developmental stages, the WT (< 20) were trimmed off from the 3′ end; 3) reads with and KO mice were returned to normoxic conditions to- length longer than 70 nucleotides and more than 70% gether with their mothers and raised for another 7 days be- bases with quality score > 25 were retained. These high- fore dissection for phenotype analysis. quality reads were mapped against the mouse genome (mm10) allowing up to two mismatches using TopHat Magnetic resonance imaging analysis software (version 2.0.13) [57, 58]. Only uniquely mapped Mice were anaesthetized with 0.75% amobarbital before reads were kept for the subsequent analysis. analysis. Mouse brain MRI was performed on a 7.0 T BioClinscan Animal MRI System (Bruker, Germany) with RNA expression analysis Siemens manipulating software. T2-weighted images of The FPKM (fragments per kilobase of transcript per million the whole brain (including the sagittal, coronal, and trans- mapped reads) values of RNAs in each sample were calcu- verse views) were acquired with a 2D-TSE (2D-turbo lated using Cufflinks toolkit (version 2.0.2) [59]. Only the spin-echo) sequence. The scanning time for each mouse transcripts with FPKM value > 0.2 were considered as was about 10 min. The scanning parameters were as fol- expressed transcripts [60]. The differentially expressed lows: sagittal view, TR = 2360 ms, TE = 41 ms, number of RNAs (DERs, P < 0.05) between any two samples were averages = 1, 15 axial slices with a slice thickness (ST) of 0. identified using Cuffdiff. The normalized FPKM value for 7 mm, a field of view (FOV) of 3.52 cm × 3.52 cm, and a each sample was used for clustering using the KMC matrix of 320 × 320 were positioned over the whole brain method of MeV software with the following parameter set- with a pixel spacing of 0.11 mm; coronal view, TR = tings: 1) Pearson correlation clustering; 2) K-means cluster- 3510 ms, TE = 53 ms, number of averages = 1, axial slices ing for six clusters; 3) other parameters set as defaults [61]. = 19, ST = 0.5 mm, image matrix size = 320 × 384, pixel To verify the reliability of the DERs, we performed a par- spacing = 0.078 mm, FOV 2.50 cm × 3.00 cm; transverse allel differential expression analysis by STAR/edgeR pack- view, TR = 3380 ms, TE = 41 ms, number of averages = 1, ages. High quality reads in each sample were mapped to axial slices = 21, ST = 0.7 mm, image matrix size = 320 × the mouse genome (mm10) using STAR software and the 320, pixel spacing = 0.094 mm, FOV 2.82 cm × 2.82 cm. uniquely mapped reads were kept for the subsequent ana- Then the signals of different brain regions were measured lysis [62, 63]. HTSeq software [64] was utilized to calculate and compared using a RadiAnt DICOM Viewer. the read count of RNAs and edgeR [65] was applied to identify the DERs (P < 0.05). All software applications were Subcellular fractionation analysis of RNA run with default parameters. Fresh cerebella of WT and KO mice at P7 were minced and incubated in Digestion Solution (30 U/mL papain, 240 μg/ m A peak calling and motif analysis 6 6 mL cystein, 400 μg/mL DNAseItype IV)at37°Cfor 1h. RNA m A-modified regions, also called m A peaks, were The reaction was stopped by adding Ovomuccoid inhibitor identified using exomePeak software (version 2.7.0) with solution (1125 μg/mL Ovomuccoid trypsine inhibitor, FDR (false discovery rate) < 0.05 [66, 67]. The consensus 525 μg/mL BSA, 400 μg/mL DNase I type IV) and incuba- sequence motifs enriched in m A peaks were identified by tion at 37 °C for 4 min. After purification, cerebellar cells HOMER [68]. For further comparison of m A modification were subjected to cellular fractionation using a PARIS™ Kit between samples, we used coverageBed of BEDToods (Thermo Fisher, AM1921) according to the manufacturer’s (version 2.26.0) [69]withthe “-s”, “-splited”, “-counts”,and protocol. Nuclear and cytoplasmic RNAs (120 ng) were re- “-F 0.50” parameters to calculate the read count of each verse transcribed using ReverTra Ace® qPCR RT Master peak. The “IP FPKM”, “input FPKM”,and “Enrichment Mix (TOYOBO, FSQ-301). Quantitative real-time PCR was score” of peaks were calculated as following methods: performed using THUNDERBIRD™ SYBR® qPCR Mix a ¼ A x10 = B xC ð1Þ i; j i; j j i; j (TOYOBO, QPS-201). Nd1 was used as a cytoplasmic marker while U1 was used as a nuclear marker. The primers b ¼ D x10 = E xC ð2Þ i; j i; j j i; j used in this study are listed in Additional file 1:Table S10. c ¼ a =b ð3Þ i; j i; j i; j RNA-seq data processing and reads mapped For each sample, single-end reads were used for bio- where a denotes the peak “IP FPKM” value of the ith i,j informatics analysis. Quality control of raw data was methylation peak in the IP sample from the jth bio- done using FastQC software (version 0.10.1). Sequencing logical sample; A denotes the total reads mapped to the i,j data were preprocessed with in-house Perl scripts using ith methylation peak in the IP sample from the jth bio- the following criteria: 1) adaptor sequence was removed logical sample; B denotes the total unique reads mapped by finding the sequence AGATCGGAAG with at most to the mouse reference (mm10) in the IP sample from Ma et al. Genome Biology (2018) 19:68 Page 15 of 18 the jth biological sample; C denotes the length (base) bin was calculated to represent the occupancy of m A i,j of the ith methylation peak in the IP sample from the peaks along the whole transcripts. jth biological sample; b denotes the peak “input FPKM” We further counted the number of m A peaks in the i,j value of the ith methylation peak in the input sample 5′ UTR, start codon region, CDS, stop codon region and from the jth biological sample; D denotes the total 3′ UTR, respectively. Specifically, a 300-nucleotide i,j reads mapped to the ith methylation peak in the input region centered on start codons or stop codons was de- sample from the jth biological sample; E denotes the fined as start codon regions or stop codon regions [36]. total unique reads mapped to the mouse reference (mm10) in the input sample from the jth biological sam- 6 Correlation analysis between RNA m A methylation level ple; c denotes the peak “enrichment score” value of the i,j and RNA expression level ith methylation peak in the IP/input sample from the jth The global clustering analysis of RNA methylation levels biological sample. (enrichment score) and expression levels (FPKM) among M A peaks that satisfied 1) exomePeak FDR value < 0. the four stages was performed using Cluster 3.0 soft- 6 6 05, 2) m A peak “IP FPKM” value > 1, and 3) m A peak ware. The enrichment scores and FPKM of each RNA in “enrichment score” value > 1.5 were used for further each sample were adjusted with “mean” parameters in comparative analysis. The continuously methylated “Normalize genes” and “Center genes”. The hierarchical RNAs (CMRs) were defined as RNAs containing at least clustering method was employed using “Complete link- one m A peak in all the four samples, while the age”. The clustered heatmap figure was generated using temporal-specific methylated RNAs (SMRs) were the TreeView-1.1.6r4-win software. RNAs with m A modification only in one sample, and In order to obtain a dynamic view of the regulatory effect 6 6 m A peaks of SMRs were defined as specific m A peaks. 6 of RNA m A methylation on RNA abundance, we By merging all the m A peaks of four developmental calculated the correlation coefficients and the associated P stages using mergeBed of BEDTools with “-s”, “-c”, and 6 valuebetween theenrichment scoresofm A peaks and “-o” parameters, we identified “ON” or “OFF” RNA m A FPKM of methylated RNAs at the four stages using switches during cerebellar development. All m A peaks Pearson test in R package. High-confidence regulatory gene absent in the former stage but present in the later stage pairs were defined by the correlation coefficients (|r| > 0.95) 6 6 were called “ON” m A switches, while all m A peaks with P value < 0.05. Cluster analysis was performed for displaying changes in the opposite direction were called those positively (r > 0.95) and negatively (r < − 0.95) corre- “OFF” m A switches. In addition, the differentially lated RNAs using MeV software, respectively. The MeV methylated RNAs (DMRs) between wild-type and cluster results were integrated and plotted using R pro- Alkbh5 knockout mice were further identified using exo- gramming language. mePeak software with default parameters, and high- confidence DMR m A peaks matching the criteria of GO analysis FDR < 0.05 and log2 (fold change) > 1 or < − 1 were The DAVID tool [74] was used for GO analysis by apply- selected for further analysis. ing default parameters, except that only those transcripts For validation analysis of the distribution patterns of 6 fulfilling the condition of FPKM > 0.2 in the input samples stage-specific m A peaks which were obtained by using were set as background. Either bubble plots or column the exomePeak software, peak calling analysis was plots were generated based on the enriched GO terms performed in parallel using the MACS2 software for the using the ggplot2 in R package [75]. Color intensity indi- P7 and P60 samples with the default parameters [54, 70]. cates the value of −log (P value), while the length of col- The integrative genomics viewer (IGV) tool was used 6 umns or the size of the circles indicate the gene counts. A for visualization of m A peaks along the whole transcript 6 full list of all selected terms of the biological process, [71, 72]. The heatmaps of the enrichment score for m A cellular components, and molecular functions category peaks were plotted using the pheatmap package in R [73]. are provided in Supplemental tables (Additional file 2-4: TableS4-S6; Additionalfile 6: Table S8). Characterization of m A peak distribution patterns Distribution of m A peaks along mRNAs were Kyoto Encyclopedia of Genes and Genomes analysis determined as previously described with minor revision All modified genes were annotated to Kyoto Encyclopedia [3, 18]. To characterize the distribution patterns of m A of Genes and Genomes (KEGG) pathways using KAAS peaks, a reference mouse transcriptome was built using (KEGG Automatic Annotation Server) [76], and the the longest transcript of each gene. Next, each of the 5′ enriched KEGG pathways for the genes encoded by P7 or UTR, CDS, and 3′ UTR regions were split into 100 bins P60 temporal-specific SMRs or CMRs were identified with equal length. The percentage of m A peaks in each using Fisher’s exact test (Additional file 4: Table S6) [77]. Ma et al. Genome Biology (2018) 19:68 Page 16 of 18 6 6 Quantification and statistical analysis Genomes; KO: Alkbh5-knockout; m A: N -methyladenosine; P7: Postnatal day 7; PCL: Purkinje cell layer; PH3: Phospho-histone 3; SMRs: Specifically The differences in distribution for m A enrichment score methylated RNAs; UHPLC-MS/MS: Ultra high performance liquid and log FPKM between samples were detected by chromatography-mass spectrum/mass spectrum analysis; TBST: Tris Buffered Wilcoxon test. Pearson test was used to perform Saline with Tween 20; WT: Wild-type correlation analysis. All statistical analysis and graphs of Acknowledgements results in Figs. 3d, 4d, 5a, 5c–5f and 6f, and Additional file We thank Prof. Xu Q and Dr. Ma K-L for their generosity in offering lentivirus 1: Figures S5c, S6a, S7d, S7e, S8a and S8g were assessed vectors. using two-tailed unpaired Student’s t-test and performed Funding using GraphPad Prism 6.0 software. Results are presented This work was supported by the National Natural Science Foundation of China as mean ± s.e.m. Immunoreactive cells in 3–5 randomly se- (31471288 to NY, 31471343 to TWM), CAMS Initiative for Innovative Medicine (2016-I2M-1-004 to NY, 2016-I2M-2-001 to TWM), and The Youth Innovation Pro- lected lobules of one cerebellum were counted using Stata- motion Association of Chinese Academy of Science (2017141 to SS.) Quest 5.0 software. For each subject, data were collected from five to seven mice. For each mouse, data were ob- Availability of data and materials The datasets have been deposited in the Sequenced Read Archive (SRA) tained from the immunostaining results of two to three under accession number PRJNA400297 (SRX3136161-SRX3136172) [78], and near-midline slices. also in the Genome Sequence Archive [79] in BIG Data Center [80], Beijing Institute of Genomics (BIG), Chinese Academy of Science under accession number CRA000472 [81]. Additional files Authors’ contributions NY and TWM conceived the project and, with SS and KA, designed the data Additional file 1: Figure S1. Dynamic RNA methylation in the analysis. MC performed phenotype analysis. CM performed m A-IP and developing mouse cerebellum. Figure S2. Comparison between phenotype analysis. LH, ZW and TX performed bioinformatic analyses. ZZW continuous and temporal-specific methylation during mouse cerebellar performed lentivirus infection experiments. HX performed mouse brain development. Figure S3. Comparison of RNA methylation between P7 dissection experiments. WG performed plasmid construction experiments. ZS and P60 based on m A peaks identified using MACS2. Figure S4. and ZY assisted with immunostaining analyses. WD assisted with western Correlation between RNA methylation and gene expression during blot analyses. LC and LQ performed MRI analysis. NY, MC and LH wrote the cerebellar development. Figure S5. Morphology analysis of mouse manuscript with input from all authors. All authors read and approved the cerebellum upon lentivirus infection for Mettl3 overexpression. Figure S6. final manuscript. Phenotype analysis of the cerebellum in the WT and KO mice under normoxic condition. Figure S7. Morphology analysis of the cerebellum in Ethics approval WT and KO mice exposed to hypobaric hypoxia and normoxia All animal experiments and euthanasia were approved and performed in successively. Figure S8. Dysregulated RNA methylation resulting from accordance with the guidelines of Animal Care and Use Committee of IBMS/ Alkbh5 deficiency in mouse cerebellum exposed to hypobaric hypoxia. PUMC. The IRB (Institutional Review Board) approval number is ACUC-A02– Table S1. Data quality and processing information of m A-seq of poly(A) 2014-001. RNA from wild-type mouse cerebellum (P7, P14, P21, and P60), the cere- bellum of wild-type (WT) and Alkbh5 knockout (KO) mice exposed to Competing interests hypobaric hypoxia (P7). Table S2. Statistics of m A peaks and expressed The authors declare that they have no competing interests. RNAs in wild-type mouse cerebellum (P7, P14, P21, and P60), the cerebel- lum of wild-type (WT) and Alkbh5 knockout (KO) mice exposed to hypo- baric hypoxia (P7). Table S3. Numbers of m A peaks located in different Publisher’sNote regions of mRNA transcripts in wild-type mouse cerebellum at P7, P14, Springer Nature remains neutral with regard to jurisdictional claims in P21, and P60. Table S9. List of antibodies and their applications used in published maps and institutional affiliations. this study. Table S10. List of primers for RT-qPCR used in this study. (PDF 14643 kb) Author details Additional file 2: Table S4. GO analysis of genes containing m AON Department of Pathology, Institute of Basic Medical Sciences Chinese and OFF switches during mouse cerebellar development. (XLSX 564 kb) Academy of Medical Science, School of Basic Medicine Peking Union Medical College; Neuroscience Center, Chinese Academy of Medical Sciences, Beijing Additional file 3: Table S5. GO analysis of genes encoded by the CMRs 100005, China. BIG Data Center, Beijing Institute of Genomics, Chinese at P7 and P60, and SMRs at the four developmental stages. (XLSX 128 kb) Academy of Sciences, Beijing 100101, China. University of Chinese Academy Additional file 4: Table S6. GO and KEGG pathway enrichment of Sciences, People’s Republic of China, Beijing 100049, China. State Key analysis of the genes encoded by SMRs and CMRs compared between Laboratory of Molecular Oncology, National Cancer Center/Cancer Institute the m A peaks at P7 and P60 that were identified by using MACS2 and Hospital, Chinese Academy of Medical Sciences, Peking Union Medical software. (XLSX 303 kb) College, Beijing 100021, China. Department of Microbiology, Oslo University Additional file 5: Table S7 List of the 839 RNAs with strong positive or Hospital, Rikshospitalet, Oslo NO-0027, Norway. Department of Molecular negative correlation between their methylation levels and expression Medicine, Institute of Basic Medical Sciences, University of Oslo, NO-0317 levels. (XLSX 82 kb) Oslo, Norway. Additional file 6: Table S8 GO analysis of the genes with altered m A Received: 3 December 2017 Accepted: 26 April 2018 levels between wild-type and Alkbh5-deficient mice exposed to hypobaric hypoxia. (XLSX 83 kb) References Abbreviations 1. Yao B, Christian KM, He C, Jin P, Ming GL, Song H. Epigenetic mechanisms BP: Biological processes; CC: Cellular components; CDS: Coding sequence; in neurogenesis. Nat Rev Neurosci. 2016;17:537–49. CMRs: Continuously methylated RNAs; DMRs: Differentially methylated RNAs; 2. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar ECL: Enhanced Chemiluminescence; EGL: External granule cell layer; L, Osenberg S, Cesarkas K, Jacob-Hirsch J, Amariglio N, Kupiec M, et al. FDR: False discovery rate; GO: Gene Ontology; HE: Hematoxylin-eosin Topology of the human and mouse m6A RNA methylomes revealed by staining; IGL: Inner granule cell layer; KEGG: Kyoto Encyclopedia of Genes and m6A-seq. Nature. 2012;485:201–6. Ma et al. Genome Biology (2018) 19:68 Page 17 of 18 3. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. 26. Li HB, Tong J, Zhu S, Batista PJ, Duffy EE, Zhao J, Bailis W, Cao G, Comprehensive analysis of mRNA methylation reveals enrichment in 3' Kroehling L, Chen Y, et al. m6A mRNA methylation controls T cell UTRs and near stop codons. Cell. 2012;149:1635–46. homeostasis by targeting the IL-7/STAT5/SOCS pathways. Nature. 2017; 4. Niu Y, Zhao X, Wu YS, Li MM, Wang XJ, Yang YG. N6-methyl-adenosine 548:338–42. (m6A) in RNA: an old modification with a novel epigenetic function. 27. Hess ME, Hess S, Meyer KD, Verhagen LA, Koch L, Bronneke HS, Dietrich MO, Genomics Proteomics Bioinformatics. 2013;11:8–17. Jordan SD, Saletore Y, Elemento O, et al. The fat mass and obesity associated gene (Fto) regulates activity of the dopaminergic midbrain 5. Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2017;18:31–42. circuitry. Nat Neurosci. 2013;16:1042–8. 6. Cao G, Li HB, Yin Z, Flavell RA. Recent advances in dynamic m6A RNA 28. Li L, Zang L, Zhang F, Chen J, Shen H, Shu L, Liang F, Feng C, Chen D, Tao modification. Open Biol. 2016;6:160003 H, et al. Fat mass and obesity-associated (FTO) protein regulates adult 7. Zhao X, Yang Y, Sun BF, Shi Y, Yang X, Xiao W, Hao YJ, Ping XL, Chen YS, neurogenesis. Hum Mol Genet. 2017;26:2398-2411. Wang WJ, et al. FTO-dependent demethylation of N6-methyladenosine 29. Walters BJ, Mercaldo V, Gillon CJ, Yip M, Neve RL, Boyce FM, Frankland PW, regulates mRNA splicing and is required for adipogenesis. Cell Res. 2014; Josselyn SA. The role of the RNA demethylase FTO (Fat mass and obesity- 24:1403–19. associated) and mRNA methylation in hippocampal memory formation. 8. Xiao W, Adhikari S, Dahal U, Chen YS, Hao YJ, Sun BF, Sun HY, Li A, Ping XL, Neuropsychopharmacology. 2017;42:1502–10. Lai WY, et al. Nuclear m(6)A Reader YTHDC1 Regulates mRNA Splicing. Mol 30. Widagdo J, Zhao QY, Kempen MJ, Tan MC, Ratnu VS, Wei W, Leighton L, Cell. 2016;61:507–19. Spadaro PA, Edson J, Anggono V, Bredy TW. Experience-dependent accumulation of N6-methyladenosine in the prefrontal cortex is associated 9. Zheng G, Dahl JA, Niu Y, Fedorcsak P, Huang CM, Li CJ, Vagbo CB, Shi Y, with memory processes in mice. J Neurosci. 2016;36:6771–7. Wang WL, Song SH, et al. ALKBH5 is a mammalian RNA demethylase that 31. Yoon K-J, Ringeling FR, Vissers C, Jacob F, Pokrass M, Jimenez-Cyrus D, Su Y, impacts RNA metabolism and mouse fertility. Mol Cell. 2013;49:18–29. Kim N-S, Zhu Y, Zheng L, et al. Temporal control of mammalian cortical 10. Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, Fu Y, Parisien M, Dai Q, Jia neurogenesis by m6A methylation. Cell. 2017;171:877-889.e17. G, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature. 2014;505:117–20. 32. Weng YL, Wang X, An R, Cassin J, Vissers C, Liu Y, Liu Y, Xu T, Wang X, 11. Ke S, Pandya-Jones A, Saito Y, Fak JJ, Vagbo CB, Geula S, Hanna JH, Black Wong SZH, et al. Epitranscriptomic m(6)A regulation of axon regeneration DL, Darnell JE Jr, Darnell RB. m6A mRNA modifications are deposited in in the adult mammalian nervous system. Neuron. 2018;97:313–25. e316 nascent pre-mRNA and are not required for splicing but do specify 33. Goldowitz D, Hamre K. The cells and molecules that make a cerebellum. cytoplasmic turnover. Genes Dev. 2017;31:990–1006. Trends Neurosci. 1998;21:375–82. 12. Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, Weng X, Chen K, Shi H, 34. Szulwach KE, Li X, Li Y, Song CX, Wu H, Dai Q, Irier H, Upadhyay AK, Gearing He C. N(6)-methyladenosine modulates messenger RNA translation M, Levey AI, et al. 5-hmC-mediated epigenetic dynamics during postnatal efficiency. Cell. 2015;161:1388–99. neurodevelopment and aging. Nat Neurosci. 2011;14:1607–16. 13. Zhou J, Wan J, Gao X, Zhang X, Jaffrey SR, Qian SB. Dynamic m(6)A mRNA 35. Frank CL, Liu F, Wijayatunge R, Song L, Biegler MT, Yang MG, Vockley CM, methylation directs translational control of heat shock response. Nature. Safi A, Gersbach CA, Crawford GE, West AE. Regulation of chromatin 2015;526:591–4. accessibility and Zic binding at enhancers in the developing cerebellum. 14. Coots RA, Liu X-M, Mao Y, Dong L, Zhou J, Wan J, Zhang X, Qian S-B. Nat Neurosci. 2015;18:647–56. m6A facilitates eIF4F-independent mRNA translation. Mol Cell. 2017;68: 36. Chang M, Lv H, Zhang W, Ma C, He X, Zhao S, Zhang Z-W, Zeng Y-X, Song 504–14. e507 S, Niu Y, Tong W-M. Region-specific RNA m A methylation represents a 15. Xiang Y, Laurent B, Hsu CH, Nachtergaele S, Lu Z, Sheng W, Xu C, Chen H, new layer of control in the gene regulatory network in the mouse brain. Ouyang J, Wang S, et al. RNA m6A methylation regulates the ultraviolet- Open Biol. 2017;7:170166 induced DNA damage response. Nature. 2017;543:573–6. 37. Barbaric I, Miller G, Dear TN. Appearances can be deceiving: phenotypes of 16. Batista PJ, Molinie B, Wang J, Qu K, Zhang J, Li L, Bouley DM, Lujan E, Haddad knockout mice. Brief Funct Genomic Proteomic. 2007;6:91–103. B, Daneshvar K, et al. m(6)A RNA modification controls cell fate transition in 38. Bunn HF, Poyton RO. Oxygen sensing and molecular adaptation to hypoxia. mammalian embryonic stem cells. Cell Stem Cell. 2014;15:707–19. Physiol Rev. 1996;76:839–85. 17. Aguilo F, Zhang F, Sancho A, Fidalgo M, Di Cecilia S, Vashisht A, Lee DF, 39. Zhang C, Samanta D, Lu H, Bullen JW, Zhang H, Chen I, He X, Semenza GL. Chen CH, Rengasamy M, Andino B, et al. Coordination of m(6)A mRNA Hypoxia induces the breast cancer stem cell phenotype by HIF-dependent methylation and gene transcription by ZFP217 regulates pluripotency and and ALKBH5-mediated m6A-demethylation of NANOG mRNA. Proc Natl reprogramming. Cell Stem Cell. 2015;17:689–704. Acad Sci U S A. 2016;113:E2047–56. 18. Chen T, Hao YJ, Zhang Y, Li MM, Wang M, Han W, Wu Y, Lv Y, Hao J, Wang 40. Jakobsen L, Vanselow K, Skogs M, Toyoda Y, Lundberg E, Poser I, Falkenby L, et al. m(6)A RNA methylation is regulated by microRNAs and promotes LG, Bennetzen M, Westendorf J, Nigg EA, et al. Novel asymmetrically reprogramming to pluripotency. Cell Stem Cell. 2015;16:289–301. localizing components of human centrosomes identified by complementary 19. Zhang C, Chen Y, Sun B, Wang L, Yang Y, Ma D, Lv J, Heng J, Ding Y, Xue Y, proteomics methods. EMBO J. 2011;30:1520–35. et al. m6A modulates haematopoietic stem and progenitor cell 41. Sarzi E, Angebault C, Seveno M, Gueguen N, Chaix B, Bielicki G, Boddaert N, specification. Nature. 2017;549:273–6. Mausset-Bonnefont AL, Cazevieille C, Rigau V, et al. The human OPA1delTTAG mutation induces premature age-related systemic 20. Zheng Q, Hou J, Zhou Y, Li Z, Cao X. The RNA helicase DDX46 inhibits neurodegeneration in mouse. Brain. 2012;135:3599–613. innate immunity by entrapping m6A-demethylated antiviral transcripts in the nucleus. Nat Immunol. 2017;18:1094–103. 42. Cui C, Chatterjee B, Lozito TP, Zhang Z, Francis RJ, Yagi H, Swanhart LM, 21. Hsu PJ, Zhu Y, Ma H, Guo Y, Shi X, Liu Y, Qi M, Lu Z, Shi H, Wang J, et al. Sanker S, Francis D, Yu Q, et al. Wdpcp, a PCP protein required for Ythdc2 is an N6-methyladenosine binding protein that regulates ciliogenesis, regulates directional cell migration and cell polarity by direct mammalian spermatogenesis. Cell Res. 2017;27:1115–27. modulation of the actin cytoskeleton. PLoS Biol. 2013;11:e1001720. 22. Lin Z, Hsu PJ, Xing X, Fang J, Lu Z, Zou Q, Zhang KJ, Zhang X, Zhou Y, 43. Jiang D, Zhao L, Clish CB, Clapham DE. Letm1, the mitochondrial Ca2+/H+ Zhang T, et al. Mettl3−/Mettl14-mediated mRNA N6-methyladenosine antiporter, is essential for normal glucose metabolism and alters brain modulates murine spermatogenesis. Cell Res. 2017;27:1216-1230. function in Wolf-Hirschhorn syndrome. Proc Natl Acad Sci U S A. 2013;110: 23. Xu K, Yang Y, Feng GH, Sun BF, Chen JQ, Li YF, Chen YS, Zhang XX, Wang E2249–54. CX, Jiang LY, et al. Mettl3-mediated m6A regulates spermatogonial 44. Sato A, Sekine Y, Saruta C, Nishibe H, Morita N, Sato Y, Sadakata T, differentiation and meiosis initiation. Cell Res. 2017;27:1100–14. Shinoda Y, Kojima T, Furuichi T. Cerebellar development transcriptome database (CDT-DB): profiling of spatio-temporal gene expression during 24. Zhao BS, Wang X, Beadell AV, Lu Z, Shi H, Kuuspalu A, Ho RK, He C. m6A- the postnatal development of mouse cerebellum. Neural Netw. 2008;21: dependent maternal mRNA clearance facilitates zebrafish maternal-to- 1056–69. zygotic transition. Nature. 2017;542:475–8. 25. Ivanova I, Much C, Di Giacomo M, Azzi C, Morgan M, Moreira PN, Monahan 45. Geula S, Moshitch-Moshkovitz S, Dominissini D, Mansour AA, Kol N, Salmon- J, Carrieri C, Enright AJ, O'Carroll D. The RNA m6A reader YTHDF2 is Divon M, Hershkovitz V, Peer E, Mor N, Manor YS, et al. Stem cells. m6A essential for the post-transcriptional regulation of the maternal mRNA methylation facilitates resolution of naive pluripotency toward transcriptome and oocyte competence. Mol Cell. 2017;67:1059–67. e1054 differentiation. Science. 2015;347:1002–6. Ma et al. Genome Biology (2018) 19:68 Page 18 of 18 46. He S, Wang H, Liu R, He M, Che T, Jin L, Deng L, Tian S, Li Y, Lu H, et al. 68. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, mRNA N6-methyladenosine methylation of postnatal liver development in Singh H, Glass CK. Simple combinations of lineage-determining pig. PLoS One. 2017;12:e0173421. transcription factors prime cis-regulatory elements required for macrophage 47. Zhang S, Zhao BS, Zhou A, Lin K, Zheng S, Lu Z, Chen Y, Sulman EP, Xie K, and B cell identities. Mol Cell. 2010;38:576–89. Bogler O, et al. m(6)A Demethylase ALKBH5 maintains tumorigenicity of 69. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing glioblastoma stem-like cells by sustaining FOXM1 expression and cell genomic features. Bioinformatics. 2010;26:841–2. proliferation program. Cancer Cell. 2017;31:591–606. e596 70. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, Liu XS. Model-based analysis of ChIP-Seq 48. Ke S, Alemu EA, Mertens C, Gantman EC, Fak JJ, Mele A, Haripal B, Zucker-Scharff (MACS). Genome Biol. 2008;9:R137. I, Moore MJ, Park CY, et al. A majority of m6A residues are in the last exons, 71. Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, allowing the potential for 3' UTR regulation. Genes Dev. 2015;29:2037–53. Mesirov JP. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6. 49. Luo GZ, MacQueen A, Zheng G, Duan H, Dore LC, Lu Z, Liu J, Chen K, Jia G, 72. Thorvaldsdottir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer Bergelson J, He C. Unique features of the m6A methylome in Arabidopsis (IGV): high-performance genomics data visualization and exploration. Brief thaliana. Nat Commun. 2014;5:5630. Bioinform. 2013;14:178–92. 50. Li Y, Wang X, Li C, Hu S, Yu J, Song S. Transcriptome-wide N(6)- 73. Kolde R: pheatmap. 2015. https://CRAN.R-project.org/package=pheatmap. methyladenosine profiling of rice callus and leaf reveals the presence of Published 11 Dec 2015. tissue-specific competitors involved in selective mRNA modification. RNA 74. Huang d W, Sherman BT, Lempicki RA. Systematic and integrative analysis of Biol. 2014;11:1180–8. large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57. 51. Tao X, Chen J, Jiang Y, Wei Y, Chen Y, Xu H, Zhu L, Tang G, Li M, Jiang A, 75. Wickham H: ggplot2: elegant graphics for data analysis. New York: Springer; et al. Transcriptome-wide N 6 -methyladenosine methylome profiling of porcine muscle and adipose tissues reveals a potential mechanism for 76. Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic transcriptional regulation and differential methylation pattern. BMC genome annotation and pathway reconstruction server. Nucleic Acids Res. Genomics. 2017;18:336. 2007;35:W182–5. 52. Linder B, Grozhik AV, Olarerin-George AO, Meydan C, Mason CE, Jaffrey SR. 77. Huang d W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: Single-nucleotide-resolution mapping of m6A and m6Am throughout the paths toward the comprehensive functional analysis of large gene lists. transcriptome. Nat Methods. 2015;12:767–72. Nucleic Acids Res. 2009;37:1–13. 53. LoVerso PR, Cui F. Cell type-specific transcriptome profiling in mammalian 78. Chunhui Ma, Mengqi Chang, Hongyi Lv, Yamei Niu, Shuhui Song, Wei-Min brains. Front Biosci (Landmark Ed). 2016;21:973–85. Tong, et al: RNA m6A methylation participates in regulation of postnatal 54. Dominissini D, Moshitch-Moshkovitz S, Salmon-Divon M, Amariglio N, development of the mouse cerebellum. SRA PRJNA400297. 2017. https:// Rechavi G. Transcriptome-wide mapping of N(6)-methyladenosine by www.ncbi.nlm.nih.gov/sra/?term=PRJNA400297. Accessed 25 Apr 2018. m(6)A-seq based on immunocapturing and massively parallel sequencing. 79. Wang Y, Song F, Zhu J, Zhang S, Yang Y, Chen T, Tang B, Dong L, Ding N, Nat Protoc. 2013;8:176–89. Zhang Q, et al. GSA: Genome Sequence Archive. Genomics Proteomics 55. Sanchez-Ortiz E, Cho W, Nazarenko I, Mo W, Chen J, Parada LF. NF1 Bioinformatics. 2017;15:14–8. regulation of RAS/ERK signaling is required for appropriate granule neuron 80. Members BIGDC. The BIG Data Center: from deposition to integration to progenitor expansion and migration in cerebellar development. Genes Dev. translation. Nucleic Acids Res. 2017;45:D18–24. 2014;28:2407–20. 81. Chunhui Ma, Mengqi Chang, Hongyi Lv, Yamei Niu, Shuhui Song, Wei-Min 56. Oomoto I, Suzuki-Hirano A, Umeshima H, Han YW, Yanagisawa H, Carlton P, Tong, et al: RNA m6A methylation participates in regulation of postnatal Harada Y, Kengaku M, Okamoto A, Shimogori T, Wang DO. ECHO-liveFISH: development of the mouse cerebellum. GSA CRA000472. 2017. http://bigd. in vivo RNA labeling reveals dynamic regulation of nuclear RNA foci in big.ac.cn/search?dbId=gsa&q=CRA000472. Accessed 25 Apr 2018. living tissues. Nucleic Acids Res. 2015;43:e126. 57. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R36. 58. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25:1105–11. 59. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7:562–78. 60. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5. 61. Saeed AI, Bhagabati NK, Braisted JC, Liang W, Sharov V, Howe EA, Li J, Thiagarajan M, White JA, Quackenbush J. TM4 microarray software suite. Methods Enzymol. 2006;411:134–93. 62. Dobin A, Gingeras TR. Optimizing RNA-Seq mapping with STAR. Methods Mol Biol. 2016;1415:245–62. 63. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21. 64. Anders S, Pyl PT, Huber W. HTSeq-a Python framework to work with high- throughput sequencing data. Bioinformatics. 2015;31:166–9. 65. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40. 66. Meng J, Cui X, Rao MK, Chen Y, Huang Y. Exome-based analysis for RNA epigenome sequencing data. Bioinformatics. 2013;29:1565–7. 67. Liu L, Zhang SW, Zhang YC, Liu H, Zhang L, Chen R, Huang Y, Meng J. Decomposition of RNA methylome reveals co-methylation patterns induced by latent enzymatic regulators of the epitranscriptome. Mol BioSyst. 2015;11: 262–74.

Journal

Genome BiologySpringer Journals

Published: May 31, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off