The role of m6A-RNA methylation in stress response regulation

1Introduction¶

Regulation of gene expression in response to stressful stimuli under healthy or pathological conditions involves epigenetic mechanisms such as DNA methylation and chromatin modifications Kloet et al., 2005McEwen et al., 2015. Elucidating the underlying molecular processes that regulate the fine-tuned transcriptional response to stress is essential for understanding stress vulnerability and the development of stress-related psychiatric disorders such as depression and anxiety.

In analogy to DNA modifications, a diverse set of covalent modifications is present on RNA nucleotides, encoding the epitranscriptome and regulating gene expression post-transcriptionally. These modifications determine how much protein is translated from available mRNA and how non-coding transcripts function Zhao et al., 2017. However, the role of this newly emerging layer of gene expression control in the central stress response and behaviour is not known yet. RNA modifications, next to epigenetic mechanisms, likely represent an yet undescribed level of transcriptional regulation highly relevant for psychiatry, and will contribute to our understanding of the combined effects of genetic and environmental factors in shaping disease risk Klengel & Binder, 2015.

N⁶-methyladenosine (m⁶A) is the most abundant internal mRNA modification, which is present transcriptome-wide in at least one fourth of all RNAs, typically located in a consensus motif (DRACH/GGACU), and enriched near stop codons and in 5’UTRs Dominissini et al., 2012Linder et al., 2015Meyer et al., 2012. Recent studies have identified mammalian m⁶A to be dynamically regulated, controlling stem cell proliferation and differentiation Klungland et al., 2017, cellular heat-shock response Zhou et al., 2015, DNA damage response Xiang et al., 2017 and tumorigenesis Cui et al., 2017. Brain RNA methylation is comparably high Meyer et al., 2012 and increases during development Meyer et al., 2012.

m⁶A is deposited co-transcriptionally Ke et al., 2017Slobodin et al., 2017 by a methyltransferase complex consisting of METTL3, METTL14 Liu et al., 2014, WTAP Ping et al., 2014, KIAA1429 (VIR) Schwartz et al., 2014, and RBM15/RBM15B (Patil et al., 2016). In contrast, it can be removed by the demethylases FTO Jia et al., 2011Mauer et al., 2017 and ALKBH5 Zheng et al., 2013. FTO further catalyses demethylation of N⁶,2′-O-dimethyladenosine (m⁶Am) with an in vitro preference for this substrate Mauer et al., 2017. m⁶Am is found at the first nucleotide of certain RNAs adjacent to the 7-methylguanosine cap and promotes transcript stability Mauer et al., 2017. The most commonly used anti-m⁶A antibody co-detects m⁶A and m⁶Am and therefore both m⁶A and m⁶Am are included in most of our measurements, henceforth both termed m⁶A without further specification. Fto has been associated with memory consolidation Walters et al., 2017Widagdo et al., 2016 and was implicated in regulation of dopaminergic brain networks Hess et al., 2013. m⁶A enzymes are expressed at different levels in different cell-types and have distinct intracellular distributions and binding motifs and thus potentially affect different subsets of target RNAs. Cellular consequences of m⁶A modifications depend on the binding of m⁶A-reader proteins (such as YTH- and HNRNP-proteins) and include RNA maturation, potentially splicing Alarcón et al., 2015Ke et al., 2017Liu et al., 2015Xiao et al., 2016, alternative polyadenylation Ke et al., 2015, RNA decay Shi et al., 2017Wang et al., 2014, as well as both promotion and inhibition of protein translation Meyer et al., 2015Shi et al., 2017Wang et al., 2015Zhou et al., 2015.

In this study, we aimed to elucidate the role of m⁶A in the context of the brain’s stress response. We delineated the effects of acute stress on m⁶A using global m⁶A measurements, m⁶A-Seq and absolute quantification of transcript-specific methylation levels. In addition, we explored the functional significance of m⁶A in the adult brain by examining conditional knockout mice for Mettl3 and Fto. Finally, we investigated m⁶A regulation in blood samples of mice and humans to determine its potential as a peripheral indicator of the central response to stress and stress-linked psychiatric disorders.

2Results¶

2.1The stress-induced m6A epitranscriptome¶

To test whether acute stress alters m⁶A, we performed m⁶A-Seq (RNA-Seq after m⁶A-immunoprecipitation) on mouse cortex 4 h following 15 min of acute restraint stress exposure. To facilitate the detection of biologically significant changes, we applied a minimum cut-off on expression level of RNA RPKM>1 (Figure S1 A). We observed a total of 25,821 m⁶A-peaks consistent across 7 biological replicates per condition, mapping to 11,534 genes (Figure 1 A, Table S1). On average, each transcript had 2.24 peaks (Figure S1 B). Of these peaks, 55.4 % overlapped with previously reported m⁶A-peaks (RMBase Sun et al., 2016) but we found an additional 3,199 novel m⁶A-methylated genes. Although 1,810 m⁶A-peaks (mapped to 1644 genes) were differentially methylated after stress (Table S1), only 179 transcripts were found to be expressed in a stress-affected manner (FDR-corrected P-Value < 0.05, absolute log2 fold change > 0.1) (Figure 1 A). Interestingly, only 16% (28) of these genes were regulated both on the RNA and m⁶A-level. Both the absolute extents of m⁶A fold changes (Figure S1 A) as well as RNA fold changes were rather small, which may reflect the cellular heterogeneity of the input material.

While the majority of peaks were located in the CDS (coding sequence) of transcripts (Figure 1 B), analysing the peak distribution along their length in respect to the landmarks of RNA transcripts Cui et al., 2016 revealed an enrichment both at the 5’UTR (5’ untranslated region) and around the stop-codon (Figure 1 C). The m⁶A consensus motif DRACH and the more conservative motif GGAC were found in 93.36 % and 68.68 % of the m⁶A-peaks, respectively, often with several occurrences per peak sequence, indicating potential multiple adjacent m⁶A-sites. Sequence motif analysis confirmed enrichment of the m⁶A consensus motif, with the top motif being GGACWB (Figure 1 D). Functional classification of m⁶A-peaks harbouring transcripts revealed over-representation of genes related to neuronal plasticity, morphogenesis and development with different biological processes and pathways enriched depending on the peak’s position (Figure S1 C).

Differential analysis reported peaks that were either unique to the control- or the stress-condition as well as m⁶A-peaks with quantitative changes between the conditions (Figure 1 E). Interestingly, both control-unique and higher-in-control peaks exhibited position preference for the 5’UTR, whereas stress-unique and higher-in-stress peaks were enriched in the 3’UTR (Figure 1 F and Figure S1 E). Differential methylation of chosen candidate transcripts was confirmed by m⁶A-RNA immunoprecipitation (RIP) followed by qPCR (Figure 1 G and Figure S1 D). Notably, differential methylation of some genes that were removed from analysis because of their low expression could also be confirmed or at least showed a similar trend of regulation (Figure S1 C).

Although only a fraction of transcripts was significantly regulated both on the m⁶A and RNA-level after correction for multiple testing, changes in m⁶A across all genes negatively correlated with the respective change in transcript levels (Figure 1 H, Table S2). This may indicate that higher methylation levels may functionally lead to reduced RNA levels. The same negative relationship was detected for transcripts with stress-altered m⁶A (Table S2). Interestingly, the degree of the inverse relationship depended on the peak’s position, with CDS peaks showing the strongest negative correlations, followed by 5‘UTR peaks (Figure S1 F, Table S2).

Since the function of altered m⁶A likely involves the action of reader-proteins, we next asked if methylation-altered transcripts are preferentially recruited to m⁶A-binding proteins. When we intersected m⁶A-regulated genes with those bound by m⁶A readers, we observed a higher than expected overlap for several but not all m⁶A readers (Figure 1 I, Q (FDR-corrected P-Value) < 0.05 for YTHDF1, YTHDF2, YTHDC1, HNRNPC), indicating that stress-regulated m⁶A-peaks are poised for downstream functional targeting by m⁶A binding proteins. Additionally, analysing the co-occurrence of the observed m⁶A-motif GGACWB to known binding motifs of RNA-binding proteins within the m⁶A-peak sequences revealed a higher than expected overlap and summit enrichment with the binding motifs of FMR1, a protein crucial for directed translocation of RNAs within neurons and critical for synaptic plasticity (Figure S1 G). Likewise, genes reported to be bound by the mouse FMR1 or the human homologue FMRP were also higher-than-likely m⁶A-methylated (mouse FMR1 Z = 9.25, Q < 0.05; human FMRP Z = 62.9, Q < 0.05; Figure S1 H), suggesting that m⁶A methylation of neuronal RNAs may regulate protein binding critical for neuronal transport and plasticity.

