Dynamic modulation of N6-methyladenosine by ionizing radiation in human cells

TL;DR

This study shows that ionizing radiation rapidly and dynamically changes RNA methylation in human cells, with potential implications for cancer treatment and radiation protection.

Contribution

The study reveals a dynamic, time-dependent modulation of m6A RNA methylation in response to ionizing radiation in human cells.

Findings

m6A methylation peaks one minute after radiation exposure and then fluctuates over time.

Two UQCR10 transcripts show stable hypermethylation following radiation exposure.

Disruption of m6A enzymes reduces RNA methylation in response to radiation.

Abstract

A cell's transcriptome is regulated through the integration of external and internal signals that activate intracellular signal pathways, epigenetic modifications and post-translational changes. Post-transcriptional regulation through RNA methylation has emerged as an important mechanism in cancer development, and informative for diagnosis and treatment. The most abundant one, N6-methyladenosine (m6A), regulates gene expression in eukaryotes. In the present study m6A RNA modifications have been characterized in response to ionizing radiation (IR) exposure in the HT1080 human cell line. Cells were exposed to a dose of 10 Gy of X-rays and harvested 1, 2, 10 min, 1 and 24 h after exposure. m6A sites were identified using long read nanopore direct RNA sequencing. A pipeline was designed using m6Anet to estimate m6A stoichiometries transcriptome-wide, which were then analysed by a beta-binomial regression model with moderated dispersion estimates and independent filtering to detect differentially methylated (DM) sites between treated and control samples. We found that IR modifies m6A sites in a dynamic way, inducing site specific increase of methylation. Remarkably, it peaks within the first minute after exposure, followed by a sharp decrease at 1 h without returning to baseline, increasing again after 24 h. Two transcripts of the nuclear encoded gene UQCR10, a subunit of the respiratory chain protein, sharing the same site presented a stable hypermethylation over time, confirmed by a modified quantitative PCR assay. Moreover, we generated Knockouts (KO) cell lines for 3 key enzymes involved in m6A methylation, a writer, a reader and an eraser namely METTL3, YTHDF2 and FTO, to better understand mechanistically IR driven m6A dynamics. Importantly, all three KOs presented a transcriptome wide decrease in RNA methylation following IR exposure. Lastly, m6A modifications were also confirmed in human skin biopsies exposed to IR, with the UQCR10 gene site also hypermethylated 24 h after a lower 2 Gy X-rays dose. To summarise, we provide evidence that IR modulates RNA m6A levels in a site-specific and dynamic way, with DM sites enriched in genes involved in bioenergetics, cell signalling/migration and apoptosis pathways, thus representing a rapid cellular response to radiation. Considering the essential role of m6A in controlling gene expression and physiological activities, this study established the basis for further studies assessing IR driving m6A with a potential role in radiation oncology and protection.

•IR dynamically modifies m6A sites, causing site-specific increases in methylation.•m6A methylation peaks at 1 min, drops at 1 h, then rises again by 24 h post-exposure.•Two UQCR10 transcripts sharing a site showed stable hypermethylation over time.•DM sites are enriched in genes of bioenergetic, signaling, apoptotic pathways.•IR-mediated RNA methylation responses require functional writers, readers and erasers.

IR dynamically modifies m6A sites, causing site-specific increases in methylation.

m6A methylation peaks at 1 min, drops at 1 h, then rises again by 24 h post-exposure.

Two UQCR10 transcripts sharing a site showed stable hypermethylation over time.

DM sites are enriched in genes of bioenergetic, signaling, apoptotic pathways.

IR-mediated RNA methylation responses require functional writers, readers and erasers.

Figures

Landscape of differentially methylated (DM) m^6^A sites after X-ray exposure. HT1080 cells were irradiated with 10 Gy X-rays and sampled at 1 min, 2 min, 10 min, 1 h, and 24 h (n = 3 per timepoint). A – E) Volcano plots of per-site differential methylation. Each dot is a candidate m^6^A site; DM sites (FDR <0.1, Benjamini–Hochberg) are red. 10 Gy: change in m^6^A stoichiometry (treated − control). y-axis: −log10(*p* value). Statistics from moderated beta-binomial regression on modified vs unmodified read counts derived from m6Anet. F) Number of DM sites detected at each timepoint. G – I) Aggre

Selection and validation of persistent DM m^6^A sites following X-ray exposure. A) Venn diagram showing the significant DM sites at each timepoint (FDR <0.05). Four sites persistently significant across all timepoints, two of which are shown in B) Genomic location of m^6^A site persistently hypermethylated across all timepoints, found in two distinct transcript isoforms of the gene *UQCR10*. Red boxes: exons; arrows: transcript orientation; vertical blue line: shared genomic coordinate of the DM site (hg38, chr22). C – D) Methylation dynamics for the site in each isoform, both of which belong

Overview of m^6^A response to IR in the HT1080 KO cell lines exposed to 10 Gy and sampled at 1 h and 24 h (n = 3 per group). WT shown for comparison purposes. A) to D) Volcano plots at 1 h for WT, and the 3 KOs FTO, METTL3, and YTHDF2 respectively. Each dot is a candidate site; red = DM (FDR <0.1). x-axis: change in stoichiometry (irradiated − matched non-irradiated control within genotype). y-axis: −log10(*p* value). Statistics from moderated beta-binomial regression on modified vs unmodified read counts derived from m6Anet. E – H) The same sampled at 24 h. I) Heatmap showing the degree of hy

Integration of transcriptomic and epitranscriptomic responses to 10 Gy X-ray exposure in HT1080 WT and KO cell lines. A) Bar plot showing the number of differentially up- (red) and down-regulated (blue) genes (DEGs) at each timepoint (1 min – 24 h) for the WT HT1080 cell line, (FDR <0.1). B) Similar to A but for the KO cell line timepoints (1 h and 24 h), WT timepoints also shown for comparison purposes. C) Relationship between Gene Ontology terms enriched in differentially methylated transcripts (p < 1e-3; columns) at any timepoint post-exposure, and differentially expressed transcripts (p <

Effect of 2 Gy X-ray exposure on m^6^A levels in human skin biopsies after 24 h post-exposure (n = 5, each paired with a non-exposed matching-donor sample). A) Volcano plot showing the site-level methylation response at 24 h post-exposure. x-axis is the change in m^6^A methylation (logit scale), y-axis is the significance (*p* value, -log10 scale). Statistics obtained using a paired moderated beta-binomial regression model on modified vs unmodified read counts derived from m6Anet. Colour indicates the m^6^A site ranking (lighter = higher rank) which was used for the site-level GSEA analysis. B

Full text

Introduction

1

