Sex-specific transcriptomic and epitranscriptomic signatures of PTSD-like fear acquisition

Summary Our understanding of the molecular pathology of posttraumatic stress disorder (PTSD) is evolving due to advances in sequencing technologies. With the recent emergence of Oxford Nanopore direct RNA-seq (dRNA-seq), it is now also possible to interrogate diverse RNA modifications, collectively known as the “epitranscriptome.”. Here, we present our analyses of the male and female mouse amygdala transcriptome and epitranscriptome, obtained using parallel Illumina RNA-seq and Oxford Nanopore dRNA-seq, associated with the acquisition of PTSD-like fear induced by Pavlovian cued-fear conditioning. We report significant sex-specific differences in the amygdala transcriptional response during fear acquisition and a range of shared and dimorphic epitranscriptomic signatures. Differential RNA modifications are enriched among mRNA transcripts associated with neurotransmitter regulation and mitochondrial function, many of which have been previously implicated in PTSD. Very few differentially modified transcripts are also differentially expressed, suggesting an influential, expression-independent role for epitranscriptional regulation in PTSD-like fear acquisition.


INTRODUCTION
Posttraumatic stress disorder (PTSD) is rooted in the acquisition of memories associated with specific traumatic events (Maren and Holmes, 2016). Evidence from clinical imaging research and studies using animal models of aspects of PTSD suggests that PTSD pathophysiology manifests from a dysregulation of the amygdala in response to trauma-related or emotional cues (Verbitsky et al., 2020;Del Casale et al., 2022;Johansen et al., 2011;Leite et al., 2022;Gonda et al., 2022;Yehuda et al., 2015).
Accumulating evidence indicates that epigenetic (dys)regulation could play a pivotal role in the pathogenesis of PTSD (Zannas et al., 2015;Howie et al., 2019). However, over the last decade, new research has shifted attention from the epigenome to the role of RNA modifications, collectively known as the epitranscriptome, in modulating brain function and behavior (Zhao et al., 2017;Madugalle et al., 2020;Livneh et al., 2020;Jonkhout et al., 2017;Murakami and Jaffrey, 2022). These modifications may impact RNA stability, localization, and splicing and regulate translation in a dynamic and reversible manner. As such, they have the potential to rapidly fine-tune responses to external stimuli and may underpin aspects of fear learning and memory (Widagdo et al., 2016;Walters et al., 2017;Engel et al., 2018). However, relatively little is currently known about the epitranscriptome, and its regulatory roles in the brain, due to technological limitations.
The recent emergence of nanopore sequencing technology enables direct sequencing of full-length native RNA molecules, circumventing the need for reverse transcription or amplification (Deamer et al., 2016). As an RNA molecule passes through a protein pore embedded in a synthetic membrane, there is a disruption in the flow of ions across the membrane, in a predictable manner depending on the composition of the bases within the pore at any given moment. This technology enables not only sequencing of full-length RNA molecules but also the detection of RNA modifications, which induce discernible shifts in current intensity. Therefore, direct RNA nanopore sequencing may be an ideal tool for comprehensive profiling of the epitranscriptome in different biological contexts (Ameur et al., 2019).
Here we employed a well-validated mouse model of aspects of PTSD to investigate transcriptional changes in the amygdala during the acquisition of conditioned fear response ( Figure 1A) ( (C) Face validity of the cued fear conditioning paradigm was verified in male and female mice, which showed strong expression of conditioned fear responses. ANOVA, ***p < 0.001 relative to non-fear-conditioned controls, **p < 0.01 for male-versus female-conditioned mice. Data are represented as mean G SEM (n = 7-10/group).  Kao et al., 2016;Assareh et al., 2020;Brzozowska et al., 2017;Dedic et al., 2018). Because we were interested in capturing the earliest molecular changes, we collected tissue for RNA extraction just 15 min after exposure to the fear conditioning cue, consistent with previous studies (Rosen et al., 1998;Huff et al., 2006). We leveraged parallel short-read RNA-seq (Illumina NextSeq) to profile gene expression and long-read Oxford Nanopore Technologies (ONT) direct RNA-seq (ONT PromethION) to build a comprehensive catalog of RNA modifications associated with conditioned fear response ( Figure 1B).