2.2Stress-regulation of m6A is brain-region-specific¶

Based on the fact that the fold-changes of stress-regulated m⁶A-peaks were rather small in the mouse cortex m⁶A-Seq, we reasoned that the true extent of the m⁶A-stress response may only be revealed when investigating more defined brain regions. Therefore, we measured the time course of RNA methylation changes in 2 regions highly involved in stress response regulation: the medial prefrontal cortex (PFC) and the basolateral and central amygdala (AMY; Figure 2 A). We found that global m⁶A was regulated in total RNA in a region-dependent manner with RNA methylation increased in the PFC and decreased in the AMY (Figure 2 B). The same regulation was observed when m⁶A was measured in mRNA using LC-MS/MS (Figure 2 C). To examine potential changes of the m⁶A-machinery related to these global changes, we measured gene expression levels of m⁶A enzymes and binding proteins. We found the demethylases Fto and Alkbh5 to be differentially regulated in a region-specific manner facilitating the effects observed on global methylation, in most cases preceding the effect observed on global m⁶A (Figure 2 D). Furthermore, Mettl3 was downregulated upon stress-exposure tissue-independently (Figure 2 D) and Wtap was regulated isoform-specifically only in the AMY (Figure S2 B). The m⁶A-reader Ythdc1 was regulated in a region-specific manner (Figure 2 D), whereas the other known enzymes and readers were not differentially expressed (Figure S2 A).

Notably, i.p. injection of corticosterone but not dexamethasone changed global m⁶A (Figure 2 E) as well as Fto and Alkbh5 expression (Figure S2 B) similarly to acute stress (Figure 2 D), demonstrating that the stress effect may be mediated by glucocorticoids (GCs). Supporting this idea, we found that the majority of m⁶A enzyme- and reader-genes contains several GC response elements in their 5’ upstream region, likewise pointing at expression regulation of those genes via GCs (Figure S2 C).

2.3Stress-regulation of m6A is gene-specific¶

m⁶A-Seq requires large amounts of input material thus limiting the capacity of the technique to investigate narrowly defined brain regions and also only measures relative enrichment changes and not absolute transcript methylation. Therefore, we performed m⁶A-immunoprecipitation of full length transcripts (m⁶A-RIP) followed by qPCR to assess absolute levels of candidate transcript methylation in narrowly defined brain areas, before and after stressful challenge. For calibration of the assay and normalization of immunoprecipitation-efficiency in downstream experiments, we designed and used an m⁶A-methylated internal spike-in RNA oligonucleotide (Figure 3 A and Figure S3 A, B). The m⁶A-RIP-qPCR detected m⁶A methylated RNA spike-in across a wide range of concentrations with low IgG-background signal and without competing with the immunoprecipitation of endogenously methylated RNAs (Figure 3 B). Using mixtures of unmethylated and methylated spike-in oligonucleotides, we confirmed that m⁶A-RIP-qPCR measured different methylation states of RNAs with high precision (Figure 3 C, r²>0.95).

Applying m⁶A-RIP-qPCR, we measured absolute methylation levels of several candidate transcripts involved in the brain’s stress response and, given the enrichment of neuronal plasticity and morphogenesis-related terms in the m⁶A-Seq, synaptic plasticity-related transcripts (Figure 3 D and Figure S3 C). Similar to the results of the m⁶A-Seq, regulation of m⁶A by stress (26/44 transcripts) was observed more often than regulation of RNA (16/44 transcripts, with 12 overlapping) in the transcripts tested. Notably, the majority of chosen candidates were either regulated or expressed in a region-specific manner, emphasizing the importance of assessing RNA-methylation in defined brain areas (Figure 3 E). The negative correlation of m⁶A- and RNA fold-changes was replicated in the m⁶A-RIP-qPCR data, with no influence of region and time point (Figure 3 F, Table S2). In detail, both PFC and AMY exhibited differential response both at 1 h and 4 h with opposite directions, paralleling the regulation observed in global m⁶A in the respective regions before (Figure 2 A). Overall, 4 h fold changes had higher effect sizes compared to 1 h fold changes (Figure 3 G). Fold changes at the 1 h time point correlated with those at 4 h for the same gene in the PFC but not AMY indicating that in the PFC 1 h m⁶A may be an intermediate state of 4 h regulation with fold-changes of regulated m⁶A increasing with time. In contrast, in the AMY for the candidate genes investigated, m⁶A regulation after 1 h and 4 h was more independent (Figure 3 H).

2.4Stress-coping behaviour is altered in mice deficient of Mettl3 or Fto¶

To explore the causal consequences of disturbed m⁶A levels on the stress response and stress-linked behaviours, we investigated mouse models of altered m⁶A-regulation in adult neurons. Since the expression of the m⁶A methyltransferase Mettl3 and the m⁶A and m⁶Am demethylase Fto were affected by acute stress, we generated conditional inducible knockout mouse models lacking these genes specifically in forebrain excitatory neurons by breeding Mettl3 or Fto flox/flox mice to tamoxifen-inducible Nex-CreERT2 mice (Mettl3-cKO and Fto-cKO). Upon induction with tamoxifen in young adults to prevent developmental effects, Mettl3 and Fto were depleted both from dorsal and ventral parts of the hippocampus, specifically in CA1 and CA3 but not in the dentate gyrus (Figure 4 A). Nex-CreERT2-induced recombination is further known to occur in small populations of principal neurons in the cortex (33). Depletion of either gene did not result in compensatory changes of gene expression of other genes involved in m6A-metabolism (Suppl. Figure S4 A), but altered transcriptome profiles as observed by mRNA-Seq of CA1 and CA3 tissue (Fig 4 B). Interestingly, in non-stressed basal animals, we observed a larger number of differentially expressed genes in Mettl3-cKOs compared to Fto-cKOs (Figure 4 B, C; Mettl3-cKOs: 205 differentially expressed genes with 164 genes with a fold change above log2=0.2; Fto-cKOs: 130 differentially expressed genes with 13 genes with a fold change above log2=0.2; Suppl. Table 3) as well as a wider spread of gene expression regulation (Supp. Figure S4 B) with no apparent preference for up- or downregulation. Although there was only small overlap of differential expressed genes between the two lines, 91 genes were differentially expressed in a knockout-specific pattern (Figure 4 C; 83 genes with a significant interaction effect of gene knockout and genotype and a fold change above log2=0.2; Suppl. Table 3) including crucial genes for regulation of neuronal activity response, and synaptic function as immediate-early genes as well as genes involved in synaptic functions (Figure 4 C). Neither Mettl3-cKO nor Fto-cKO mice showed altered anxiety-like behaviour or locomotion (assessed by the Open Field, Elevated Plus Maze and Dark-Light-Box tests; Figure S5 A) but we observed significant changes in spontaneous digging behaviour (Marble Burying Test, Figure S5 A). Both mouse models exhibited increased cued fear memory long-term maintained during memory extinction (Figure 5 A), suggesting that both perturbations in the m⁶A-system lead to long-term increased fear memory. Additionally, Fto-cKO mice exhibited increased contextual fear memory (Figure 5 A). Importantly, we observed no differences in non-fear-related memory or short-term working memory when tested in the object recognition test and Y-Maze spontaneous alternation test (Figure 5 A). Next, we investigated the transcriptional response patterns 24 h after fear conditioning comparing fear conditioned animals (“FC”) to control animals that experienced the same handling but no foot shock (“Box”). For both Mettl3-cKOs and Fto-cKOs we observed a large number of genes differentially expressed after fear conditioning in a genotype-dependent manner, implying a widely altered transcriptional response pattern after stress in animals with disturbed m6A-system (Figure 5 B). Thereby, significant gene-expression regulation was more extended in fear conditioned animals compared to non-fear conditioned animals (supp. Figure S4 C; Suppl. Table 3) with low overlap between both mouse models (Figure 5 C). In contrast to basal animals, Fto-cKOs showed more genotype-dependent expression changes after the stressful fear conditioning event than Mettl3-cKOs (Figure 5 B; Suppl. Table 3), implying that despite the major lack of Fto-cKO genotype effect in steady-animals, Fto is crucial for the regulation of the stress and fear response. Upon fear-conditioning stress, genotype-dependent transcriptomic changes involve genes crucial for neuronal systems like neurotransmitter receptors and transporters as well as transcription factors (Figure 5 C), pointing at a role of m⁶A in regulating synaptic function in the hippocampus after fear conditioning. Consequently, investigating the effects of Fto and Mettl3 depletion on electrophysiological correlates of network plasticity and brain function, we found that CA1 long-term potentiation was impaired in Fto-cKO but not in Mettl3-cKO mice (Figure 5 D) with no effect on paired-pulse facilitation (Figure 5 D) nor basal neurotransmission (Figure S5 B).