One of the earliest DNA epigenetic modifications identified in humans was methylation. DNA methylation modifications were brought to light initially by Hotchkiss (1948) who identified modified cytosines [1], but it wasn't until the 1980s that these DNA modifications were linked to the regulation of gene expression and cell differentiation [2,3]. These modifications affect the DNA activity or function without changing the sequence. Gene methylation is mainly a silencing epigenetic mark and errors in this process can lead to the development of different human diseases [4]. These modifications are not exclusive to DNA, and RNA goes also through post-transcriptional methylation. m^6^A RNA modifications in eukaryotic mRNA and lncRNA were first reported in the 1970s [5]. Due to the difficulties to detect these modifications, the RNA modifications research field was kept on hold for many years [6,7]. More recently, numerous modifications have been identified such as 5-methylcytidine (m^5^C), 5-methyluridine (m^5^U), 3-methyluridine (m^3^U), N^6^,2′-O-dimethyladenosine (m^6^Am), 1-methyladenosine (m^1^A) and 1-methylguanosine (m^1^G) [8]. Their role in RNA metabolism such as nuclear export, translation, decay and alternative splicing has started to emerge [[9], [10], [11], [12], [13], [14]] following the development of m6A-seq and MeRIP-seq [7] techniques in 2012. Currently, there are several validated methods available, such as immunoprecipitation techniques, enzyme-dependent and chemical methods, but all of them have limitations [15]. Nanopore sequencing has arisen as a promising technology for detecting RNA modifications at single-nucleotide resolution, with the main advantage of being able to sequence full-length RNA molecules in their native state [16]. This technology improves the identification of sites in different isoforms, the presence of m^6^A clusters and the simultaneous presence of different modifications on the same RNA molecule [17]. With these methodological advances, the regulatory effect of RNA modifications has expanded into all the different aspects of RNA metabolism [18,19]. The function of these modifications depends on the nucleoside modified but also on the type of RNA (mRNAs, tRNAs, lncRNAs, miRNAs, circRNAs, or rRNAs) and these modifications are added to exons at pre-mRNA level before splicing [20]. Deposition can also happen within introns and it has been linked to regulation of alternative splicing events [21].The most abundant and one of the few reversible modifications in eukaryotes is m^6^A [22] and it is located in 3′ untranslated regions (UTRs) but also enriched in long exons and near stop codons [7]. m^6^A dynamic equilibrium is orchestrated by 3 types of enzymes, methyltransferases (“writers”), m^6^A demethylases (“erasers”) and m^6^A -binding proteins (“readers”). Alterations on this complex equilibrium can lead to diseases and cancer development [[23], [24], [25]]. As a result, studying m^6^A modifications is of great interest, not only mechanistically but also as emerging prognostic biomarkers for cancer detection [26].

Deregulation of m^6^A modifications can be induced by factors which affect the writers, erasers and readers such as cellular stress, post-translational modifications or miRNAs [27,28]. Genotoxic factors like IR, are known to not only damage DNA but also RNA by inducing strand breaks and oxidative damage which will lead to errors in protein synthesis or dysregulation of gene expression [[29], [30], [31]]. Recent studies have presented evidence of IR and non-ionizing ultraviolet (UV) radiation induced biological effects through RNA epitranscriptomic modifications [32] such as activation of DNA damage response [33,34], DNA repair [[35], [36], [37]] and programmed cell death [38,39]. For instance, some studies have found a correlation between transcriptional changes of the “reader” and “eraser” enzymes and radiotherapy resistance [40,41] suggesting a role of epitranscriptomics on cellular response to radiation. Besides, m^6^A modifications in radiation-responsive mRNAs have been identified in blood from mice exposed to acute doses of IR [42]. These modifications showed dose- and time-dependent effects suggesting the potential use of RNA m^6^A modifications for biodosimetry [43]. However, differences in regulation of m^6^A levels between gamma or UV irradiated in U2Os cells have been observed, presenting a lack of induction after gamma irradiation [33]. A recent study has also demonstrated that ionizing radiation decreases RNA m^6^A methylation by downregulating METTL3 expression through regulation of its promoter via PARP1 [44]. METTL3 also regulates DNA damage repair and p53 signaling pathways [45,46]. This evidence demonstrates that factors which induce the DNA damage response can have a strong effect on modulating the m^6^A mRNA patterns impacting in the mRNA biology and metabolism probably in a tissue, cell type, dose- or radiation quality-dependent manner.

In the present study, we explored the modulation of m^6^A patterns at very early time points (inside 24 h) after IR exposure in human cells. We provided a comprehensive characterization of the modified sites, along with an evaluation of the affected genes and pathways. These new findings on the regulation of m^6^A dynamics by IR establish the basis for further mechanistic studies in RNA methylation biology and in radiation biomarker research.

Material and methods

2

Cell culture and irradiation

2.1

HT1080 human cell line was obtained from ATCC (ATCC CCL-121) and maintained in Minimum Essential Medium (MEM) containing 10% FBS (fetal bovine serum) and 1% penicillin/streptomycin, at 37 °C and 5% CO2. All cell culture media were purchased from Thermo Fisher Scientific (Paisley, Scotland, UK). HT1080 cells are a human fibrosarcoma cell line derived from connective tissue from a patient with fibrosarcoma. Cell authentication was performed by the supplier (ATCC), and experiments were conducted within two years of obtaining the cells.

Cells in T-25 flasks were exposed to a 10 Gy X-rays dose (1.7 Gy/min) using A Harmonized System (HS) x-ray system (AGO X-Ray Ltd., Aldermaston, Reading RG7 4 PW, UK) (output 13 mA, 250 kV peak). Once irradiated, the samples were incubated at 37 °C with 5% CO_2_ for 24 h, 1 h, 10, 2 and 1 min. The experiment was performed 3 times from the same cell passage. After the different incubation periods, the media was replaced with cold PBS and the flasks were placed on top of ice to stop any enzymatic activity. Then, the media was replaced with lysis buffer with β-mercaptoethanol, and the lysates were collected into 15 ml tubes and stored at −80 °C until RNA extraction.

HT1080 FUCCI cell line was kindly donated by Dr James Orth from the University of Colorado [47]. This cell line was used to generate a METTL3 and YTHDF2 and HT1080 wild type (WT) for FTO KO. CRISPR Cas9 mediated knockout cell clones of METTL3, YTHDF2 and FTO in HT-1080 cells were generated by EditCo Bio, Inc. (Redwood City, CA, USA) (Table S1). Mycoplasma testing was performed by EditCo Bio upon receipt of the cell lines.

To generate these cells, ribonucleoproteins containing the Cas9 protein and synthetic chemically modified sgRNA were electroporated into the cells using EditCo's optimized protocol. Editing efficiency is assessed upon recovery, 48-72 h post electroporation. Genomic DNA is extracted from a portion of the cells, PCR amplified and sequenced using Next-Generation Sequencing (Table S2). The resulting sequencing was analysed using EditCo's ICE analysis software.

To create monoclonal cell populations, edited cell pools are seeded at 1 cell/well using a single-cell printer into 96 or 384 well plates. All wells are imaged every 3 days to ensure expansion from a single-cell clone. Clonal populations are screened and identified using the PCR-ICE-Analysis genotyping strategy described above.

The information of the mutations present in each cell line can be found in Table S1 together with the demonstration of complete loss of each specific protein by western blot (Fig. S1).

RNA extraction

2.2

RNA was extracted using the RNeasy Mini kit (Qiagen, UK) following the manufacturer instructions. The quantity of isolated RNA was determined by spectrophotometry with a ND-1000 NanoDrop (Thermo Fisher Scientific, USA) and quality was assessed using a Tapestation 2200 (Agilent Technologies, USA).

Quantitative real-time RT-PCR validation of specific m6A sites

2.3

A specific m^6^A methylated site identified with the m^6^A single nucleotide resolution analysis was validated using a qPCR based protocol described elsewhere [48]. Briefly, 150 ng of RNA was reverse transcribed (RT) with the following specific primers surrounding a methylation site on UQCR10 transcript: UQCR10 reverse positive primer: 5′ caaacacagcagaaggtaaagag 3’; UQCR10 reverse negative primer 5′ tggtccatgaagagtaccttgagca 3’. The positive primer will include the methylated site but not the negative primer. Two RT mixes were prepared for each reverse transcription enzyme: BstI mix: 50uM dNTPs, 0.1 U enzyme, 1ul 10x buffer, up to 6ul water and MRT mix: 50uM dTNPs, 0.8 U enzyme, 2ul 5X RT buffer, up to 6ul water.