RESULTS
Our PTSD mouse model was developed using a modified auditory, cued-fear conditioning protocol developed by Siegmund and Wotjak (2007) and has been validated for its face, construct, and predictive validities in male mice (Kao et al., 2016;Assareh et al., 2020;Dedic et al., 2018;Smith et al., 2019;Brzozowska et al., 2017;Li et al., 2005). Accordingly, for our first experiment, we investigated if cued-fear memory was also observable in adult female mice 24 h after fear conditioning. We observed a significant increase in freezing behavior in both male and female conditioned (CS tone, shocked) mice compared with nonconditioned (CS tone, nonshocked) controls, with females freezing significantly more than males ( Figure 1C). Our findings were consistent with one other study that also reported greater cued-fear expression in females than males, suggestive of sex differences in discriminative ability to the cue (ter Horst et al., 2012).
These findings led us to investigate whether this significant sex difference in expression of PTSD-like fear memory could be associated with differences in transcriptional and epitranscriptional processes associated with the initial acquisition of the fear memory. To address this, a separate group of mice was sacrificed within 15 min of the fear conditioning procedure, capturing the earliest molecular signatures of fear acquisition ( Figure 1A). mRNA was extracted and isolated from the amygdala of each individual, pooled in groups of 3-4 individuals per condition, and analyzed by short-read RNA-seq to assess differential gene expression responses to the PTSD-like stressor ( Figure 1B). We obtained on average 37.9 million read pairs per sample (SD = 7.23 x 10 6 ; Table S1). In each library, most of the reads were successfully mapped to the reference genome (84.72% G 0.56) and subsequently assigned to an annotated gene (65.83% G 1.03; Table S2).
We first determined differences in gene expression in male and female mice exposed to the CS tone only (i.e., the control group). These mice did experience a novel environment, and their gene expression profile should therefore be considered not as ''baseline'' but as mildly stressed. There were 129 differentially expressed genes in the amygdala of females compared with males, including a number of genes previously identified as sexually dimorphic (Figures S1A and S1B) (Yang et al., 2006;Deschepper, 2020). Importantly, 5.1% of dimorphic expression changes observed within this female versus male control comparison were also observed in response to fear conditioning (Figures S1C, Tables S3, and S4).
We identified 190 upregulated and 401 downregulated genes in conditioned male mice (fold-change > 0.2 and FDR <0.05; Figures 1D and Table S5). We compared our list of differentially expressed genes in males with those found in another study that examined the amygdala transcriptome 2 h following auditory fear conditioning (Lori et al., 2019). There were 53 out of the 591 genes that overlapped; despite differences in methodology (e.g., foot-shock current 0.6 mA versus 1.5 mA), we speculate that changes in the expression of these genes may reflect a general, temporally independent response to the fear conditioning paradigm. In males, the expression of the top 10 differentially expressed genes (DEGs) included Zinc finger and BTB domain containing 16 (Zbtb16), LIM homeobox protein 9 (Lhx9), fibrinogen C domain containing 1 (Fibcd1), serum/glucocorticoid regulated kinase 1 (Sgk1), and Isl1 transcription factor, LIM/homeodomain (Isl1), all of which were qPCR validated using the same samples that were sequenced ( Figures 1E and 1F).
Zbtb16 was found to be significantly upregulated in conditioned male mice and is a transcriptional regulator involved in a myriad of processes including those linked to neurodevelopment (Lin et al., 2019;Gaber et al., 2013). Zbtb16 knockout mice show several behavioral impairments relevant to neurodevelopmental iScience Article disorders such as autism spectrum disorder and schizophrenia. In addition to impaired cortical morphology, these mice also show dysregulation in genes associated with myelination and neurogenesis (Usui et al., 2021). Functional expression of this gene in the adult hypothalamus has been recently linked to regulation of metabolism (Cheng et al., 2020). Lhx9 and Isl1 are transcription factors that also have roles in neurodevelopment (Garcia-Calero and Puelles, 2021; Zhuang et al., 2013). Fibcd1 belongs to a class of proteins known as fibrinogen-related domains (FreD) with multiple functions that include innate immunity and are expressed in glial cells in the human brain (von Huth et al., 2018). To our knowledge, the functional expression of these four genes in the adult amygdala and in response to PTSD paradigms such as the one used in our study has heretofore not been characterized. They therefore represent interesting targets for future investigations.
Sgk1 was downregulated in expression. It is involved in several functions, which include regulating ion channel activity, interfering with transcription, and neuronal excitability (Lang et al., 2006). In the amygdala, an earlier microarray study reported an increase in Sgk1 expression following cued-fear conditioning, which contrasts with our findings; however, these experiments were carried out using a different conditioning protocol, and analysis of amygdala RNA was at later time points (Andero et al., 2014). It is worth noting that Sgk1 is also expressed in other brain regions including the hippocampus, where it has a role in contextual fear memory formation (Lee et al., 2007) and in the prefrontal cortex (PFC). A downregulation in expression has been reported in postmortem tissue analyses of PFC tissue from patients who had suffered from PTSD and is associated with increased expression of contextual fear memory in rats (Licznerski et al., 2015).
Independent Gene Set Enrichment Analysis (GSEA) for conditioned male mice indicated a significant overrepresentation of genes implicated in synaptic activity (Table S6), a process that is well known to be implicated in the acquisition of memories (Takeuchi et al., 2013;Silva, 2003).
Epidemiological evidence suggests that females are at higher risk (2:1) of developing PTSD following trauma with sex-specific central and peripheral patho-transcriptional signatures (Breen et al., 2018;Christiansen and Berke, 2020;Girgenti et al., 2021). However, to our knowledge, few PTSD-stress paradigms, including the one we have employed here, have been applied in females (Verbitsky et al., 2020). To our surprise, despite the increased expression of cued-fear memory, fear conditioning induced a modest transcriptional response in conditioned female mice. We observed a significant decrease in expression for only a single gene, Rps27A; however, this was not validated using qPCR. There remains debate as to whether the amygdala is a sexually dimorphic structure (Marwha et al., 2017;Wolf et al., 2021;Chatzinakos et al., 2021;Montalvo-Ortiz et al., 2022). Evidence suggests there are intrinsic differences linked to sex hormones that may explain our sex-specific observations in the formation and expression of emotional memories (Blume et al., 2017;Morris et al., 2008;Cooke and Woolley, 2005). Overall, our gene expression data suggest significant distinctions exist during the early transcriptional events of fear acquisition in males and females, warranting much further investigation.
To explore epitranscriptome dynamics during fear acquisition, we performed dRNA sequencing on the same amygdala RNA samples as described above. Each sample was sequenced on a single PromethION flow cell, yielding a median of 3.9 million QC-passed reads per sample, each of which represents a native mRNA transcript (Table S7). This enabled a broad survey of expressed transcripts (average read count >3), with 12,792 out of 32,604 (39.2%%) of protein coding genes reaching an average read count greater than 15, the minimum threshold required for detection of putative RNA modifications (Table S8).
An overview of our analysis workflow used to profile RNA modifications is shown in Figure 2A, exemplified using the transcript for Neurensin-2 (Nrsn2), a known modulator of emotional behavior and putative biomarker for PTSD (Umschweif et al., 2021;Glatt et al., 2013). Briefly, we used a combination of two recently developed software packages, Xpore (Pratanwanich et al., 2021) and Nanocompore (Leger  Figure 2B). Nanocompore identified 43,172 candidate modification sites across 1,547 transcripts for males and 33,216 candidates across 1,482 genes for females ( Figure 2B). The tools use alternative metrics and statistical approaches to identify modified k-mers, and concordant results provide the best evidence of modified sites. To maximize stringency, we retained only those candidates who were identified by both Nanocompore and Xpore and collapsed any candidates within 10 nt into a single modification site ( Figure 2A). For males, this approach identified 7,487 candidate modifications between conditioned versus control conditions that were collapsed into 3,397 unique differentially modified sites ( Figure 2B). In females, 4,808 candidates were collapsed into 2,331 unique differentially modified sites ( Figure 2B). These high confidence sites showed mean 1.7-fold higher modification rates, as measured by displacement of the current signal, under the stress condition in both males (stress = 0.48, control = 0.27) and females (stress = 0.49, control = 0.29; Figure 2C).
Most high-confidence modified sites were found within 3 0 untranslated regions (UTR) of mRNA transcripts, with an enrichment close to the boundary with the coding sequence (CDS) observed in both males and females ( Figures 2D and 2E). A smaller proportion of modified sites were found within 5 0 UTRs, which are generally shorter than 3 0 UTRs, but 5 0 UTRs exhibited the highest density of sites per nucleotide compared with the 3 0 UTR and CDS in both sexes ( Figures S1D and S1E). These observations are consistent with positional distributions previously reported for m6A and other known RNA modifications (Meyer et al., 2012(Meyer et al., , 2015Dominissini et al., 2012).
Although Nanocompore and Xpore can identify differentially modified bases, neither tool determines the specific type of RNA modification at a given site. To further characterize and validate our modification sites, we assessed the enrichment of sequence k-mers consistent with the canonical m6A DRACH motif (Wang et al., 2021a). We found an enrichment of DRACH motif k-mers within high-confidence modification sites, compared with 10 sets of random nonmodified sites, and no enrichment for k-mers of a randomly selected motif ( Figure 2F). We also intersected our set of high-confidence modification sites with RMBase, a database of known RNA modifications (Xuan et al., 2017). We found 506 annotated RMBase modifications overlapping differentially modified sites in males and 421 in females, with almost all of them identified as m6A ( Figure 2B). We retrieved the consensus sequence around these m6A sites and consistently observed the DRACH motif for both sexes ( Figure 2G). Concordance with previously annotated RNA modifications and adherence to known patterns of transcript position and sequence context confirms the reliability of our catalog of differentially modified sites during fear acquisition.
Next, we functionally characterized genes and transcripts identified as being differentially modified. Out of 3,397 and 2,331 modified sites detected in males and females (Tables S9 and S10), respectively, we found 580 common modifications in 186 genes (average = 3.11 sites/gene; Figure 2G; Figure S1F), suggesting a degree of both shared and sexually dimorphic epitranscriptome dynamics. GSEA for males and females indicated enrichment of equivalent and/or related gene ontology terms ( Figure 2I). Within the top 20 most significantly enriched terms, 60% were found in both males and females (Tables S11, S12, and S13). iScience Article Mitochondrial function and neurotransmitter production/secretion were the predominant functional categories identified by GSEA ( Figure 2I). The former is particularly interesting in the context of evidence implicating mitochondrial function in the stress response (Daniels et al., 2020;Manoli et al., 2007;Picard et al., 2018) and mitochondrial dysfunction in PTSD pathology (Zhang et al., 2015;Flaquer et al., 2015). Mitochondrial metabolic processes are components of the cell danger response (CDR) that is activated by stressors. Negative feedback mechanisms exist to ''switch off'' the CDR once the stress has passed but if they fail, this may contribute to oxidative stress and inflammation (Naviaux, 2014).
Among the top differentially modified gene candidates were multiple synaptic genes previously implicated in PTSD, including neuropeptide Y (Npy) (Sah and Geracioti, 2013). Evidence suggests the Npy system promotes resilience-to-stress in rodents and reduced NPY in humans is associated with PTSD (Sah and Geracioti, 2013). We observed differentially modified sites in the 3 0 UTR and CDS regions of Npy in males and females. In each case, the modification frequency was upregulated during fear conditioning, whereas expression of Npy mRNA was not found to be upregulated, suggesting transcription-independent regulation of Npy by RNA modifications may occur during fear conditioning. The lack of concordance between differential modification and expression in Npy was similarly observed across the transcriptome, with just 150/3,397 (4.41%) of high-confidence modifications being harbored within 42 differentially expressed genes ( Figure 2J; Table S14).

DISCUSSION
Gaining an insight into the molecular mechanisms that are triggered immediately following trauma is important to comprehend how they are implicated in the neuropathological cascades that culminate in the onset of PTSD. Our study provides a framework for exploring acquisition-induced transcriptomic and epitranscriptomic events underpinning PTSD-like phenotypes in mice. Our results show that any effect of RNA modifications does not result in obvious gene expression changes at an early stage of fear conditioning. Therefore, it is possible that RNA modifications do not have any role at all in that process, but their enrichment in neurological/PTSD-related genes is an interesting finding that warrants further investigation. There is broad scope for future work, but one major issue must be addressed: how does one interrogate the functional roles of RNA modifications (Schaefer et al., 2017).
For some modifications such as m6A, the majority of known m6A sites are unmethylated at baseline (He and He, 2021). The functional relevance of constitutive versus regulated m6A sites is unknown but such sub-stoichiometric levels would indicate a large margin in which to regulate RNA metabolism and gene expression. Furthermore, although the life cycle of modifications is thought to be transient, this may vary from transcript to transcript. Our data are also consistent with recent work showing that RNA modifications can form clusters on the same transcript , exerting cumulative or combinatorial effects on RNA metabolism. One low-resolution approach to providing answers is by manipulating the expression and activity of epitranscriptomic machinery enzymes-''writers,'' ''erasers,'' and ''readers.'' This has been used by numerous studies in the context of learning, memory, and the stress response (Walters et al., 2017;Widagdo et al., 2016;Engel et al., 2018). There have also been several pharmacological compounds developed to inhibit or activate these enzymes Selberg et al., 2021). One higher resolution strategy involves the use of CRISPR-Cas-inspired RNA-targeting (CIRT) system (Rauch et al., 2018(Rauch et al., , 2019 but this is yet to validated in vivo. In conclusion, our study opens a new and exciting frontier in molecular psychiatry research that has potential to reshape our understanding of stress-induced psychiatric disorders. RNA modifications could be among the earliest molecular events associated with fear conditioning, but more studies are needed at additional time points, using proteomics and ribosome profiling among other techniques to disentangle the many possible (and nonmutually exclusive) impacts of specific m6A sites to the pathophysiology of PTSD (Engel et al., 2018). RNA-based therapies have found success in treatment of diseases from cancer to COVID19 (Damase et al., 2021). Perhaps, there is scope to treat psychiatric disorders such as PTSD in a similar manner.

Limitations of the study
Each sample analyzed in our study was a pool of 3-4 independent individuals, meaning a total of 12 different individuals per condition were represented. Sequencing each individual as a separate biological replicate rather than pooling would increase the statistical power of our study; however, due to the high cost and large RNA input quantity required for direct-RNA sequencing, adding more samples to our study was unfeasible.

OPEN ACCESS
iScience 25, 104861, September 16, 2022 7 iScience Article Sample multiplexing is not supported for direct-RNA sequencing, so each sample must be sequenced on a separate ONT PromethION flow cell, resulting in prohibitively expensive costs for large experimental designs. Furthermore, although we were focused on capturing the earliest transcriptional changes associated with fear conditioning, additional time points will illuminate the temporal relationships between transcript abundance and transcript modification to the pathophysiology of PTSD. In addition, the use of proteomics techniques, such as ribosome profiling, will enable the understanding of the real impact of transcriptomic and epitranscriptomic changes to the acquisition of traumatic memories at the molecular level.

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

ACKNOWLEDGMENTS
AG wishes to thank Dr. Nicholas Everett and Dr. Erin Lynch for technical assistance and the Laboratory Animal Services team at the Brain & Mind Center for assistance with animal husbandry. We thank our colleagues Derrick Lin and Warren Kaplan for providing technical support and oversight of high-performance computing at the Garvan Institute. AG also wishes to thank Professor Tim Bredy (University of Queensland) for the many intellectually challenging conversations associated with this study.
We would also like to acknowledge that AG was supported by a University of Sydney Research Fellowship.

DECLARATION OF INTERESTS
IWD manages a fee-for-service nanopore sequencing facility at the Garvan Institute of Medical Research that is a customer of Oxford Nanopore Technologies but has no further financial relationship. The authors declare no other competing financial or nonfinancial interests. iScience Article Data and code availability Raw sequencing data has been uploaded to the NCBI SRA (PRJNA779821) and are publicly available as of the date of publication. Accession numbers are listed in the key resources table.
All original code has been deposited at Zenodo and is publicly available as of the date of publication. DOIs are listed in the key resources table.
Any additional information required to re-analyse the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
The study used 22 adult male and 24 adult female (10 weeks of age) C57BL/6JArc mice purchased from the Animal Resource Centre, Western Australia. All experiments with mice conformed to the regulatory standards approved by the University of Sydney, Animal Ethics Committee (2018/1425). As per the approved ethics protocol, all animals were initially group-housed for 1 week followed by single housing in individually ventilated cages with environmental enrichment and ad-libitum access to food and water. Cages were in racks in a room maintained in a temperature/humidity-controlled environment (21 C, 55.5%) with reverse-phase lighting (lights ON: 21h00, OFF: 09h00). In our experiments, we included the use of randomly cycling female mice (Shansky and Woolley, 2016).

METHOD DETAILS
Cued fear conditioning and cued fear memory testing Mice were fear conditioned in a wooden sound-proof chamber, transparent front and rear-facing walls, opaque side walls with a white light and a metal grid floor (Context A). The conditioning session started with a 198 s acclimation period which was followed by a 30 s 80dB tone (9KHz sinewave, conditioned stimulus, CS) that co-terminated with a 2s 1.5mA foot-shock delivered through the metal grid as a constant current (unconditioned stimulus, US). The conditioning session ended with a 60 s after-shock interval. Nonshock control mice were exposed only to the CS tone. The conditioning session from start to finish is just under 5 min. Fear conditioning chambers were cleaned using 70% ethanol spray solution. Within 15 mins of being fear conditioned, one group of stress and control mice were culled, brains extracted and snap frozen. Cued fear memory was tested in another group of stress and control mice 24hrs later in a novel context (Context B: wooden sound-proof chamber, red light, red perspex flooring, peppermint oil scent). Following a 2 min baseline period, all mice were presented with 3 x 30 sec CS with an inter-tone interval of 1 min. The memory retrieval session ended with a 60 s interval after the last tone was presented. Chambers were cleaned using F10 disinfectant spray solutions. Freezing was measured in response to the CS using the CleverSys FreezeScanâ video tracking system and software.

Tissue collection & RNA extraction
Mice were sacrificed by cervical decapitation within 15 min of being fear conditioned. Brains were extracted, snap-frozen in liquid nitrogen and stored at -80 C. Bilateral punches (0.5 to 1mm) of the amygdala (AMG: central and basolateral amygdala) were made using the Paxinos and Watson Atlas as a guide (Paxinos, 2004;Zapala et al., 2005). Tissue was collected in RNA-free tubes and RNA was extracted using Qiagen RNeasy TM Micro kit (Qiagen, MD, USA) according to the manufacturer's instructions. RNA concentrations were quantified using a NanoDropä One spectrophotometer (Thermofisher Scientificâ, MA, USA) and only samples with 260/280 ratios of greater than 1.7 were used for downstream analyses. Samples from 3-4 mice were randomly pooled for subsequent analyses by Illumina NextSeq and Promethion, meaning that a total of 12 different individuals were represented in each condition. mRNA was purified from total RNA pools using the Dynabeadsâ mRNA Purification Kit (Thermofisher Scientificâ, MA, USA) according to the manufacturer's instructions.

Bioinformatics -gene expression profiling
Reads were trimmed using Trim Galore (Galaxy version 0.6.3) and subsequently aligned to the reference Mus musculus genome (GRCm38) using HISAT2 (Galaxy version 2.1.0). Reads were annotated (Gencode v.M25) and counted using featureCounts (Galaxy version 2.0.1). Biostatistics and visualization were run in R (version 3.6.3) with the Rstudio GUI (version 1.2.5033). For each comparison -male: stress(3) v control(3), female:stress(3) v control(3) -genes were filtered to have at least 5 reads in at least 2 samples for each gene and have gene biotype of protein-coding, long non-coding or microRNA. RUVSeq (Risso et al., 2014) and DESeq2 (Love et al,. 2014) were then used to quantify differential gene expression between groups. We set a threshold criteria of a minimum 0.2 expression fold-change and FDR<0.05. Gene ontology and enrichment analysis was performed using EnrichR with a significance threshold set at FDR <0.05 .

Bioinformatics -Epitranscriptome profiling
To generate a comprehensive catalogue of high-confidence RNA modification sites present in the amygdala of male and female mice exposed to conditioned fear, we used two recently released software Nanocompore and Xpore on our direct RNA-seq libraries. Nanocompore collects current intensity and dwell time at each position, as the native RNA goes through the pore, and uses those variables to perform pairwise statistical comparisons between control and stress samples. The tool reports positions with a significant statistical difference for each of the individual variables, but also models a logistic regression that takes both current intensity and dwell time into account in search for significant differences. To minimise the false-positive discovery rate, we applied strict filtering to the results, only retaining putative modifications with a significant p-value (alpha=0.01) in all three statistical tests. Xpore, on the other hand, models the current intensity as a mixture of two Gaussian distributions, for modified and unmodified RNA, using prior information to help guide parameter estimation. Then given two conditions, such as control and stress, the tool is able to identify differential modification, but also quantitate them. We only selected putative modifications with significant differences in modification rates (t-test, p-value < 0.01) between the conditions and that were increasingly modified in the stress mice. To further increase the confidence in our dataset, we only retained modifications identified by both Nanocompore and Xpore, and that were at most 10 nt apart from each other and we also collapsed modifications less than 10 nt apart into a single modification site. We compared the coordinates of the modification sites identified using this approach against modification sites previously annotated in RMBase v2.06. We also performed gene set enrichment analysis using EnrichR as described above.

Statistical analysis
Cued-fear expression data was analysed using ANOVA with significance threshold set at p<0.05 (**p<0.01 & ***p<0.001). For all bioinformatics analyses (see supplemental information), a significance threshold was set at FDR <0.05.

ll
OPEN ACCESS