2.5Regulation of m6A is impaired in MDD patient blood¶

To evaluate the potential of blood m⁶A as a peripheral proxy of the central m⁶A stress response, we measured global m⁶A methylation levels in mouse and human blood after an acute stressful challenge and GC stimulation. Global methylation was transiently decreased in whole blood of mice after acute stress (Figure 6 A), with gene expression of Mettl3 and Alkbh5 altered in accordance with the global m⁶A change and Wtap being upregulated (Figure 6 B). Similarly, global m⁶A was decreased in mouse blood 4 h after i.p. injections of both corticosterone and dexamethasone (Figure 6 C). Comparably, blood from healthy human volunteers, drawn before and after intake of 1.5 mg dexamethasone, showed both reduced global m⁶A levels (Figure 6 D) and changes in gene expression of the m⁶A-machinery enzymes 3 h after dexamethasone intake (Figure 6 E) Arloth et al., 2015.

Since dysregulation of the stress response may be an important feature of psychopathologies like Major Depressive Disorder (MDD), we next investigated whether m⁶A-regulation in response to dexamethasone differs between healthy individuals and MDD patients. In contrast to healthy subjects, down-regulation of m⁶A in response to dexamethasone was neither observed in male nor female MDD patients (Figure 7 A). To prevent potential contamination of results by antidepressant treatment present in blood of MDD patients, we confirmed the lack of response to GC-stimulation using dexamethasone (Figure S5 A) and cortisol (Figure 7 B total RNA, Figure 7 C mRNA by LC/MS-MS) in B-lymphocyte cell lines (BLCLs) donated by each 6 healthy volunteers and 6 MDD-patients propagated in absence of antidepressants. NR3C1 (GC Receptor) mRNA and protein expression as well as transcriptional response to GC stimulation was unchanged in BLCLs of MDD donors (Figure S6 B-d). To estimate the transcriptome-wide distribution of m⁶A after GC stimulation, we performed m⁶A-Seq on BLCLs of one healthy donor treated for 1 h with 100 nM cortisol or mock conditions. We found over 17,000 m⁶A-peaks in around 9,000 genes in the healthy donor BLCLs (89% of peaks with GGAC motif, 99% of peaks with DRACH motif), with 12 % of those being cortisol-treatment-responsive (Figure 7 D, regulated peaks with absolute log2 fold change > 1, FDR-corrected P-Value < 0.05). Similar to the mouse m⁶A-Seq, we found the m⁶A consensus motif was enriched in the m⁶A-peaks (top motif GKACW) and m⁶A-peaks were distributed along the transcript-length with preference for both m⁶A-peaks in the 5’UTR and around the stop codon. Interestingly, in m⁶A-peaks that were found to be cortisol-responsive, we observed less preference for the 5’UTR with higher contribution of CDS and 3’UTR peaks (Figure 7 E,F). In parallel to the global demethylation observed in BLCLs of healthy donors after cortisol-stimulation, differential m⁶A-peaks were found to be almost exclusively hypomethylated with higher fold changes than observed in mouse brain (Figure 7g). To precisely quantify regulation of m⁶A-levels in BLCLs from healthy and MDD-donors, we performed m⁶A-RIP-qPCR testing for GC-responsive genes in BLCLs after stimulation with cortisol. We observed specific downregulation of m⁶A in FKBP5, IRS2 and TSC22D3 in cells from healthy but not from MDD individuals (Figure 7h). In line with the general trends observed before, methylation of tested candidates in cells derived from healthy, but not MDD donors, was significantly decreased (Figure 7i).

2.6Discusion¶

Here, we have identified m⁶A and m⁶Am as a widely regulated epitranscriptomic mark responsive to acute stress. Using m⁶A-Seq in mouse cortex, we observed over 25,000 m⁶A-peaks with around 7% being stress-responsive, thus far outnumbering the changes observed on the RNA expression level. Peaks were enriched at both the 5’UTR and around the stop-codon with a higher peak contribution in the 5’UTR than previously reported Meyer et al., 2015. In contrast to previous reports indicating changes in peaks primarily at the 5’UTR in cellular heat shock response Zhou et al., 2015, we observed no preference in location for differential peaks. However, we detected a preference for stress-downregulated peaks for the 5’UTR and for stress-upregulated peaks for the 3’UTR. Given the wide-variety of mechanisms of how RNA methylation may regulate RNA levels and translation, the molecular consequence of differential m⁶A regulation after stress is not yet clear. Stress-m⁶A regulated genes were significantly enriched for target genes of several m⁶A-readers as well as for the neuronal RNA-binding and cell-transport-regulating protein FMR1/FMRP recently shown to bind m⁶A (Edupuganti et al., 2017), suggesting that the effects of altered methylation may be equally broad. Overlap of FMRP-binding sites with synaptic transcripts bound by NSun2, a methyltransferase for another RNA modification, has been suggested before Hussain & Bashir, 2015. Although we observed a general negative correlation between m⁶A-change and RNA-abundance, most m⁶A-changes were not accompanied by significant transcript changes neither in m⁶A-RIP-qPCR nor in m⁶A-Seq. This may imply that differential m⁶A may act both by regulating RNA decay as well as location and translation control. Our m⁶A-Seq reported relatively small fold changes likely due to the large cellular heterogeneity of the material used with only a small fraction of cells being responsive to the treatment. Further, we propose that m⁶A regulation is highly specific to smaller brain areas with even potential opposite regulation in different areas leading to an underestimation of the actual regulation in m⁶A-Seq. This was supported by region-specific investigation of global and gene-specific m⁶A levels in the PFC and AMY. These 2 areas regulate behavioural and hormonal stress responses, fear and anxiety McEwen et al., 2015, with the PFC exhibiting top-down control of the AMY in anxiety and fear in mice Adhikari et al., 2015. These changes were accompanied by matching regulation in the demethylase Fto and Alkbh5 expression, as well as regulation of the methyltransferase Mettl3. Interestingly, previous reports also showed transcriptional regulation of Fto after acute stress by fear conditioning Walters et al., 2017Widagdo et al., 2016. Further, we observed regulation of the direct m⁶A-reader Ythdc1 after stress. In Drosophila, the Ythdc1 homologue YT521-B regulates neuronal function with behavioural defects in knockout-flies Lence et al., 2016. Notably, we observe that corticosterone i.p. injection in mice causes similar effects on m⁶A and enzyme-expression like acute stress, pointing towards a potential signalling mechanism via centrally acting GCs. Additional work is needed to unravel the pertinent signalling cascades involved. Single-cell RNA-Seq data showed that all m⁶A enzymes and readers are expressed in all major brain cell types including glia Zeisel et al., 2015. Therefore, it is unclear which of these cell types drive the observed expression changes. To investigate the mechanisms of m⁶A-methylation in neurons only, we specifically depleted Mettl3 and Fto from adult excitatory neurons of the hippocampus and cortex. We observed increased fear memory for both cued and contextual fear (latter in Fto-cKO mice only), which may be indicative of impaired stress-coping. Increased fear expression was previously reported upon knock-down of Fto in the dorsal hippocampus Walters et al., 2017 and in the PFC Widagdo et al., 2016. Extending the previous findings, we found fear memory to be enhanced upon disruption of both m⁶A methyltransferase Mettl3 or demethylase Fto, indicating that fear-memory requires fine-tuned regulation of m⁶A-levels rather than being directly regulated by m⁶A levels. Anxiety-like behaviour was not altered in either mouse line but may be regulated by m⁶A affected by different enzymes or in different brain areas and cell types other than those investigated here. Importantly, we did not observe other impairments in cognition or memory in these mice, pointing towards a potential importance of m⁶A for fear-specific memory. Mechanistically, we found that both Mettl3 and Fto depletion not only alter the steady-state transcriptome in adult hippocampal neurons but also the transcriptomic response to the fear conditioning stress, including regulation of several genes involved in neuronal circuits function pointing out a function of m⁶A in regulating neuronal circuits. Consequently, we describe that network plasticity is specifically altered in the CA1, a brain region crucial for contextual fear, in Fto-cKO, but not Mettl3-cKO mice. This may reflect a neuronal correlate of altered m⁶A underlying the altered contextual fear memory observed in Fto-cKO mice. Several arguments point at both enzymes targeting largely non-shared methylation sites and/or genes in the specific case investigated, either by different m⁶A sites affected or by Mettl3 targeting predominantly m⁶A and Fto predominately targeting m⁶Am: First, we did not observe oppositely directed behavioural or electrophysicological effects in both lines. Second, the gene sets affected in transcription by Mettl3 and Fto depletion were widely non-overlapping. Lastly, both enzymes showed different contribution to steady-state transcriptomic regulation (mainly Mettl3-cKOs) and fear conditioning-regulated gene expression (dominantly Fto-cKOs).