After the reverse transcription, the following primers were designed for qPCR analysis UQCR10 forward 5′ actgtgaagaaacatggcgg 3′ and reverse 5′ atgatggtgagggcgaagg 3’. Quantitative real-time RT-PCR was performed in a Quantstudio^TM^ 7 Pro system (Thermo Fisher Scientific, UK) using SYBR Green master mix (A46109, Thermo Fisher Scientific, UK). The samples were run in duplicates in 10 μl reactions with 1 μl of the diluted RTs (1/2 in sterile H_2_O) with the sets of primers at 300 nM concentration each. The reactions were performed with the following cycling conditions: 95 °C 2min, 40 cycles x (95 °C-15 s, 60°C-1min) followed by the melting curve analysis using the following conditions: 95 °C-15 s, 60 °C to 95 °C with a 0.5 °C increment. The calculation of the relative m^6^A levels followed the calculations described by Olazagagoitia-Garmendia & Castellanos-Rubio [48,49].

ELISA assay

2.4

Total RNA from the samples used for sequencing (n = 3 per time point) was poly(A) + enriched using Oligo(dT) 25 Dynabeads (61005, Thermo Fisher Scientific, UK) according to manufacturer's protocol before performing the Elisa assay. The detection of relative m^6^A levels from 50 ng of poly(A) enriched RNA samples was performed using the EpiQuik m^6^A RNA Methylation Quantification Kit (P-9005-96, EpiGentek, USA) following the manufacturer's instructions.

Cell proliferation

2.5

A cell proliferation assay using the IncuCyte live-cell imaging (Sartorius, IncuCyte S3, Germany) was carried out to evaluate and compare the rate of proliferation between irradiated and non-irradiated HT1080 wild-type and KO cell lines. Cells were seeded in 6-well plates at a density of 5 × 10^4^ cells per well in α-MEM GlutaMAX medium and incubated for 24 h at 37 °C in an atmosphere containing 5% CO_2_. Each 6-well plate was arranged as such that the upper 3 wells contained the WT cells, and the lower three wells contained cells from a KO cell line. This setup was repeated separately for each of the 3 KO cell lines. The medium was refreshed after 24 h, and the cells were subsequently exposed to X-ray doses of 2 and 10 Gy. Plates were placed inside the IncuCyte system directly after irradiation, where image recording started. Images were taken every 30 min for 6 days. During the 6 days experiment, the medium had to be changed once, on day 3.

Confluence measurements over time were extracted from the images using the IncuCyte software version 2020B and imported into R [50] for all further analysis. This included filtering out timepoints with missing values, identifying the plateau phase and removing timepoints beyond it, and fitting a generalised logistic curve to model the cell growth kinetics at each condition. The plateau phase was identified in supervised automated way, by fitting a smooth spline with 15 degrees of freedom, predicting its first derivative and finding where it crosses 0, and visually inspecting the start of the plateau was correctly identified. Modelling was done using the 'nlme' package [51], where the 3 wells of each genotype were modelled as a random variable. Note that the fitting was done separately for each 6-well plate.

Western blotting

2.6

Cells were lysed in RIPA buffer (Sigma-Aldrich, UK) containing a protease inhibitor cocktail (Thermo Fisher Scientific, UK). Protein content was quantified using a BCA protein assay (Thermo Fisher Scientific, Pierce BCA Protein Assay Kits, UK). 40 μg of protein was separated on a precast 4–20% gradient polyacrylamide gels (4561094, Bio-Rad, UK) and then transferred onto a 0.2 μm PVDF membrane using a Turbo Mini PVDF Transfer packs (Bio-Rad) in a Trans-blot turbo transfer system (Bio-Rad, UK). Membranes were incubated with antibodies against METTL3 (rabbit, ab195352, secondary antibody anti Rb ab205718, Abcam, UK), YTHDF2 (rabbit, ab220163, secondary antibody anti Rb ab205718 Abcam, UK), and FTO (rabbit, ab195352, secondary antibody anti Rb ab205718, Abcam, UK). β-actin (C4 HRP SC-47778, Santa Cruz, USA) was used as loading control. Chemiluminescence from the blots was detected with a ChemiDoc imaging system (Bio-Rad, UK).

Human skin biopsies

2.7

Adult skin biopsy punches from the abdominal region were used for the experiments. Skin samples from 5 donors were obtained from plastic surgery procedures (purchased from Biopredic International, France), with consent for research use assured by the supplier. Biopredic International complies strictly with the ethical rules for donation and collection of human body products, in view of research use only, according to the French Law L.1245-2 CSP. Biopredic International holds all regulatory permits for the collection, processing, transfer and export of human body products to be used exclusively for research.

The donors were Caucasian women with ages ranging from 30 to 55, body mass index between 19 and 36, and skin phototype I and II. From the skin samples, we obtained two 4 mm biopsies per sample using a biopsy punch tool (Selles Medical, UK). Biopsies were transferred into 12 well plates with 3 ml of DMEM:F12 1:1 media and 1% of antibiotic and antimycotic solution (Sigma-Aldrich, UK) per well. Then, 5 skin biopsies were exposed to 2 Gy dose of X-rays (dose rate 0.5 Gy/min) and the other 5 were mock-irradiated (controls). After exposure, biopsies were incubated for 24 h at 37 °C with 5% CO2 and immediately after, preserved in RNA later at – 80 °C until the RNA extraction. This study was approved by the UKHSA Research Ethics and Governance Group (Reference: R&D 471).

RNA was extracted using the RNeasy fibrous kit (Qiagen, Manchester, UK) following the manufacturer instructions.

Direct RNA sequencing

2.8

HT1080 cells: independent sequencing libraries per replicate (n = 3 per time point) were prepared using 500 ng of total RNA per sample in accordance with the protocol provided for the Direct RNA Sequencing Kit SQK-RNA002 (Oxford Nanopore Technologies). Library quantification was carried out using the Qubit dsDNA High Sensitivity Assay Kit (Thermo Fisher Scientific) to ensure accurate input concentrations. Sequencing was performed on a PromethION 2 Solo device (Oxford Nanopore Technologies) with R9 flow cells, which were primed using the Flow Cell Priming Kit EXP-FLP002 (Oxford Nanopore Technologies) to optimise performance.

Data acquisition was managed with MinKNOW software, with sequencing runs performed for up to 72 h to maximise data yield and throughput.

Human skin: Sequencing libraries were prepared using 1 μg of total RNA per sample in accordance with the protocol provided for the Direct RNA Sequencing Kit SQK-RNA004 (Oxford Nanopore Technologies). Library quantification was carried out using the Qubit dsDNA High Sensitivity Assay Kit (Thermo Fisher Scientific) to ensure accurate input concentrations. Sequencing was performed on a PromethION 24 (Oxford Nanopore Technologies) with PromethION RNA flow cells (FLO-PRO004RA). Data acquisition was managed with MinKNOW software, with sequencing runs performed for up to 72 h to maximise data yield and throughput.

Bioinformatic and statistical analysis

2.9

m6A differential methylation analysis

2.9.1