Finally, we propose that regulation of m⁶A and its cellular machinery in blood may represent a peripheral proxy for part of the brain’s m⁶A response, similar to DNA methylation changes Ewald et al., 2014Provençal et al., 2012. Both mice and humans showed global blood demethylation after stress or GC intake, respectively. The m⁶A-response to GCs is impaired in blood and blood cells obtained from MDD patients, which may be a consequence of the altered GC receptor reactivity and downstream signalling reported in MDD Kloet et al., 2005. Interestingly, genetic variants in FTO Milaneschi et al., 2014Samaan et al., 2013 and ALKBH5 Du et al., 2015 have been reported to associate with risk for MDD before but are yet to be replicated in larger cohorts. There is growing evidence that fine-tuned transcriptional regulation is especially relevant for psychiatric disorders including disease-associated SNPs in enhancer regions (Cross-Disorder Group of the Psychiatric Genomics Consortium, 2013) and epigenetic changes Klengel & Binder, 2015 including regulation of chromatin conformation Won et al., 2016, histone modifications The Network and Pathway Analysis Subgroup of the Psychiatric Genomics Consortium, 2015 as well as short and long ncRNAs Issler & Chen, 2015Parikshak et al., 2016. Here, we reveal RNA modifications as a novel layer of regulation of gene expression that is potentially crucial for these disorders. Integrating all layers of transcriptional regulation will be highly relevant for understanding psychiatric disorders.

In summary, m⁶A and m⁶Am methylation constitute a novel layer of complexity in gene expression regulation following stress exposure, which is pivotal for the adaptation of stress-responsive circuits to acute challenges. The exciting finding of m⁶A dysregulation in MDD opens the possibility for the development of novel diagnostic biomarkers and eventually to better treatments for anxiety disorders, depression and other stress-related diseases.

3Author contributions¶

M. En., and A.C. conceived and designed the experiments and wrote the manuscript. M.En. performed, and analysed most experiments. M. En. and S.R. performed bioinformatics analysis. C. E., P. M. K, L.T., and M.R.-H. assisted in experiments. M. J. performed preliminary experiments with no data used. M. U. performed and analysed mass spectrometry experiments. M. Ed. performed and analysed electrophysiological studies. J.A. and E. B. B. provided human microarray data. P.W., C. T. W., M. V. S., J. M. D., E. B. B., and A.C. conceived and designed the project. The project was supervised by A.C.

5Tables¶

None.

6STAR Methods¶

6.1Contact for reagent and resource sharing¶

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Alon Chen, (alon_chen@psych.mpg.de).

6.2Experimental model and subject details¶

6.2.1Animals¶

All experiments were approved by and conducted in accordance with the regulations of the local Animal Care and Use Committee (Government of Upper Bavaria, Munich, Germany and Weizmann Institute of Science, Rehovot, Israel).

For all experiments characterizing m⁶A changes after stress, 10-12 w old adult C57 BL/6 male mice were used (Charles River, Sulzfeld, Germany). Mettl3-cKO and Fto-cKO mice were generated by breeding Mettl3^{tm1 C(KOMP)Wtsi} lox/lox mice Geula et al., 2015 and Fto^{tm1 C(EUCOMM)Wtsi} lox/lox mice obtained from EMMA (EM:05094) to Nex-CreERT2 mice Agarwal et al., 2012, respectively. Experimental mice were homozygous floxed Nex-CreERT2-positive and Nex-CreERT2-negative littermates generated by breeding of homozygous floxed mice negative and hemizygous for the CreERT2-allele. All experimental animals of the 2 knockout lines were fed with tamoxifen-containing chow (Genobios LASCR diet Cre Active TAM 400) starting at the age of 4-6 w. Animals were housed in groups until being single housed 7 d before the experiments started in standard plastic cages and maintained in a temperature-controlled environment (21 ± 2°C) on a 12 h light/dark cycle with food and water available ad libitum. Restraint stress was performed for 15 min in ventilated 50 ml falcon tubes, starting at 2 h post lights on. For pharmacological studies, mice were injected with vehicle solution (saline), 250 μg/kg corticosterone (corticosterone-HBC complex, Sigma) or 10 mg/kg dexamethasone (Ratiopharm Dexa-ratiopharm) i.p. 2 h post switching the lights on.

6.2.2Sample collection¶

Whole mouse cortex for m⁶A-Seq was collected at designated time points by manual dissection of fresh brains on ice. For each sample, 3 animals randomly selected from the same group were pooled. For investigation of regions-specific effects in PFC and AMY, brains were immediately flash-frozen after dissection and defined tissue punches of medial prefrontal cortex (PFC; consisting of infralimbic and prelimbic cortex) and amygdala (AMY; consisting of central and basolateral amygdala) were collected using a 1 mm round tissue punch while sectioning brains on a cryostat. Mouse whole blood was collected in EDTA tubes, aliquoted and flash-frozen.

6.2.3Cell culture¶

Human immortalized BLCLs derived from age-matched (33-53 y) male subjects either healthy or diagnosed with MDD were cultured in RPMI-1640 medium (Merck KGaA, Darmstadt, Germany) supplemented with 10 % fetal calf serum at 37 ^°C with 5 % CO2. The cells were tested to be free of mycoplasma. Cells were treated with cortisol (Sigma-Aldrich, St. Louis, MO, in ethanol, final concentration 0.1% v/v) or dexamethasone (Ratiopharm Dexa-ratiopharm, in saline), or ethanol or saline mock control, respectively.

6.2.4Human blood¶

Human whole blood was collected using PAXgene Blood RNA Tubes (PreAnalytiX, Hombrechtikon, Switzerland) either unstimulated or after oral administration of 1.5 mg dexamethasone and processed as described previously Menke et al., 2012. Age-matched healthy Caucasian male and females subjects were selected from the “MPIP” and “MARS” cohorts described previously Arloth et al., 2015Menke et al., 2012.

6.3Method Details¶

6.3.1RNA isolation¶

Total RNA from tissue, mouse blood and BLCL cells was purified using Trizol (Invitrogen, Life Technologies, Carlsbad, CA) according to the manufacturer’s instructions followed by isopropanol precipitation. For mouse whole blood, RNA was isolated using a 1:10 ratio of blood to Trizol.

6.3.2Global m⁶A measurements¶

Global m⁶A in total RNA was quantified by the EpiQuik m⁶A RNA Methylation Quantification Kit (Epigentek Group Inc., Farmingdale, NY) following manufacturers' specifications and using 100-300 ng input (in duplicates or triplicates). Comparing total RNA global m⁶A measurements with LC-MS/MS data from the same conditions, we observed high correlation of stress-changes, suggesting that the total RNA colorimetric assay represents an appropriate tool to detect global m⁶A regulation patterns. Brain global methylation in PFC and AMY is not regulated by circadian rhythm (data not shown).