m6Anet [52] was used to estimate the m^6^A modification ratio of sites across the transcriptome. The required input files were created based on the authors’ recommendations, using Samtools [53], minimap2 [54], and nanopolish [55] (Fig. S2). The resulting sites were filtered using the probability_modified metric, setting the cutoff to at least one sample >0.9. The effect of IR on the estimated m^6^A methylation level at each site was modelled using beta-binomial regression in the aod R package [56], with moderated dispersions which were calculated using the DRIMSeq package [57] (see supplementary methods for details). Multiple test correction was done using the Benjamini-Hochberg procedure [58], with a cutoff set at FDR <0.1. Independent filtering of sites was done based on mean coverage to reduce the statistical burden and thus improve sensitivity. Additional details included in the Supplementary Text (Differential Methylation Analysis).

Over-representation analysis

2.9.2

Hypergeometric tests were done to assess the over- or under-representation of transcript region, 5-mer motif and transcript biotype in IR-responding DM sites using the hyper function from the stats package in R [59].

Clustering of m6A dynamics

2.9.3

Clustering of sites was performed using the divisive hierarchical algorithm ‘DIANA’ via the cluster package in R [60], as implemented in the DEGreport package [61].

Gene ontology

2.9.4

Gene ontology (GO) enrichment analysis was performed using the goana function from the limma package [62]. The number of significant ontology terms was reduced using the minimal_set function of the ontologyIndex library [63]. The proximity of pairs of ontology terms was measured using the Jaccard index between their respective gene sets.

Mapping sites to coordinates and classification

2.9.5

Esembl's REST API, biomaRt [64], and a custom R script were used to identify the chromosomal coordinates of m^6^A-modified sites and classify them by transcript region.

KO genotype interactions with IR m6A response

2.9.6

To explore the role of METTL3, FTO and YTHDF2 on the effect of IR on m^6^A methylation, significant interactions between genotype and time post-IR were identified using a beta-binomial regression model with moderated dispersion on the WT and KO samples including the control, 1 h and 24 h conditions, followed by Benjamini-Hochberg multiple test correction.

Comparison of global shift in methylation across KOs

2.9.7

To test whether the fraction of positive and negative m^6^A responses to IR were different between each of the KOs and the WT we used McNemar's test for paired nominal data in R.

Differential gene expression analysis (DGE)

2.9.8

Direct RNA sequencing reads were aligned to the human transcriptome GRCh38 (release 112) using minimap2. Transcript counts were obtained using salmon. Transcripts were mapped to genes using the R library txdbmaker [65]. Differential gene expression analysis was performed using DESeq2 [57]. Significant differences in gene expression responses to IR between KO genotypes and WT were obtained by filtering genes for significant genotype:treatment interaction terms.

Gene Set Enrichment analysis

2.9.9

To test whether the IR responding sites in the human skin biopsies were enriched for DM sites detected in the WT HT1080 cell line, the fgsea [66] algorithm was used. Beta-binomial regression with moderated dispersions was performed in the same way as for the HT1080 samples, and the resulting effect-sizes and p-values were used jointly to rank the skin biopsy sites and passed to fgsea. The WT HT1080 sites which were DM at 24 h post IR were used as the set being tested for enrichment in the skin biopsy ranked list.

Additional details of the analysis can be found in supplementary material and methods.

Results

3

Characterization of RNA m6A modification dynamics in HT1080 human cell line after IR exposure

3.1

RNA m^6^A modifications were identified in HT1080 cells transcriptome using long-read direct RNA nanopore sequencing technology in a PromethIOn P2 Solo instrument after 1, 2, 10 min and 1 and 24 h after exposure to a 10 Gy dose of X-rays. We obtained an average of 5.5 million reads per sample with a AvrQ score of 19 and AvrN50 of 1.412 kb. Data were then processed through minimap2 to align the reads, followed by nanopolish and m6anet to identify and quantify m^6^A sites. This yielded stoichiometries (the fraction of RNA molecules methylated at a given site) for over 1.5 M unique sites, 223,480 of which were present across all samples. This number was reduced to 9432 after filtering. To identify radiation responsive sites, the change in stoichiometry between treatment conditions was modelled using beta-binomial regression with shrunken dispersions. Most of the DM sites identified were hypermethylated as seen in the volcano plots (Fig. 1A–E). There was a rapid increase in m^6^A methylated sites in response to IR. This was particularly strong at short time points, 1-2 min, with 615 and 528 hypermethylated sites respectively (Fig. 1F, FDR <0.1). A decrease was visible over time as seen at 1 h post-exposure with only 74 hypermethylated sites. 24 h following exposure, the numbers of hypermethylated sites significantly increased again (349 hypermethylated sites) (Fig. 1F). Looking at the global levels of m^6^A methylation with an ELISA assay (Fig. S3), we could not see any major changes after IR over time apart from a slight decrease at 2 min post-exposure. Comparable results were observed in the m6Anet data when averaging the methylation levels across all detected sites (Fig. S4), differing only by a slight increase in m^6^A levels after 1, 2 and 10 min after exposure.Fig. 1. Landscape of differentially methylated (DM) m^6^A sites after X-ray exposure. HT1080 cells were irradiated with 10 Gy X-rays and sampled at 1 min, 2 min, 10 min, 1 h, and 24 h (n = 3 per timepoint). A – E) Volcano plots of per-site differential methylation. Each dot is a candidate m^6^A site; DM sites (FDR <0.1, Benjamini–Hochberg) are red. 10 Gy: change in m^6^A stoichiometry (treated − control). y-axis: −log10(p value). Statistics from moderated beta-binomial regression on modified vs unmodified read counts derived from m6Anet. F) Number of DM sites detected at each timepoint. G – I) Aggregated counts of DM sites across all timepoints, showing (G) the distribution of DM sites across transcript regions (5′UTR, CDS, 3′UTR), (H) a histogram of the number of transcript isoforms on which each DM site occurs, and (I) a histogram of the number of DM sites per transcript. J – K) Over-representation of DM sites by (J) transcript biotype and (K) local 5-mer sequence context. Enrichment is relative to all tested sites as background, using a hypergeometric test with FDR correction; bars show log2 fold-enrichment and significance (∗∗∗∗p < 1e−4, ∗∗∗p < 1e−3, ∗∗p < 1e−2, ∗p < 0.05). L – M) Temporal clustering of DM sites showing distinct methylation dynamics. Stoichiometry per site was z-scored across timepoints and hierarchically clustered; the two largest clusters are shown (Groups A and B). Points show site-level z-scores; boxes indicate median and interquartile range; cyan and green splines trace the cluster means. N) Gene Ontology summary for genes whose transcripts contain ≥1 DM site. Cell numbers are the number of genes with methylated transcripts per term at each timepoint; colour encodes the row-scaled z-score of those counts. Ontology categories: BP (Biological Process), CC (Cellular Component), MF (Molecular Function). Enrichment tested against all expressed genes.Fig. 1

When looking at all the DM sites in our data, the sites are located mainly in the 3′ region (78%) with some presence in the coding region (21.2%) and barely presence in the 5’ region (Fig. 1G). This distribution of the sites within transcript regions is very similar when looking at individual time points for 1, 2 min and 24 h (Fig. S5), and is not significantly different to the distribution of the non-DM sites.

Regarding the uniqueness of the sites identified, 60% of DM sites are found in 2 or more isoforms and only 268 out of 737 DM sites are unique across all transcripts (Fig. 1H). Most transcripts with a DM m^6^A site carry just one modified site (Fig. 1I).

When comparing the distribution of DM sites on transcript types with respect to the background distribution (all detected sites), there was a significant increase in the number of sites on transcripts with retained introns (p = 0.0001) and those targeted for nonsense mediated decay (p = 0.059), while there was a decrease of sites on protein coding transcripts (p = 0.00005) (Fig. 1J and Fig. S6). There was also a significant increase in the number of sites with DRACH motif consensus sequence GGACT (p < 0.0001), and a significant decrease for sequences GAACT (p < 0.0001) and GGACC (p < 0.0001), see Fig. 1K and Fig. S7. While m6Anet is known to have certain biases regarding the surrounding sequence of m^6^A sites, these biases are constant across experimental conditions. Therefore, the inclusion of non-irradiated controls allows us to control for these biases, which in turn allows us to isolate the effects of IR on the types of m^6^A sites affected most strongly.

Sites found to be DM in at least one time point were clustered into subcategories based on their m^6^A methylation levels over time (Fig. S8). The two largest clusters (representing 35% of DM sites) were characterised by a fluctuating pattern over time (Fig. 1L) and a gradual decrease in methylation over time without returning to the background levels of the non-exposed samples (Fig. 1M).

We looked at the gene ontology analysis of all genes with m^6^A significantly modified sites after IR exposure; it revealed that the most affected pathways are involved in bioenergetics, cell signalling, migration and apoptosis (Fig. 1N and Fig. S9).

IR induced differently methylated sites with sustained hypermethylation over time

3.2

Comparing the methylated sites across the different time points, we observed that the short time points after exposure, 1 and 2 min, shared the most DM sites, with a total of 280 sites (FDR<0.05). When comparing all the time points except 1 h post exposure, the rest of the time points only share 35 sites. When combining all the time points, only 4 common m^6^A sites were identified which remained hypermethylated over time as shown on the Venn diagram (Fig. 2A). Two of these 4 sites were found in two different isoforms of the same gene, UQCR10; variant 1 UQCR10-201 (ENST00000330029.6 position 374) and variant 2 UQCR10-202 (ENST00000401406.3 position 428) (Fig. 2B–D). These sites are in the non-coding part of the second exon in both variants and followed the methylation dynamic patterns of group A (Fig. 1L). The results obtained with the sequencing for those two methylated sites located in two UQCR10 variants were validated with a RT-QPCR based approach [48,49] (Fig. 2E). Both methods showed a significant hypermethylation at 1 h and 24 h post irradiation time points.Fig. 2. Selection and validation of persistent DM m^6^A sites following X-ray exposure. A) Venn diagram showing the significant DM sites at each timepoint (FDR <0.05). Four sites persistently significant across all timepoints, two of which are shown in B) Genomic location of m^6^A site persistently hypermethylated across all timepoints, found in two distinct transcript isoforms of the gene UQCR10. Red boxes: exons; arrows: transcript orientation; vertical blue line: shared genomic coordinate of the DM site (hg38, chr22). C – D) Methylation dynamics for the site in each isoform, both of which belong to cluster Group A (Fig. 1L). Boxplots show per-replicate m^6^A stoichiometry (n = 3) across timepoints; the line traces the mean trajectory. E) qPCR validation of the DM sites in the UQCR10 transcripts. Asterisks denote significance vs control (paired t-test, p ≤ 0.05). Asterisks denote significance vs control (paired t-test, p ≤ 0.05).Fig. 2

RNA m6A profile in METTL3, YTHDF2 and FTO KO HT1080 cell lines after IR exposure

3.3

To better characterize the mechanisms involved in the effect of radiation on the increase in the deposition of methylated sites, we generated 3 HT1080 cell lines deficient in the expression of specific enzymes involved in m^6^A RNA methylation metabolism using CRISPR-CAS technology. To cover the 3 steps of the m^6^A cycle, we selected a writer, METTL3, an eraser, FTO and a reader YTHDF2, (Table S1). We confirmed the absence of these specific enzymes by Western blot (Fig. S1) and compared the proliferation capacity at 2 Gy and 10 Gy doses (Fig. S10). A joint non-linear mixed-effects model was used to compare WT vs KO within each dose condition. The results show that the FTO KO cell line grows significantly slower than the WT (p < 0.001) across all conditions, including in the absence of irradiation (0 Gy control). YTHDF2 KO cells grow faster compared to WT cells under no exposure conditions (p < 0.001), while dividing slightly slower following a 2 Gy dose (p < 0.1). We note however that each condition was modelled separately, so we cannot test whether the differences in growth rate observed at 0 Gy and 2 Gy are significant. METTL3 KO cells grew slightly faster than WT cells at 0 Gy and 2 Gy (p < 0.1) and showed no difference in growth rate compared to WT at 10 Gy.

The overall abundance of m^6^A methylated sites in the KO cell lines was measured using an Elisa assay (Fig. S11). The results did not show any differences for METTL3 at 1 h and 24 h after exposure, but for YTHDF2 and FTO, a decrease was observed at 24 h. However, when aggregating the modified reads across all sites discovered by m6Anet to get an equivalent score from the sequencing data, no differences were found between the exposed and non-exposed samples (Fig. S12).

We also performed direct RNA nanopore sequencing analysis in the KOs cell lines exposed to 10 Gy at 1 h and 24 h after exposure. We obtained an average of 6.8 million reads per samples with a AvrQ score of 18.07 and AvrN50 of 1.4 kb. We found that, regardless of the enzyme involved in m^6^A metabolism being knocked out, the lack of one of them induced a significant transcriptome-wide decrease in radiation induced hypermethylation levels (Fig. 3A–H), with an increase in hypomethylation for all of them, to a greater extent for YTHDF2 (Fig. S13). A clear decrease of methylation levels for all DM sites identified in the WT after exposure was found in all three KO clones (Fig. 3I), with a more pronounced effect in YTHDF2 followed by METTL3.Fig. 3. Overview of m^6^A response to IR in the HT1080 KO cell lines exposed to 10 Gy and sampled at 1 h and 24 h (n = 3 per group). WT shown for comparison purposes. A) to D) Volcano plots at 1 h for WT, and the 3 KOs FTO, METTL3, and YTHDF2 respectively. Each dot is a candidate site; red = DM (FDR <0.1). x-axis: change in stoichiometry (irradiated − matched non-irradiated control within genotype). y-axis: −log10(p value). Statistics from moderated beta-binomial regression on modified vs unmodified read counts derived from m6Anet. E – H) The same sampled at 24 h. I) Heatmap showing the degree of hypermethylation in each of the different genotypes at 1 h and 24 h post-exposure. Colour intensity is the scaled (row-wise z-score) change in stoichiometry with respect to the genotype-matched control condition. Only sites found to have a significantly different methylation response to IR between each of the KOs and the WT shown (significant genotype × treatment interaction, FDR <0.1). J) Boxplots showing the per-replicate change in m^6^A stoichiometry of the site found on transcript UQCR10-201, shown for each genotype and for each timepoint. Control conditions at baseline (dashed horizontal line) added for comparison purposes.Fig. 3

The radiation-induced hypermethylation reported in the WT cell line for two UQCR10 transcripts were indeed unresponsive to radiation in the three KO cell lines at 1 h and 24 h (Fig. 3J and Fig. S14).

Comparing transcriptomic analysis with m6A modifications in WT and METTL3, YTHDF2 and FTO KO HT1080 cells

3.4