6.3.3LC-MS/MS¶

Samples were pooled from 4 mice randomly selected from the same group. Residual genomic DNA was removed using the TurboDNA-free kit (Ambion, Life Technologies, Carlsbad, CA). RNA integrity and absence of DNA was confirmed by Bioanalyzer RNA Nano chips (Agilent Technologies, Santa Clara, CA, RIN >9) and Qubit DNA High sensitivity kit (Thermo Fisher Scientific, Waltham, MA), respectively. PolyA+ RNA was prepared using 2 rounds of the Genelute mRNA Prep Kit (Sigma-Aldrich, St. Louis, MO) with rRNA depletion confirmed by Bioanalyzer RNA Nano chips (Agilent Technologies, Santa Clara, CA, mRNA mode). 250 ng PolyA-RNA per sample and a N6-methyladenosine/adenosine standard curve were mixed with deuterated N6-(methyl-d3)-adenosine as an internal spike-in calibrator and processed and measured as reported before Jia et al., 2011. Quantification was performed by comparison with the standard curve obtained from pure nucleoside standards normalized by the deuterated spike-in calibrator run within the same experiment.

6.3.4m⁶A-Seq¶

For mouse m⁶A-Seq, whole mouse cortex samples were used pooling 3 individuals each, since m⁶A-Seq on PolyA-RNA of smaller regions did not result in sufficient enrichment quality. For m⁶A-Seq of human BLCLs, RNA from one cell line (healthy male donor) 1 h after treatment with 100 nM cortisol or mock condition was used with 1 technical replicate per condition and no IgG control. m⁶A-Seq was preformed using the previously published m⁶A-Seq-protocol Dominissini et al., 2013 (BLCL samples) or the following slightly modified version of it (mouse brain): Residual genomic DNA was removed using the TurboDNA-free kit (Ambion, Life Technologies, Carlsbad, CA). RNA integrity and absence of DNA was confirmed by Bioanalyzer RNA Nano chips (Agilent Technologies, St. Louis, MO, RIN > 9.5) and Qubit DNA High sensitivity kit, respectively. PolyA+ RNA was prepared using 1 round of the Genelute mRNA Prep Kit (Sigma-Aldrich, St. Louis, MO) with less than 5 % residual rRNA as confirmed by Bioanalyzer RNA Nano chips (Agilent Technologies, St. Louis, MO, mRNA mode). RNA was fragmented using fragmentation reagent (Life Technologies, Carlsbad, CA). mRNA fragments were precipitated with ethanol and used for m⁶A-immunoprecipitation, IgG control and input samples. m⁶A-immunoprecipitation (10 ug mRNA fragments, 10 ug rabbit polyclonal anti-m⁶A 202 003, Synaptic Systems, Göttingen, Germany) or IgG control (10 ug mRNA fragments mixed from all samples, 10 μg IgG 2729, Cell Signalling Technology, Beverly, MA) was performed in precipitation buffer (50 mM Tris, pH 7.4, 100 mM NaCl, 0.05% NP-40, 1 ml total volume) with 1 μl RNasin Plus (Promega, Madison, WI) rotating head over tail at 4 °C for 2 h, followed by incubation with washed 30 μl Protein A/G beads (Thermo Fisher Scientific, Waltham, MA) rotating at 4 °C for 2 h. Bead-bound antibody-RNA complexes were recovered on a magnetic stand and washed twice with immunoprecipitation buffer, twice with high-salt buffer (50 mM Tris, pH 7.4, 1 M NaCl, 1 mM EDTA, 1 % NP-40, 0.1 % SDS), and twice with immunoprecipitation buffer. Fragments were eluted by Proteinase K treatment (300 μl elution buffer: 5 mM Tris-HCL pH 7.5, 1 mM EDTA pH 8.0, 0.05 % SDS, 4.2 μl 20 mg/ml proteinase K). RNA was recovered from the eluate using Trizol LS (Invitrogen Life Technologies, Carlsbad, CA) following manufacturers' recommendations. Sequencing libraries were prepared using the Illumina TruSeq non-stranded mRNA protocol following the standard protocol starting from mRNA fragments (mouse brain) or the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, Ipswich, MA). Libraries were quality-checked using Bioanalyzer DNA High Sensitivity chips (Agilent Technologies, St. Louis, MO) and quantified using the KAPA Library Quantification Kit (KAPA Biosystems, Boston, MA). Sequencing was performed on 4 lanes of an Illumina HiSeq4000 PE 2x100 (mouse brain) or on 1 lane of an Illumina MiSeq PE 2x 75 (BLCLs, Illumina, San Diego, CA) multiplexing all m⁶A-, IgG- and input samples.

6.3.5mRNA-Seq¶

Brains were collected from 5 of each of the following: Mettl3-cKO and WT as well as Fto-cKO and WT mice 24 h after fear conditioning (“FC”, details in “Animal behaviour testing”) or comparable handling without fear induction (“Box”: handling and exposure to context as in “FC” in “Animal behaviour testing” but without foot shock and tone/CS and US). The entire CA1 and CA3 was cryo-punched using 0.7 and 1 mm punching tools from snap-frozen brains sliced at 250 μm using a cryostat and RNA isolated. Residual genomic DNA was removed using the TurboDNA-free kit (Ambion, Life Technologies, Carlsbad, CA). RNA integrity and absence of DNA was confirmed by Bioanalyzer RNA Nano chips (Agilent Technologies, St. Louis, MO, RIN > 8.5) and Qubit DNA High sensitivity kit, respectively. mRNA-Seq libraries were prepared from 4 μg total RNA using the llumina TruSeq stranded mRNA protocol HT (Illumina, San Diego, CA) following the standard protocol starting using Superscript III and 11 cycles of PCR. Libraries were quality-checked using Bioanalyzer DNA High Sensitivity chips (Agilent Technologies, St. Louis, MO) and quantified using the KAPA Library Quantification Kit (KAPA Biosystems, Boston, MA). Sequencing was performed on 4 lanes of an Illumina HiSeq4000 PE 2x100 (Illumina, San Diego, CA) multiplexing all samples.

6.3.6Gene expression¶

Gene expression of m⁶A-related enzymes was done by SYBR-green-based qPCR. RNA was reverse-transcribed using the SuperScript III VILO cDNA Synthesis Kit (Invitrogen Life Technologies, Carlsbad, CA) and QuantiFast SYBR Green PCR Kit (QIAGEN, Hilden, Germany) on a Quantstudio 7 (Applied Biosystems, Waltham, MA) with the following primers: Mettl3NM_019721(ATTGAGAGACTGTCCCCTGG, AGCTTTGTAAGGAAGTGCGT), Mettl14 NM_201638 (AGACGCCTTCATCTCTTTGG, AGCCTCTCGATTTCCTCTGT), Wtap_consensus (GTTATGGCACGGGATGAGTT, ATCTCCTGCTCTTTGGTTGC),Wtap_short NM_001113532 (CTAGCAACCAAAGAGCAGGA, AGTCTTGACTGGGGAGTATGA), Wtap_long NM_001113533 (GGCAAAAAGCTAATGGCGAA, GCTGTCGTGTCTCCTTCAAT), Fto NM_011936 (CTGAGGAAGGAGTGGCATG, TCTCCACCTAAGACTTGTGC), Vir-Kiaa1429NM_001081183(CATTACGGCCGCTTAGTTCT, TACCACTGCCTCCACTAACA), Alkbh5 NM_172943 (ACAAGATTAGATGCACCGCG, TGTCCATTTCCAGGATCCGG), Ythdf1 NM_173761 (CATTATGAGAAGCGCCAGGA, AGATGCAACAATCAACCCCG), Ythdf2 NM_145393 (ACCAACTCTAGGGACACTCA, GGATAAGGAGATGCAACCGT), Ythdf3 NM_172677 (TGCACATTATGAAAAGCGTCA, AGATGCGCTGATGAAAACCA), Ythdc1 NM_177680 (TTCATAACATGGGACCACCG, TCATAGTCATGTACTCGTTTATCTC),HnrnpcNM_016884 (CAAACGTCAGCGTGTTTCAG, TGGGGATGAGAAGGACAAGT), Hnrnpa2 B1 NM_016806 (GTGGAGGGAACTATGGTCCT, TGAAGGCACCAACAAGAACT). Each qPCR assay was performed in duplicates or triplicates with a standard dilution curve of a calibrator and using assay efficiency for calculations. Expression levels were quantified by the ddCT method normalizing to an average of 4-5 housekeeping genes chosen based on maximum stability between conditions from the following: Hprt NM_013556 (ACCTCTCGAAGTGTTGGATACAGG, CTTGCGCTCATCTTAGGCTTTG), Rpl13 A NM_009438 (CACTCTGGAGGAGAAACGGAAGG, GCAGGCATGAGGCAAACAGTC), Atp5j NM_001302213(TATTGGCCCAGAGTATCAGCA, GGGGTTTGTCGATGACTTCAAAT), Polr2B NM_153798 (CAAGACAAGGATCATATCTGATGG, AGAGTTTAGACGACGCAGGTG), Rn18s NR_003278 (CAGGATTGACAGATTGATAGC, ATCACAGACCTGTTATTGCTC), Ubc NM_019639 (CTGCCCTCCCACACAAAG, GATGGTCTTACCAGTTAAGGTT), Hmbs NM_001110251(TCTGAAAGACAGATGGAATGCC, CCACACGGAAAGAGAAGAGGC). For human samples, the following primers were used: NR3C1 NM_000176 (CAGCAGTGAAATGGGCAAAG, TCGTACATGCAGGGTAGAGT), NR3C2 NM_000901 (GATCCAAGTCGTGAAGTGGG, TGAAGGCTGATTTGGTGCAT), FKBP5 NM_004117 (CGGCGACAGGTTCTCTACTT, TCTCCAATCATCGGCGTTTC), TSC22D3 NM_004089 (TCCGTTAAGCTGGACAACAG, TTCAACAGGGTGTTCTCACG) with housekeeping genes TBP NM_003194 (GGGAGCTGTGATGTGAAGTT, GAGCCATTACGTCGTCTTCC), RPL13 A NM_012423 (GCGTCTGAAGCCTACAAGAA, CCTGTTTCCGTAGCCTCATG), and SDHA NM_004168 (CAGGGAAGACTACAAGGTGC, CAGTCAGCCTCGTTCAAAGT).