Transcriptomic analysis of the WT cell line showed an increase of differentially expressed genes (DEGs) over time after a 10 Gy irradiation dose (Fig. 4A). The transcriptomic signature at 24 h included genes known as radiation responsive such as CDKN1, FDXR and MYC among others (Fig. S15). However, the transcriptional levels of METTL3, YTHDF2 and FTO were not affected by IR exposure in this cell line (Fig. S16). Comparisons between gene ontology on the DM sites (Fig. 1N) and on the gene expression at 24 h (Fig. 4C) presented a very modest overlap, bringing forward several pathways: apoptosis, locomotion, wound healing and animal organ development.Fig. 4. Integration of transcriptomic and epitranscriptomic responses to 10 Gy X-ray exposure in HT1080 WT and KO cell lines. A) Bar plot showing the number of differentially up- (red) and down-regulated (blue) genes (DEGs) at each timepoint (1 min – 24 h) for the WT HT1080 cell line, (FDR <0.1). B) Similar to A but for the KO cell line timepoints (1 h and 24 h), WT timepoints also shown for comparison purposes. C) Relationship between Gene Ontology terms enriched in differentially methylated transcripts (p < 1e-3; columns) at any timepoint post-exposure, and differentially expressed transcripts (p < 1e-3; rows) at the 24 h timepoint. Heatmap colour shows the Jaccard index between gene sets of the corresponding Gene Ontology terms. Only terms with at least one Jaccard index greater than 0.25 are shown. D – E) Genes with a significantly different transcriptional response to IR between WT and KO cell lines at 1 h and 24 h respectively. Heatmap colour intensities (and numbers for the 1 h) correspond to the expression log2 fold-change. Genes are hierarchically clustered in 3 and 6 main groups, for 1 h and 24 h respectively. The black lines on the left of each gene indicates which KO cell line has a significantly different response in expression with respect to the WT. The red lines on the left of each gene indicate whether the gene was hypermethylated in the WT at any timepoint after the treatment. Gene Ontology analysis of each of the clusters shown in Table S3.Fig. 4

When comparing the effect on the transcriptomic response to radiation of METTL3, YTHDF2 and FTO KO clones, similar numbers of DEGs were found compared to the WT both after 1 h and 24 h post exposure, apart from METTL3 at 1 h which presented a decrease in DEGs (Fig. 4B). Interestingly, each KO presented a specific expression profile, differing from the WT at 1 h and 24 h (Fig. 4D–E and Fig. S15). However, looking into individual genes, some were commonly identified in the KOs and in the WT cells in response to radiation. This is the case for BTG2 1 h post exposure for YTHDF2 and FTO, and at 24 h for SUGCT gene for METTL3 and YTHDF2, CDKN1 and FDXR genes for YTHDF2 and FTO (Fig. S15). Differential gene expression responses to the treatment between the WT and KO cell lines after 1 h identify a signature of 18 genes. Lack of YTHDF2 presented the strongest effect on these genes’ expression compared to the WT (Fig. 4D).

Pathway analysis of the differentially expressed genes at 24 h after exposure revealed that the lack of FTO affects more specifically pathways associated with inflammation, cell cycle, cell proliferation and DNA replication and cell communication (Fig. 4E and Table S3). For METTL3, it affects DNA repair and replication, cell proliferation, differentiation and communication. YTHDF2 KO affects more specifically cell survival and function, cell division and cell differentiation and cell communication (Fig. 4E). Looking at the distribution of genes with transcripts carrying ≥1 DM m^6^A site in WT (DM transcripts), the majority of them were found in 3 of the 6 clusters (Fig. 4E, red ticks). In cluster 2, we found 14% of genes within the cluster were associated with DM transcripts (Fig. S17). This cluster is characterised by genes with a differential transcriptional response of FTO to IR with respect to the WT (Fig. 4E, black ticks to the left), suggesting that FTO and m^6^A deposition may be involved in inflammation related pathways (Table S3). For cluster 4, we found 14% of genes within the cluster were associated with DM transcripts (Fig. S17). This cluster is characterised by the differential transcriptional response of METTL3 to IR with respect to the WT, showing a clear decrease in the overexpression of genes linked to the PCNA-p21 complex, pathways involved in DNA replication and repair, particularly in response to DNA damage (Table S3). In cluster 3, with 8% of genes within the cluster associated with DM transcripts (Fig. S17), all KOs had a differential transcriptional response with respect to the WT for the majority of the genes within the cluster. These were enriched for genes involved in the haemostasis and response to injury pathways (Table S3). Namely, the downregulation of these pathways seen in WT following IR exposure was largely attenuated across all KO cell lines, and a subset of genes was overexpressed in the FTO KO.

RNA m6A profile in human skin biopsies exposed to IR

3.5

Finally, we wanted to validate our findings in non-cancerous cells and carried out experiments with human tissue. We therefore used human skin biopsies to test whether any of the sites found to be DM in the HT1080 samples were also hypermethylated in skin following IR exposure. Human skin biopsies from 5 donors were exposed to 2 Gy X-rays and incubated for 24 h at 37 °C, and paired-donor non-irradiated samples were included. The sequencing analysis performed with a new version of the direct RNA nanopore sequencing chemistry and flow cells presented an increase of reads per sample to an average of 9.7 million reads, with AvrQscore of 13.60 and AvrN50 of 1.396. The resulting sequencing data were analysed in the same way as for HT1080 cells, and the identified sites were ranked by the strength of their response to IR (Fig. 5A and B). Importantly, using a site-level GSEA implementation, we found that sites previously found to have a hypermethylated response to IR in HT1080 were significantly enriched in the most strongly responding sites in the skin (Fig. 5C). Furthermore, the hypermethylated site on UQCR10-202 (ENST00000401406.3_428 variant) was among the highest ranking (Leading Edge) sites in skin biopsies (Fig. 5D).Fig. 5. Effect of 2 Gy X-ray exposure on m^6^A levels in human skin biopsies after 24 h post-exposure (n = 5, each paired with a non-exposed matching-donor sample). A) Volcano plot showing the site-level methylation response at 24 h post-exposure. x-axis is the change in m^6^A methylation (logit scale), y-axis is the significance (p value, -log10 scale). Statistics obtained using a paired moderated beta-binomial regression model on modified vs unmodified read counts derived from m6Anet. Colour indicates the m^6^A site ranking (lighter = higher rank) which was used for the site-level GSEA analysis. B) Site-level Gene Set Enrichment plot, showing that the m^6^A sites in human skin with increased methylation at 24 h post-exposure are significantly enriched for sites found to be DM in HT1080 following 10 Gy exposure (NES = 1.62, p value = 0.008). x-axis shows the rank of each site in the human skin, y-axis is the running sum statistic, dashed red line denotes the GSEA score; blue ticks mark the positions of the HT1080 DM m^6^A sites in the skin biopsy ranking; the green area denotes the leading-edge (LE) subset (DM m^6^A sites contributing most to the enrichment). C) Heatmap showing the scaled m^6^A stoichiometry (row-wise z-score) before and after X-ray exposure across the 5 donors, for the genes with m^6^A sites found in the LE subset of the site-level GSEA. D) Boxplots showing the m^6^A stoichiometry of the site found on the UQCR10-202 transcript in the 5 donor paired skin biopsies, before and after X-ray exposure.Fig. 5

Discussion

4

The deposition of m^6^A modifications in specific adenines have a posttranscriptional regulatory role which regulates multiple aspects of the RNA metabolism. These chemical modifications have been associated with the integration of responses to external stimuli and environmental changes such as genotoxic factors including IR [32,42,44]. Given our expertise in IR transcriptional response [[67], [68], [69]], we aimed to characterize this epigenetic mark in response to IR exposure using direct RNA sequencing. In the present study, we have examined the dynamic response of m^6^A modifications to IR in a human cell line, HT1080, and further validated some of our findings in human skin biopsies. To our knowledge, we reported for the first time that the response is characterised by a remarkably rapid (∼1 min) site specific increase in m^6^A methylation. More specifically, we have demonstrated that IR regulates RNA m^6^A methylation metabolism by increasing the deposition of methyl groups in N6 position of adenosine mainly on GGACT motifs located in the 3′UTR region and at a lesser extend in the coding region of certain transcripts. IR exposure mainly induces hypermethylation, and after accounting for variable transcript abundance, it is biased toward intron-retaining transcripts (and, to a lesser extent, nonsense-mediated decay targets), while protein-coding transcripts are proportionally under-represented. Using KO cell lines, we showed that these regulatory effects were drastically affected by the depletion of specific enzymes involved in m^6^A methylation demonstrating an interplay between IR and RNA methylation.

High-throughput long-read direct RNA nanopore sequencing technology was used to generate the data which provided information not only on RNA methylation at single-base resolution, but also in gene expression at isoform level without amplification biases [70]. This technology has been used for detection of m^6^A levels in multiple scenarios such as blood from astronauts [71], MOLM13 leukemia cell line [72] and HepG2 [73] cells. There are several tools available for estimating m^6^A stoichiometries from direct RNA sequencing data and m6Net is currently the best option for human data [16,74]. However, currently there is no established method for comparing the output of m6Anet between different experimental conditions, in other words for performing differential methylation analysis. While others have used simple logistic regression for this purpose [71,75], we showed that this fails to capture the technical and biological variability present in the data, therefore underestimating the false discovery rate of DM sites. For that reason, we employed an alternative framework based on beta-binomial regression [56] and robust estimation of dispersion using an empirical Bayes approach [76]. This allowed us to better estimate the true variability in the data and to detect DM sites with greater precision.

RNA DM analysis revealed that IR triggers an immediate activation of the m^6^A methylation machinery, reflected by an increase in DM sites in HT1080 cells within the first minute of exposure. Similar rapid responses have been previously reported in U2OS exposed to UV light, linked to DNA repair [33]. Since genotoxic factors introduce single and double DNA strand breaks which activate the DNA damage response in a matter of seconds, it is not surprising that IR and UV regulation of m^6^A methylation could both be associated with the DNA damage response [33,36]. Since the induction of transcription by IR exposure does not occur within the first minute after exposure [68,77,78], we hypothesize that the deposition of m6A sites is associated with pre-existing mRNA pools. Post-transcriptional m6A modifications of pre-existing mRNA transcripts in response to oxidative stress has already been observed within minutes after cells were treated with hydrogen peroxide [79]. This study supports our observations of an early increase in m6A deposition, consistent with the known rise in reactive oxygen species following IR exposure [80]. After the initial peak in methylation within the first minute post-exposure, the number of hypermethylated m^6^A sites decreased drastically over the first hour post-exposure. The addition of the methyl group is reversible through the action of methyltransferases and demethylases, and these chemical modifications have been previously described as transient [33]. However, at 24 h post-exposure, we observed another surge in methylation. This late increase is probably not related to the initial DNA damage response and activation of the different DNA repair pathways. We hypothesize that this may be associated with the known disruptions of the cell cycle [81] and cell proliferation [82] at 24 h after IR exposure.

Bulk m^6^A levels assessed with an ELISA assay were somewhat different to the sequencing data presented; differences are probably due to their fundamentally different approach to detect and quantify m^6^A methylation. An ELISA assay is a relatively simple but rapid method to evaluate the overall level of m^6^A methylation, but it provides a global measurement of the total m^6^A levels [83]. On the other hand, direct RNA nanopore sequencing can identify sites at a single-nucleotide resolution. ELISA's limitations such as issues with antibody specificity, cross-reactivity, and potential contamination of rRNA, could be the source of a certain degree of bias in the data [83].

Looking into individual DM sites, multiple dynamic patterns of methylation in response to IR over time were observed. m^6^A modifications are not static, they are dynamic and adaptive to modulate gene expression depending on the context and linked to different functional outcomes [84]. m^6^A deposition can occur co-transcriptionally; therefore, fluctuations in m^6^A levels over time are likely associated with newly synthesized transcripts rather than pre-existing mRNAs [85]. However, m^6^A modifications can also be installed on pre-existing transcripts [86]. Accordingly, we hypothesize that the observed dynamics in m^6^A reflect a cellular response to the stimulus [68,69,78], although the mechanisms by which mRNAs are regulated at early versus late time points after exposure and whether the DM sites have been deposited in nascent or pre-existing mRNAs remain unclear; in our case, these modifications might be associated to different metabolic processes cells undergo after radiation exposure, such as cell cycle arrest, DNA repair or cell death. Indeed, gene ontology analysis, performed on the transcripts carrying DM sites, identified specific pathways involved in DNA repair and apoptosis such as protein disulfide isomerase, mitochondrial functions or cell-matrix adhesion, connecting the dynamic responses observed with potential biological responses.

In this study, most of the transcripts carry only one DM site but a small proportion carry multiple sites. The location of the sites and the number of sites per transcript can affect the post-transcriptional regulation of the mRNA. For instance, mRNA transcripts containing multiple sites act as scaffolds, gathering different YTHDF proteins leading them to different cell compartments and affecting their fate [87]. Around 78 % of the DM sites after radiation exposure were found in 3′ UTR while 21.2 % were located inside coding regions. Previous studies have pointed out that sites situated near stop codons and 3’ UTRs affect mRNA stability [88,89], turnover and translation efficiency [90], but have a different role in coding regions promoting mRNA degradation [91]. Looking at the transcript biotypes enriched in DM sites, IR promoted deposition on alternative transcripts, specifically non-sense mediated decay and intron retention variants, with a significant decrease in protein coding variants. It is known that IR induces alternative splicing [92,93], probably as an adaptative response to cell stress which might be mediated through m^6^A deposition since this chemical modification is known to contribute on regulating alternative splicing [21].

Exploring the DMs across all the time points, we identified 4 sites which remained hypermethylated over time after IR exposure. Two of them are in variants of the gene UQCR10. Both sites are located in the last exon, which is a predominant region for m^6^A methylation in mRNAs regardless of the location in the coding or non-coding part of the sequence [89]. UQCR10 is a subunit of mitochondrial complex III, part of the mitochondrial electron transport chain and its overexpression has been associated with a decrease in apoptosis [94]. Considering that the location of the methylated site could be associated with degradation [89], this could possibly have a role in regulating apoptosis, a pathway well known to be triggered by doses of IR high enough to trigger cell death. It would be interesting to assess whether hypermethylation of this site can be detected at sublethal doses where cell death via apoptosis is minimal. Moreover, IR is known to disrupt the electron transport chain equilibrium through the increase in reactive oxygen species (ROS) [95,96]. ROS induces mRNA m^6^A methylation through the upregulation of mRNA m^6^A writers at protein and transcript level and by inhibition of ALKBH5 demethylase activity [79]. Therefore, this evidence suggests that the effect of IR on mRNA m^6^A methylation could be mediated through IR-induced ROS production in the cells. However, this regulatory effect might be independent of the transcriptional regulation of METTL3, YTHDF2 or FTO, since the expression of these enzymes is not affected by IR in HT1080 cells.