6.3.7Upstream GRE prediction¶

10 kb upstream sequences of m⁶A-related genes were retrieved using Biomart Smedley et al., 2015. GC response elements were predicted by the JASPAR vertebrate core transcription factor binding site prediction Mathelier et al., 2016 querying NR3C1 motifs MA0113.1 (mammalian), MA0113.2 (mmu), and MA0113.3 (hsa) with a conservative relative profile score threshold of 90 %.

6.3.8Spike-in Oligo¶

The spike-in RNA oligo was designed with the following specifications: 100 bp length, 3 internal m⁶A sites within GGAC motif flanked by the most frequent nucleotides 5’ U/A, 3’ A/U, not complementary to hsa or mmu RefSeq mRNA or genome, secondary structure exposing m⁶A sites, mean % GC = 51. The sequence is GCAGAACCUAGUAGCGUGUGGmACACGAACAGGUAUCAAUAUGCGGGUAUGG mACUAAAGCAACGUGCGAGAUUACGCUGAGGmACUACAAUCUCAGUUACCA.

Fully m⁶A-methylated or unmethylated RNA oligos were purchased from Sigma (Sigma-Aldrich, St. Louis, MO). m⁶A site prediction was performed using SRAMP Zhou et al., 2016 (full transcript mode, generic predictive model) confirming that the motif sequence context is similar to those occurring in real m⁶A data. Structure prediction was performed using RNAstructure Reuter & Mathews, 2010 (Fold mode, Version 5.8.1).

6.3.9Candidate m⁶A-RIP-qPCR¶

To validate m⁶A-Seq experiments, candidates were chosen from the list of differentially methylated transcripts selecting for transcripts with only 1 or few m⁶A-peaks, in the latter case with either only 1 of them stress-regulated or all regulated in a similar manner. For investigation of candidate transcript methylation in small brain areas, candidate lists were constructed by intersecting microarray results of mouse brain PFC, AMY and hippocampus after acute stress and GC stimulation Arloth et al., 2015 with genes known to be methylated in mouse brain Hess et al., 2013Meyer et al., 2012 and functional annotation GO-terms. For investigation of candidate transcript methylation in BLCL cell lines, dexamethasone-responsive genes from human blood microarray data Arloth et al., 2015 were intersected with BLCL m⁶A-Seq data (unpublished).

15 μg PolyA-RNA (m⁶A-Seq validation) or 3 μg total RNA (candidate m⁶A-RIP-qPCR validation) or 1.5 μg total RNA (brain area/cell line candidate m⁶A-RIP-qPCR) was mixed well with 30 fmol or indicated amount of spike-in or 3 fmol spike-in, respectively, and equally split into 3 conditions: m⁶A-RIP, IgG control and input. For m⁶A-Seq validation and brain are candidate m⁶A-RIP-qPCR only fully methylated spike-in was used. Input samples were flash-frozen during the course of the experiments. m⁶A-RIP and IgG control samples were incubated in parallel to m⁶A-Seq with 1 ug anti-m⁶A antibody (rabbit polyclonal 202 003, Synaptic Systems, Göttingen, Germany) or 1 ug IgG (rabbit polyclonal IgG 2729, Cell Signalling Technology, Beverly, MA) in immunoprecipitation buffer (0.5 ml total volume) with 1 μl RNasin Plus (Promega, Madison, WI) rotating head over tail at 4 °C for 2 h, followed by incubation with washed 25 μl Dynabeads M-280 (Sheep anti-Rabbit IgG Thermo Fisher Scientific, Waltham, MA, 11203 D) rotating head over tail at 4 °C for 2 h. Bead-bound antibody-RNA complexes were recovered on a magnetic stand and washed twice with immunoprecipitation buffer, twice with high-salt buffer, and twice with immunoprecipitation buffer. RNA was eluted directly into Trizol and input RNA was also taken up in Trizol. RNA from all conditions was purified in parallel using the miRNeasy micro RNA isolation kit (QIAGEN, Hilden, Germany) including a 3-time repeated elution 15 μl H2O to ensure the complete elution of all RNA. The entire eluate was transcribed to cDNA using the SuperScript III VILO cDNA Synthesis Kit (Invitrogen, Life Technologies, Carlsbad, CA). Gene expression was quantified using TaqMan Fast Advanced Master Mix (Applied Biosystems, Waltham, MA) on a Quantstudio 7 (Applied Biosystems, Waltham, MA) by the following Taqman gene expression assays: Actb NM_007393 (Mm01205647_g1), Akt1 NM_009652 (Mm01331626_m1), Arc NM_018790 (Mm01204954_g1), Atp1 B1 NM_009721 (Mm00437612_m1), Bsn NM_007567 (Mm00464452_m1), Camk2 A NM_009792 (Mm00437967_m1), Camk2n1 NM_025451 (Mm01718423_s1), Cited1 NM_007709 (Mm01235642_g1), Cnr1 NM_007726 (Mm01212171_s1), Crh NM_205769 (Mm04206019_m1), Crhbp NM_198408 (Mm01283832_m1), Crhr1 NM_007762 (Mm00432670_m1), Ctsb NM_007798 (Mm01310508_g1), Cyfip2 NM_133769 (Mm00460148_m1), Dlg4 NM_007864 (Mm00492193_m1), Dnmt1 NM_001199433 (Mm01151063_m1), Dusp1 NM_013642 (Mm00457274_g1), Egr3 NM_018781 (Mm00516979_m1), Fkbp5 NM_010220 (Mm00487406_m1), Fscn1 NM_007984 (Mm00456046_m1), Fth1 NR_073181 (Mm04336020_g1), Gabbr1 NM_019439 (Mm00444578_m1), Gabbr2 NM_001081141 (Mm01352561_m1), Gadd45g NM_011817 (Mm01352550_g1), Grm1 NM_001114333 (Mm00810219_m1), Grm3 NM_181850 (Mm01316764_m1), Homer1 NM_011982 (Mm00516275_m1), Htra1 NM_019564 (Mm00479887_m1), Mllt11 NM_019914 (Mm00480176_m1), Nlgn2 NM_198862 (Mm01245481_g1), Nodal NM_013611 (Mm00443040_m1), Notumos AK028718 (Mm00845023_s1), Nr3c1 NM_008173 (Mm00433832_m1), Nr4a1 NM_010444 (Mm01300401_m1), Nrcam NM_176930 (Mm00663607_m1), Nrxn1 NM_020252 (Mm03808856_m1), Nrxn2 NM_020253 (Mm01236844_g1), Onecut1 NM_008262 (Mm00839394_m1), P2ry13 NM_028808 (Mm00546978_m1), Plekhg3 NM_153804 (Mm00770086_m1), Plin4 NM_020568 (Mm00491061_m1), Pomc NM_008895 (Mm00435874_m1), Prkcb NM_008855 (Mm00435749_m1), Prkcg NM_011102 (Mm00440861_m1), Pvrl3 NM_021495 (Mm01342993_m1), Rgs4 NM_009062 (Mm00501392_g1), Rhou NM_133955 (Mm00505976_m1), Sgk1 NM_001161850 (Mm00441387_g1), Sgk2 NM_013731 (Mm00449845_m1), Sirt2 NM_022432 (Mm01149204_m1), Spats1 NM_027649 (Mm01270591_m1), Sumo1 NM_009460 (Mm01609844_g1), Syn1 NM_013680 (Mm00449772_m1), Syngap1 NM_001281491 (Mm01306145_m1), Tec NM_013689 (Mm00443230_m1), Tsc22d3 NM_001077364 (Mm00726417_s1). Mouse housekeeping genes: Hprt1 NM_013556 (Mm03024075_m1), Rpl13a NM_009438 (Mm01612987_g1), Tbp NM_013684 (Mm01277045_m1), Ubc NM_011664 (Mm02525934_g1), Uchl1 NM_011670 (Mm00495900_m1). Human gene expression assays: ID3 NM_002167 (Hs00171409_m1), DUSP1 NM_004417 (Hs00610256_g1), DDIT4 NM_019058 (Hs01111686_g1), GPER NM_001505 (Hs01922715_s1), IRS2 NM_003749 (Hs00275843_s1), FKBP5 NM_004117 (Hs01561006_m1), NR3C1 NM_000176 (Hs00353740_m1), TSC22D3 NM_004089 (Hs00608272_m1). Human housekeeping genes: RPL13A NM_012423 (Hs04194366_g1), TBP NM_003194 (Hs00427620_m1). The spike-in was quantified using a custom Taqman expression assay (primers TCAATATGCGGGTATGGACTAAAGC, TGAGGACTACAATCTCAGTTACCA and probe AACGTGCGAGATTACG).

6.3.10Human microarray data¶

Gene expression of m⁶A-related genes was extracted from microarray expression data of human whole blood published previously Arloth et al., 2015.

6.3.11Animal behaviour testing¶

All behavioural assessments were performed during the light phase. The experimenter was blinded to the genotype of the animals. Retesting followed the order of least-to-most stressful with 2-3 days’ rest in between tests.

Anxiety-like behaviour was assessed using the Open Field Test (OF, 10 min, 10 lux, grey plastic box 50 × 50 × 50 cm, centre defined as the inner 25 x 25 cm area), Elevated Plus Maze (EPM 5 min, 10 lux on closed arms, 100 lux on open arms, grey plastic maze 50 × 50 cm elevated 25 cm above the floor), Dark Light Box-Test at baseline (DLB basal) and 4 h post 15 min restraint stress (DLB 4 h post stress) (5 min, 100 lux in lit compartment), each with automated tracking (ANY-maze, Stoelting Co., Wood Dale, IL). The Marble Burying Test was performed by placing the mice in a fresh cage with 5 cm flattened fresh bedding with 15 black, clean marbles spaced evenly across (20 min, 10 lux, counting the number of buried marbles every 5 min). Cognitive function was assessed using the Y maze alternation task for working memory (5 min, 10 lux, Y-shaped 3-arm apparatus with 25 cm arm length and distinguishing visual cues on the walls and at the end of each arm, with automated tracking). The proportion of spontaneous non-repeated subsequent entries into each of the 3 arms (alternations) from the total number of 3-arm entries (including repeat entries) was used as the readout. Non-fear-related memory was assessed using the Object-Recognition-Task (ORT, 2x 5 min with 1 h intertrial interval, 10 lux, grey plastic box 50 × 50 × 50 cm, training trial: 2 identical objects with 1 out of 2 objects without object preference randomly assigned to all mice, test trial: 1 known, 1 novel object). The object discrimination ratio DI was determined by DI = (Time with novel object–Time with familiar object) / (Time with novel object + Time with familiar object) within the test trial.

Fear-related memory was assessed by conditional fear learning. Mice were fear conditioned (FC) within the same session for both contextual and cued fear by 180 s of baseline exposure to context A (a metallic/plastic cubic chamber with metal grid conditioned with 70% ethanol smell), followed by a 20 s 80 dB tone (9 kHz sine-wave, conditioned stimulus, CS), which co-terminated with an electric foot shock (unconditioned stimulus, US, 0.7mA, 2s, constant current delivered through the metal grid) and a 60 s after-shock interval. Memory was assessed by measuring freezing in response to the different cues by a highly experienced observer blind to the genotype. Auditory cued fear memory was tested 1 day after FC in context B (cylindrical plastic chamber with bedding conditioned with 1 % acetic acid) by presenting a 3 min CS after 180 s baseline recording and followed by a 60 s post-tone recording. Freezing across the 180 s tone exposure was binned in 60 s intervals to assess short-term stimulus-habituation. Context memory was tested 2 days after FC in context A without presenting US or CS. Fear extinction was achieved by 10*20 s CS presentations (variable inter-trial-interval of 20-60s) in context B on 3 consecutive days 2 w after FC with freezing assessed across the first 3 tone presentations. Fear extinction memory retention was measured 1 w after the extinction by presenting 3*20 s CS and measuring the freezing across those 3 presentations. Animals with generalized fear response (over 50% freezing in any of the baseline recordings of the extinctions trials) were excluded from the analysis for extinction memory.

6.3.12Electrophysiology¶

Recordings were conducted blind to the animal genotype. Preparation of dorsal hippocampal slices and electrophysiological measurements were performed according to standard procedures as we described previously Schmidt et al., 2011. From every animal, 2 slices were used for the experiments.

6.3.13In situ hybridization¶

Expression quantification of Mettl3 mRNA and Fto mRNA in Mettl3-cKO and Fto-cKO animals was performed by in situ hybridization using S-35 labelled antisense probes targeting the floxed exon as described previously Refojo et al., 2011. Probes were designed for Mettl3 NM_019721 exon 4 (probe cloned using TCAGTCAGGAGATCCTAGAGCTATT and CTGAAGTGCAGCTTGCGACA) and Fto NM_011936 exon 4 (probe cloned using TGGCAGCTGAAATACCCTAAACT and ATAGCTGTACACTGCCACGG). Slides were exposed to Kodak Biomax MR films (Eastman Kodak Co., Rochester, NY), developed, and autoradiographs digitized and quantified by optical densitometry of 2 slides each averaging the signal across both hemispheres and slides utilizing ImageJ (dorsal: Bregma −1.82, −1.94; ventral Bregma −3.16, −3.28).

6.3.14Western Blot¶

Cells were lysed on ice in RIPA buffer (150 mM NaCl, 1 % NP-40, 0,5 % Sodium deoxycholate, 0,1 % SDS, 50 mM Tris-HCl pH 8 with cOmplete, EDTA-free Protease Inhibitor Cocktail Mini, Roche Applied Science, Roche Diagnotics, Indianapolis, IN) for 30min. 25 μg total protein as determined by the Bio-Rad Quick Start Bradford Kit (Bio-Rad Laboratories Inc., Hercules CA) was heated for 10 min in SDS/PAGE sample buffer (final concentration 62.5 mM Tris-HCl pH 6.8, 2% SDS, 10% glycerol, 5% b-mercapto-ethanol), separated on a Tris-Glycine SDS–PAGE (Bio-Rad Laboratories Inc., Hercules, CA) and transferred to a nitrocellulose membrane (Amersham Protran, Millipore, Billerica, MA). Membranes were blocked for 1 h in TBST containing 5% non-fat milk, followed by incubation with primary antibodies overnight at 4 °C (anti-GR monoclonal rabbit, ab109022, 1:50000; anti-BTUB polyclonal rabbit, ab6046, 1:10000; Abcam, Cambridge, UK) in TBST with 3% non-fat milk. After incubation with horseradish-peroxidase-coupled secondary antibody (Cell Signalling Technology, Beverly, MA; 7074) at room temperature for 2 h, immunoblots were visualized using enhanced chemiluminescence (ECL Plus, GE Healthcare Life Sciences, Freiburg, Germany). Band intensity was quantified using ImageJ.