While the sustained methylation level of these two sites on UQCR10 variants over the first 24 h period after exposure is mechanistically intriguing; it does also evidence the potential use of these sites as biomarkers of radiation exposure. The use of m^6^A methylation levels for triage and biodosimetry have been recently explored [42,43]. m^6^A methylation levels in the nuclear receptor coactivator 4 (NCOA4) gene exhibit and increase in a dose- and time-dependent manner after IR exposure in mouse blood [42]. However, this methylated transcript was not confirmed in our study, suggesting that specific m^6^A site deposition is cell line, tissue and/or species dependent [97]. Given that m^6^A levels were also elevated in UQCR10 in human skin biopsies exposed to a much lower dose (2 Gy) as in HT1080 cells, and that NCOA4 was also detected in human endothelial and blood cells [42], the discrepancy could be due to different experimental conditions, in particular the different time scale sampling used in both studies.

Transcriptomic analysis of the HT1080 cell line showed an increase in differential gene expression over time. At 24 h post-exposure, several known radiation responsive genes were identified such as CDKN1A and FDXR as expected [98]. Next, we investigated the potential regulatory role of m^6^A on the gene expression response to IR. m^6^A methylation does not regulate every single transcript in a cell, instead, it is deposited in a specific subset of transcripts [19]. Therefore, we decided to correlate genes with DM transcripts and differentially expressed genes following IR. To this end, we compared the results of Gene Ontology analysis of DEGs at 24 h after exposure with the Gene Ontology of genes with transcripts carrying at least one DM site. While we are not demonstrating a causal link, in both cases we found differential regulation of genes associated with apoptosis, cell motility, response to injury and cell fate in organ development.

Furthermore, as we observed that IR-induced cell stress leads to m^6^A methylome regulation, we decided to explore the potential regulatory mechanisms involved by knocking out 3 major m^6^A methylation enzymes involved in the 3 phases of m^6^A metabolism: METTL3 (“writer”), YTHDF2 (“reader”) and FTO (“eraser”). The initial observation was a decrease in hypermethylation following IR exposure in all the 3 KOs, suggesting that cells require an intact m^6^A enzymatic machinery to respond effectively at the epitranscritomic level. We observed that this was also occurring at individual m^6^A sites such as UQCR10 sites in all KO cell lines, at least at 1 h and 24 h after exposure. Surprisingly, depleting a writer, reader or an eraser had a similar effect on the IR response, irrespectively of its role. A decrease in methylation in cells lacking METTL3 has been previously reported [33]. However, the lack of FTO seems to increase m^6^A after UV exposure in U2OS cells [33]. Therefore, the unexpected decrease in m^6^A methylation observed after FTO depletion might be related to indirect downregulation of the m^6^A methyltransferase complex. Similar, counterintuitive results were observed after YTHDF2 depletion. Since this enzyme is a reader whose function relies on binding to m^6^A sites and targeting them for degradation, its loss would be expected to increase the abundance of certain mRNA transcripts [99]. However, in our case, the loss of YTHDF2 appears to affect the RNA turnover, possibly by promoting alternative decay mechanisms, which could lead to a decrease in the abundance of m^6^A-marked RNAs. There is evidence suggesting that loss of this reader can destabilise m^6^A-modified transcripts, resulting in their degradation [100].

The transcriptional changes induced by IR in these three KO cell lines lead to different signatures linked to different regulatory pathways. Although the overall m^6^A methylation response to radiation in these KOs was similar, the consequences at the gene expression level and in the pathways involved did differ. Investigating the transcriptional signatures at 24 h after exposure in the KOs, we observed that in the case of METTL3 KO, well characterized radiation responsive genes such as FDXR and CDKN1 were not significantly up-regulated. This response may reflect the established regulatory role of METTL3 in controlling p53 target mRNAs [101] after DNA damage, leading to the expression of target genes through different downstream effects. Therefore, in the absence of METTL3, the IR-induced activation of p53 downstream targets might appear to be disrupted. Pathway analysis indicated the potential regulatory effect of METTL3 on DNA repair, proliferation and differentiation which is strongly supported by other studies. For example, METTL3 is known to regulate homologous recombination repair [45] and loss of METTL3 causes accumulation of unrepaired double strand breaks [36]. For example, in liver cancer, inhibition of METTL3 leads to the suppression of tumour cell growth by inducing apoptosis [102] and in neural progenitors cells, an intact METTL3 is required for proliferation [103].

In the case of FTO, pathway analysis indicated that loss of FTO could have a potential role in modulating inflammation. Previous studies in FTO-deficient models have demonstrated its role in inflammation with either a pro- or anti-inflammatory effect depending on the stimulus and cell type [104,105].

Under YTHDF2 deficiency conditions in our cell model, cell survival, division and differentiation seem to be affected. YTHDF2 regulatory effect on cell survival has been observed in neuroblastoma cells after its silencing, leading them to apoptosis and inhibition of proliferation [106]. Depletion of YTHDF2 have also shown to affect cell division for example in HeLa cells [107] and lung squamous cell carcinoma cells [108], and cell differentiation in neural stem/progenitor cells [109].

Although this study revealed changes in the distribution and abundance of m^6^A marks in response to radiation in two different models, there are some limitations related to the functional impact of the identified modifications in specific transcripts. Future research will allow to understand how these modifications might affect RNA stability, translation, splicing and degradation. This may be achieved by applying methods such as RNA immunoprecipitation followed by sequencing, polysome profiling, ribosome profiling, metabolically labelling RNA or computational tools like FLAIR2 [110].

In conclusion, we reported for the first time the early dynamics of m^6^A following IR exposure in a human cell line and further confirmed in human skin using direct RNA long read sequencing. Our findings have demonstrated that IR regulates m^6^A metabolism, mainly inducing hypermethylation in a dynamic way. We have also revealed a regulatory axis between IR and m^6^A methylation machinery, demonstrating that functional writers, reader and erasers (respectively METTL3, YTHDF2 and FTO) are necessary to integrate the response of IR into regulation of m^6^A depositions. Finally, taking advantage of the single-nucleotide resolution offered by the sequencing technology used, we identify a specific site with sustained methylation over time making it a good candidate as a biomarker of radiation exposure. All our results demonstrates that IR interacts with the m^6^A regulatory layer to rapidly impact on the modulation of gene expression. The results in this study pave the way for further m^6^A research following radiation exposure; the speed of response is fascinating and requires further exploration. In radiation oncology, it would be of particular interest to explore the potential roles in both tumour response and normal tissue toxicity. In the context of individual sensitivity, m^6^A may become a prognostic biomarker.

CRediT authorship contribution statement

Lourdes Cruz-Garcia: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Writing - Original Draft, Writing - Review & Editing, Supervision, Project administration, Resources, Visualization.

Philip Davies: Conceptualization, Methodology, Validation, Formal analysis, software, Investigation, Writing - Original Draft, Writing - Review & Editing, Data Curation, Visualization.

Veronika Goriacha: Formal analysis, Writing - Original Draft.

Mustafa Najim: Investigation, Writing - Original Draft.

Stanislav Polozov: Methodology, software, Formal analysis, Data Curation.

Maria Polozova: Investigation, Writing - Original Draft.

Christophe Badie: Conceptualization, Original Draft, Writing - Review & Editing, Supervision, Project administration, Resources.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.