6.4Quantificationandstatisticalanalysis¶

6.4.1m⁶A-Seq analysis¶

Sequencing data quality control was performed by FASTQC Andrews, 2010. Genomic alignment was performed for the mouse brain samples using the STAR aligner Dobin et al., 2013 (to Gencode M11/Ensemble 86 build, mm10, at default settings using only those reads mapping uniquely to the genome and with a Phred quality score R 20) and for BLCL samples using TopHat2 Kim et al., 2013 (to UCSC GRCh37/hg19 at default settings). We used exomePeak Meng et al., 2014 for peak calling for both stress and basal condition or cortisol and mock condition, computing the enrichment of m⁶A-Seq signal in relation to the baseline RNA-Seq (input) signal at a given position, for mouse samples using high-confidence peaks called across all biological replicates and working in transcriptome space. Peaks were annotated using ChIPpeakAnno Zhu et al., 2010. Upon inspection, condition-unique peaks had some enrichment in the m⁶A-Seq in both stress- and basal conditions arguing for condition-unique peaks being rather caused by peak-detection thresholds than being true present/absent peaks. Therefore, we continued our analyses merging both peak sets (GRanges Lawrence et al., 2013). Differential m⁶A expression was analysed using DRME Liu et al., 2016, a tool specifically developed to assess methylation changes on RNA. IgG control libraries had a low percentage of alignment with no specific enrichments present. Differential gene expression was evaluated using the EdgeR package Robinson et al., 2010. Calculations and plots were done using R R Development Core Team, 2011 and ggplot2 Wickham, 2009. For mouse samples, only genes with an average estimated RPKM>1 in the input samples (Reads Per Kilobase of transcript per Million mapped reads; estimated by EdgeR rpkm function using the longest annotated isoform from mm10) were considered for both DGE (13,090 of 16,031 detected expressed genes) and m⁶A-peak calling (25,821 mapping to 11,534 genes of unfiltered 28,725 detected peaks mapping to 13,536 genes). Both peaks and genes expressed were considered differential with an absolute fold change > 0.1 and a Benjamini-Hochberg corrected P-value < 0.05 for mouse samples and fold change > 1 and a Benjamini-Hochberg corrected P-value < 0.05 for human samples. Distribution-plots of m⁶A across the transcript length were evaluated using the Guitar plots Cui et al., 2016 package. Our data revealed a higher peak abundancy within the 5’UTR than previously reported which may be due to higher representation of full length 5’UTRs in our input RNA compared to earlier studies Meyer et al., 2012 and a high contribution of m⁶Am peaks Linder et al., 2015. GO-term overrepresentation was calculated using the PANTHER Overrepresentation Test Mi et al., 2013 for “GO biological process complete” and “PANTHER Pathways” with the list of all detected genes as background. Motif search was performed by DREME Bailey, 2011 and CentriMo Bailey & Machanick, 2012 using all detected m⁶A peaks as input (mouse: 200nt sequences centred on peak summit). Motif across the whole manuscript are presented with R=A/G, W=A/T, K=G/T, B=G/C/T, H=A/C/T, D=A/G/T, V=A/G/C. For comparison of the detected mouse m⁶A-Seq GGACWB motif with known motifs, we employed Tomtom Gupta et al., 2007Ray et al., 2013 and CentriMo Bailey & Machanick, 2012. Comparison of peaks to known m⁶A was done using m⁶A data from RMBase Sun et al., 2016 (data set of 2015-10-20, GSE29714 mouse brain Meyer et al., 2012, GSE71154 mouse brain and liver Ke et al., 2015, SRA280261 mouse embryonic fibroblasts Zhou et al., 2015, GSE63753 mouse liver Linder et al., 2015) using GRanges Lawrence et al., 2013. Peak position was annotated with Biomart Smedley et al., 2015 and ChIPpeakAnno Zhu et al., 2010. Peaks on 5’UTR-CDS border and CDS-3’UTR border are counted as 5’UTR and 3’UTR, respectively. Sequencing tracks were visualized with the UCSC browser Kent et al., 2002. Overlap with existing m⁶A-reader and FMRP/FMR1 PAR-CLIP and HITS-CLIP data was done using data sets for: YTHDF1 Wang et al., 2015, YTHDF2 Wang et al., 2014, YTHDF3 Shi et al., 2017, YTHDC1 Xu et al., 2014, HNRNPC Liu et al., 2015, HNRNPA2 B1 Goodarzi et al., 2012, FMRP Ascano et al., 2012, FMR1 Darnell et al., 2011. To compare human and mouse m⁶A-peaks and binding sites gene symbols were used, first intersecting all datasets with a background list consisting of only human homologues with genes expressed in our m⁶A-Seq Smedley et al., 2015 and genes expressed in the cell lines from which the CLIP data is derived from (Higareda-Almaraz et al., 2013; Sultan et al., 2008). For random models, an appropriate number of genes were randomly selected from this background list.

6.4.2mRNA-Seq¶

Sequencing data quality control was performed by FASTQC Andrews, 2010. Reads were quality filtered (Q 20) and adapter trimmed using Cutadapt (Martin, 2011) and pseudoaligned to the mouse transcriptome (Gencode M15 transcripts GRCm38.p5) using kallisto (Bray et al., 2016; using an index with k-mer length 31, paired end strand-specific alignment, sequence based bias correction and 100 bootstraps). Differential gene expression analysis was performed using DESeq2 (Love et al., 2014) after importing the counts and summarizing to gene level with tximport (Soneson et al., 2015; considering genes with a minimum of 25 raw counts in all samples only). ENSMUSG00000019768 and ENSMUSG00000037984 were excluded from the differential gene expression analysis. For analysis of basal expression patterns, as presented in Figure 4, only the subset of unstressed “Box” animals was used.

6.4.3Gene expression¶

Statistics were performed on log2 normalized data using a 2x2 MANOVA in SPSS and were multiple testing-corrected by the Benjamini-Hochberg test (cut-off Q<0.05) in R and a cut-off by effect size (η²> 0.01) and post-hoc testing (Tukey HSD).

6.4.4Candidate m⁶A-RIP-qPCR¶

RNA abundance levels were quantified from the input samples using the ddCT method normalizing to the average of all housekeeping genes. Immunoprecipitation efficiency for each biological sample was assessed using the measured abundance of spike-in per m⁶A-RIP/IgG control sample. Because all conditions per sample were equally split at the beginning, % methylation or IgG signal was calculated as follows:

Statistics were performed on log2 normalized data (RNA) or absolute values (m⁶A) using a 2x2 MANOVA in SPSS with multiple testing performed using the Benjamini-Hochberg method (cut-off Q<0.05) in R and a cut-off by effect size (η²> 0.01) and post-hoc testing (Tukey HSD).

6.5Statisticalanalysis¶

Statistical tests were performed using SPSS (IBM SPSS Statistics, Armonk, NY: IBM Corp.) and R R Development Core Team, 2011 as indicated in Figure legends with n and statistical results indicated in Figure legends und Supplemental Tables. Plots were produced with R R Development Core Team, 2011 ggplot2 Wickham, 2009 with definition of presented measurements indicated in the Figure legends. For animal experiments, sample size was estimated a priori using G*Power Faul et al., 2007 using α = 0.05 and an experience-based β. Animals and samples within experiments were randomized using stratified randomization assisted by random number generation.

7Suplementary Information¶

7.1Suplementary Tables¶

Table S1: m⁶A peaks. Related to Figures 1 and S1.

Table S2: Statistics (GLMs & ANOVAs). Related to Figures 1-3 and S1-3.

Table S3: mRNA-Seq of Mettl3-cKO and Fto-cKO mice. Related to Figures 4-5 and S4-5.

7.2Suplementary Figures¶

Figure S1-S6