The Alteration of M6A-Tagged Transcript Profiles in the Retina of Rats After Traumatic Optic Neuropathy

Messager RNA (mRNA) can be modified in a variety of ways, among which the modification of N6-methyladenosine (m6A) is one of the most common ones. Recent studies have found that the m6A modification in mRNA could functionally regulate the splicing, localization, translation, and stability of mRNA, which might be closely related to multiple diseases. However, the roles of m6A modification in traumatic optic neuropathy (TON) are unknown. Herein, we detected the expression of m6A-related genes via quantitative real-time PCR (qRT-PCR) and performed methylated RNA immunoprecipitation sequencing (MeRIP-seq) as well as RNA-sequencing to analyze the alteration profiles of m6A modification after TON. The results showed that the expression of m6A-related genes (METTL3, WTAP, FTO, and ALKBH5) were all upregulated after TON. In all, 2,810 m6A peaks were differentially upregulated and 689 m6A peaks were downregulated. In addition, the hypermethylated and hypomethylated profiles of mRNA transcripts were also identified. To sum up, our study revealed the differentially expressed m6A modification in the early stage of TON, which may provide novel insights into the mechanism and treatment of TON.


INTRODUCTION
Traumatic optic neuropathy (TON) refers to a common complication of traumatic brain injury (TBI), and the incidence of TON ranges from 1.5 to 4% (Guy et al., 2014). The sufferers of TON mainly showed vision impairment and even ablepsia, bringing huge burdens to their family and the whole of society (Sosin et al., 2016;Bastakis et al., 2019). Although numerous clinical trials and fundamental research has been conducted, therapeutic strategies for TON are still limited, which may be owing to the insufficient perception of the pathophysiological mechanisms in the progression of TON (Chaon and Lee, 2015;Sosin et al., 2016). Nowadays, corticosteroids and optic nerve decompressive surgery are major treatments for TON, though neither of them could reach definitive curative effects (Emanuelli et al., 2015;Oh et al., 2018). Therefore, it is urgent to explore more effective ways to improve the vision of TON patients. More recently, noncoding RNAs (ncRNAs) and specific RNA modification, including circular RNAs (circRNAs), long-noncoding RNAs (lncRNAs), microRNAs (miRNAs), and N6-methyladenosine (m6A) modifications have been found to be expressed differentially after Frontiers in Genetics | www.frontiersin.org 2 February 2021 | Volume 12 | Article 628841 TBI, which provided new directions for elucidating the physiological process of TBI (Wang et al., 2019;Qu et al., 2020). However, there is no study on the relationship between m6A modification and TON. The roles of m6A-tagged transcripts in TON need to be elucidated. Generally, there are two covalent modifications at the two termini of mRNA:7-methylguanosine cap structure and polyadenylate tail (Ramanathan et al., 2016;Stewart, 2019). Besides, a large number of nucleotide modifications have been found to be abundant in eukaryotic mRNAs, such as m6A, 5-methylcytosine (m5C), and N1-methyladenosine (m1A; Maity and Das, 2016). m6A modification is the most common one in eukaryotic mRNAs, accounting for approximately 80% of the total methylation modifications (Li et al., 2019). The development of methylated RNA immunoprecipitation sequencing (MeRIP-seq) and other technologies have increased our understanding of the role of m6A modifications in mRNAs. Regulated by methyltransferases (METTL3, METTL14, WTAP, etc.) and demethylases (FTO,ALKBH5,etc.), m6A modification is capable of participating in various fundamental biological mechanisms, for example, the regulation of gene expression, self-renewal of neural stem cells, the formation of metabolic diseases, the maintenance of biological circadian rhythms, and the development of cancers (Fischer et al., 2009;Fustin et al., 2013;Zhou et al., 2015Zhou et al., , 2018Maity and Das, 2016;Hsu et al., 2017).
In view of the inseparable relationship between m6A modification and protein translation, it is advisable to speculate that the dysregulation of m6A modification may be correlated with the pathophysiological processes following TON. The combination of quantitative real time-PCR (qRT-PCR), MeRIP-seq, RNA-sequencing and bioinformatic analysis were used to figure out the alteration profiles of m6A modification in the rat retina after TON. Furthermore, this study also revealed some potential roles of m6A-modified transcripts in TON.

Animals
The protocols for animal studies were approved by the Animal Ethics Committee of the Naval Medical University (Shanghai, China), Ethic Certificate No: 81371382. All experiments were performed under the guidance of the National Institute of Health Guide for the Care and Use of Laboratory Animals. Male Sprague-Dawley (SD) rats (220-250 g, 7-8 weeks of age) were housed for at least 7 days in a temperature-regulated (22-25°C) and humidity-controlled (50% relative humidity) animal faculty with a 12-h light/dark cycle. During the experiments, all rats had free access to water and food, except that food was withheld overnight before surgery.

Establishment of the Experimental Traumatic Optic Neuropathy Model
In this study, the TON model was established by clamping the optic nerve of rats, as previously described (Jiang et al., 2013;Ibrahim et al., 2018;Oku et al., 2019). Briefly, adult male SD rats were intraperitoneally anesthetized with pentobarbital sodium (1%, 35 mg/kg body mass) and placed on an electric heating blanket to ensure the stabilization of body temperature at 37°C (Zhou et al., 2019). After anesthesia, microsurgical instruments were used to cut the upper eyelid skin vertically under a microscope. The optic nerve was then exposed by cutting the bulbar conjunctiva and conducting the blunt separation of retrobulbar adipose tissue. To cause the injury, the optic nerve was directly clamped at 2 mm behind the eyeball for 10 s by forceps with a 40 g constant holding force. After injury, the surgical incision was sutured, chlortetracycline eye ointment was used for anti-infection. In the sham-operated group, rats received identical surgical procedure, but without clamping applied to the optic nerve. A total of 40 rats were assigned to the TON model group and 36 to the sham group.

The Separation and Purification of Retina
The extraction of retinal tissue was in accordance with the protocol provided by Nie et al. (2018). In brief, 12 h after operation, rats were decapitated under sterile conditions. Under the microscope, the eyeball was cut off along the limbus of the cornea, and the cornea, lens, and vitreous body were separated to obtain the retinal tissues. The separated retinal tissue was then detached by 0.125% trypsin and maintained at 37°C for 20 min. The supernatant was removed by centrifuging the digested retinal tissue at the rate of 200×g, followed by adding 0.25% trypsin inhibitor (T9003; Sigma-Aldrich Chemical Company, MO, United States) for neutralization and re-centrifugation. After the supernatant was discarded, the retinas were rinsed with Kreb's solution comprising Mg 2+ and 1% bovine serum albumin (BSA) three times. After the culture medium was added, tubularis was used to triturate and disperse the cells. The single-cell suspension was then transferred to culture dishes covered by 0.1 mg/ml of 192 poly-L-ornithine (P4538; Sigma-Aldrich Chemical Company, MO, United States) and 1 μg/ml laminin (l6274; Sigma-Aldrich Chemical Company, MO, United States). Then, Basal medium eagle (BME) was added comprised of 25 μmol/L glutamine, 10% fetal bovine serum (FBS), and 0.1 mg/ml gentamicin (SBJ-ME1268, SenBeiJia Biological Technology Co., Ltd., Jiangsu, China) and incubated at 37°C with 5% CO 2 . At last, The suspension was placed in BME medium comprised of rat Thy-1.1 antibody (ab44898; 1:1000; Abcam, Cambridge, United Kingdom) and goat anti-rat immunoglobulin G (IgG) antibody (ab150077, 1:1000; Abcam, Cambridge, United Kingdom) with incubation for 30 min at 37°C.

Quantitative Real-Time PCR Assay
The expression levels of m6A-related genes FTO, ALKBH5, METTL3 and WTAP in the retina of injured and sham-operated rats were analyzed by using qRT-PCR. Total RNAs were extracted from the purified retina using Trizol reagent (Invitrogen, CA, United States) under the guidance of the manufacturer's protocol. The concentration and purity of total RNA were measured by NanoDrop ND-1000 (Thermo Fisher Scientific, MA, United States). Reverse transcription was conducted using PrimeScript RT Master Mix (Perfect Real Time; Takara, Shiga, Japan) to synthesize complement DNA. The amplification of Frontiers in Genetics | www.frontiersin.org 3 February 2021 | Volume 12 | Article 628841 qRT-PCR was performed through pre-denaturation, 40 PCR cycles and establishing the melting curve. qRT-PCR was carried out in ProFlex 96-well PCR System (Thermo Fisher Scientific, MA, United States) with SYBR Green master mix (Applied Biosystems, CA, United States). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal parameter to normalize the data. Comparative threshold method (2 −ΔΔCT method) was used to acquire relative expression levels, which were calculated according to the formula (Livak and Schmittgen, 2001): Where Ct, replicate-averaged threshold cycle. ∆Ct, the Ct of the target gene subtracted from the Ct of the reference gene. Each qPCR reaction was repeated three times.
The full base sequences of METTL3, ALKBH5, WTAP, FTO, and GAPDH were downloaded from NCBI in line with GenBank sequence number, and the primers were designed by Primer premier V6.0 software (Premier Biosoft International, United States). The primers used in this procedure are presented in Table 1.

RNA-Seq
Twelve hours after operation, rats were anesthetized with pentobarbital sodium and were transcardially perfused with 20 ml of 4°C PBS. The retinal tissues were rapidly dissected and stored in liquid nitrogen. According to the manufacturer's protocol, total RNAs from the retina of each rat were isolated using TRIzol reagent (Invitrogen). Three replicates were conducted for the sham-operated and TON model group. To attain 60 μg of total RNA, retinal samples from four rats were integrated into one centrifuge Tube. The concentration and quality of total RNAs were assessed using NanoDrop ND-1000 (Thermo Fisher Scientific). The integrity of RNA was evaluated by denaturing agarose gel electrophoresis. The extraction of mRNA was conducted using NEBNext rRNA Depletion Kit (New England Biolabs, Hertfordshire, United Kingdom) following

Name
Sequence the manufacturer's protocol. Bioanalyzer 2,100 system (Agilent Technologies, CA, United States) was used for quality control and quantification of the RNA libraries, and Illumina Hiseq4000 platforms (Illumina, CA, United States) was used for doubleended sequencing of RNA. Paired-end reads were harvested from Illumina HiSeq 4,000 sequencer and were quality controlled by Q30. After 3' adaptortrimming and low quality reads removing by cutadapt software (v1.9.3), the high quality clean reads were aligned to the reference genome (UCSC hg19) with hisat2 software (v2.0.4; Kechin et al., 2017). Then, guided by the Ensembl gtf gene annotation file, cuffdiff software (part of cufflinks) was used to get the gene level FPKM as the expression profiles of mRNA, and fold change and p-value were calculated based on FPKM, differentially expressed mRNA were identified (Trapnell et al., 2010).

MeRIP-Seq, Data Analysis, and Bioinformatics
At 12 h after operation, rats were anesthetized with pentobarbital sodium and were transcardially perfused with 20 ml of 4°C PBS. The retina was then rapidly dissected and stored in liquid nitrogen. Total RNAs from the retinal tissue of each rat were isolated using TRIzol reagent (Invitrogen) according to the manufacturer's protocol. In this study, three biological replicates were conducted in the sham-operated and TON model group, and retinal tissue samples from four rats were combined into one sample to attain 60 μg of total RNA. mRNA extraction was performed using NEBNext rRNA Depletion Kit (New England Biolabs) following the manufacturer's protocol. Briefly, m6A RNA immunoprecipitation was performed with the GenSeq ™ m6A RNA IP Kit (GenSeq Inc., shanghai, China) by following the manufacturer's instructions. Both the input sample without immunoprecipitation and the m6A IP samples were used for RNA-seq library generation with NEBNext ® Ultra II Directional RNA Library Prep Kit (New England Biolabs). The quality control (QC) of the libraries was evaluated using BioAnalyzer 2,100 system (Agilent Technologies). Library sequencing was performed on an Illumina Hiseq4000 platforms under 150 bp pairended resolution.
Cutadapt software (V1.9.3) was used to remove low quality reads and obtain clean reads (Kechin et al., 2017). The clean reads from all samples were then aligned to Ensembl reference genome by HISAT2 software (v2.0.4; Kim et al., 2015). Modelbased Analysis of ChIP-Seq(MACS) software (v1.4.2) was used to identify methylated genes in each sample (Zhang et al., 2008). DiffReps software was used for identification of differential methylated genes (Shen et al., 2013). The differentially methylated peaks (fold changes ≥ 2.0 and p < 0.05) between the shamoperated group and TON model group were analyzed by exomePeak and annotated accordingly by the latest Ensembl database. Gene ontology (GO) database 1 and the latest Kyoto Encyclopedia of Genes and Genomes (KEGG) database 2 were used to perform GO analysis and pathway enrichment analysis on differentially methylated mRNA genes.

Statistical Analyses
Data are expressed as mean ± standard deviation (SD). Statistical analyses were conducted using GraphPad Prism Version 8.3 software (GraphPad Software LLC, CA, United States). Student's t-test (paired) was used for comparing the statistical significance between two groups. For each analysis, p < 0.05 was considered as statistically significant and marked with *. 1 www.geneontology.org 2 www.genome.jp/kegg

m6A-Related Genes Were Upregulated in Retina After TON in Rats
Through qRT-PCR assay, the expression levels of four major m6A-related enzymes including METTL3, WTAP, FTO, and ALKBH5 in the retina were assessed at 12 h after the establishment of the TON model. Twelve hours after TON, mRNA levels of methyltransferases (METTL3 and WTAP) and demethylases (FTO and ALKBH5) were both significantly upregulated in the retina (Figures 1A,B; METTL3, FTO, and ALKBH5 p < 0.05; WTAP p < 0.01).

Profiles of m6A Modification in Retina After TON
The alterations of m6A-tagged mRNA were assessed via genomewide screening after TON in rats (Supplementary Table S1). Twelve hours after TON, a total of 2,810 m6A peaks were differentially upregulated, while 689 m6A peaks were downregulated (fold changes ≥ 2.0 and p < 0.05). The top 20 differentially expressed m6A peaks are listed in Table 2. These m6A peaks were enriched in coding sequence (CDS; Figures 2A,B). The differentially expressed m6A peaks can be found on all chromosomes, mainly on chr1, chr2, chr3, and chr10 ( Figure 2C). Moreover, the m6A peaks were both mainly characterized by GGACU motif (Figure 2D).

GO Analysis and Pathway Analysis of Differentially m6A-Modified mRNA
GO analysis and KEGG analysis were conducted to investigate the pathophysiological significance of m6A modification. According to the GO analysis, the upregulated peaks in the TON retina were significantly associated with protein binding (ontology: molecular function), nervous system development (ontology: biological process), and intracellular (ontology: cellular component; Figure 3A). The downregulated peaks were mainly related to unfolded protein binding (ontology: molecular function), cellular localization (ontology: biological process), and intracellular part (ontology: cellular component; Figure 3B). KEGG analysis showed that upregulated m6A peaks were remarkably associated with MAPK signaling pathway, NF-κB signaling pathway, and TNF signaling pathway ( Figure 4A). The downregulated m6A peaks were significantly correlated with ribosome pathway (Figure 4B). The enrichment of genes in the four major pathways is shown in Figure 5.

Overview of mRNA Expression Profiles and Conjoint Analysis of MeRIP-Seq & RNA-Seq
Through RNA-seq, the expression profiles of altered genes were determined (Supplementary Table S2). A total of 520 up-regulated genes and 258 down-regulated genes were detected (fold change ≥ 2.0, p ≤ 0.05; Figures 6A,B). The top 20 genes were shown in Table 3. The degree of differentially expressed genes between the two groups were analyzed by hierarchical cluster analysis (Figure 6C). By conducting conjoint analysis of MeRIP-seq & RNA-seq, 176 hypermethylated m6A peaks were identified in 161 upregulated mRNA transcripts (hyperup) and 15 downregulated transcripts (hyper-down); 63 hypomethylated m6A peaks were identified in 13 upregulated mRNA transcripts (hypo-up) and 50 downregulated transcripts (hypo-down; Figure 6D).

DISCUSSION
m6A is the most common methylation modification in eukaryotic cells, accounting for more than 80% of methylated sites in RNA (Bodi et al., 2010). By regulating the eukaryotic transcriptome, m6A modification can affect the stability, splice, export, and translation of mRNA transcription. A variety of studies have shown that dysregulation of RNA methylation could be associated with many biological processes, including neuron development and neurodegenerative diseases . Besides, m6A modification also plays an important role in the process of self-renewal and differentiation of neural stem cells, which may promote stem cell treatment and gene therapy for neurological diseases (Boles and Temple, 2017;Zhou et al., 2018). A recent study conducted genome-wide comparison of m6A-tagged transcript in the hippocampus of mice after TBI, and a total of 922 m6A modification sites were detected, among which 370 were up-regulated and 552 down-regulated (Wang et al., 2019). However, the role of m6A modification in TON has not been characterized. To explore the potential function of differentially expressed m6A modifications in TON, we performed MeRIP-seq to screen the profiles of m6A-tagged transcript in the retina of TON at the early stage. First, using qRT-PCR, we found that m6A-related genes (METTL3, WTAP, FTO, and ALKBH5) were both upregulated at 12 h in the retina of rats after TON. Then, MeRIP-seq was performed to determine the profiles of genome-wide m6A-tagged transcript at 12 h after TON. Briefly, 2,810 m6A peaks were differentially upregulated and 689 m6A peaks were downregulated. Furthermore, GO analysis and pathway analysis were conducted to predict the potential role of m6A modification in the pathophysiological process after TON. Finally, we applied conjoint analysis of MeRIP-seq &RNA-seq. The "write" of m6A methylation, namely, the addition of methyl groups to the N6 position of adenosine, is catalyzed by a large methyltransferase complex whose core region contains two proteins, the heterodimer of METTL3 and METTL14. There is mounting evidence to implicate the importance of METTL3 in the nervous system. METTL3 knockout may lead to early embryonic death and impair the formation of mature neurons in embryoid bodies (Angelova et al., 2018). METTL3 mediated m6A modification may be involved in cerebellar development by regulating the stability of concerned mRNA . As a subunit of m6A methyltransferase, WTAP regulates the catalytic activity of methyltransferase complex by interacting with METTL3 and METTL14 (Ping et al., 2014). WTAP is over-expressed in the glioblastoma stem cells, controlling the invasion of glioblastoma cells via regulating epidermal growth factor (EGF) signaling (Jin et al., 2012;Xi et al., 2016). FTO and ALKBH are two major m6A RNA demethylases. FTO was firstly recognized as a member of Fe (II)-and oxoglutarate-dependent AlkB dioxygenase family and related to adipogenesis (Zhao et al., 2014). The deficiency of FTO could restrain the differentiation of neural stem cells in vivo, impairing learning and memory (Li et al., 2017). Decreasing FTO via microinjection of trained herpes simplex virus (HSV) vector could improve contextual fear memory (Walters et al., 2017). Li et al. (2018) have found that FTO is involved in insulin defects-related Alzheimer's disease mice through activating TSC1-mTOR-Tau signaling. In the neurons of mice, m6A RNA demethylase ALKBH5 was recently found to be widely expressed (Du et al., 2020). During the development of the brain, ALKBH5 decreased dramatically. In brief, these studies indicated that m6A modification methyltransferases and demethylases (METTL3, WTAP, FTO, and ALKBH5) have much to do with the development of nervous system and neurological function such as learning and memory. In the present study, METTL3, WTAP, FTO, and ALKBH5 showed increased expression in the retina of mice at 12 h after TON. However, the exact roles of methyltransferases and demethylases in the pathophysiological process still need further investigation. Nowadays, m6A modification profiles in the neurological diseases have drawn more and more attention. Liu et al. (2020) have characterized the landscape and distribution patterns of m6A and m6Am modification in the tissues of humans and mice. They identified that 594 out of 21,480 m6A peaks were tissue specific, among which the brain-specific m6A signals were closely related to the functions of head development. In transient focal ischemia model of Adult C57BL/6J mice, m6A levels were found to be increased significantly at 12 h and 24 h of reperfusion. m6A peaks were upregulated in 139 transcripts (122 mRNAs and 17 lncRNAs) and downregulated in8 transcripts (5mRNAs and 3 lncRNAs) at 12 h after reperfusion (Chokkalla et al., 2019). Yu et al. (2020) found that 2,165 m6A peaks were significantly changed in the cerebral cortex of rats after TBI, of which 1,062 were upregulated and 1,103 were downregulated. According to their results, functional FTO is of great necessity to the repair of TBI-related neurological damage. In the present study, m6A modification peaks were screened by MeRIP-seq, 2,810 m6A peaks were upregulated and 689 m6A peaks were downregulated. GO analysis showed that the altered m6A peaks mainly correlated with nervous system development, indicating that m6A RNA methylation may be crucial to the pathophysiological process after TON. On the basis of pathway analysis, m6A-tagged transcripts were mainly related to MAPK signaling pathway, NF-κB signaling pathway, protein processing in endoplasmic reticulum and spliceosome. The MAPK pathway has been proven to be related to nerve injury after TBI, and the Sarm1-MAPK pathway can disturb the energy balance of axons, leading to the exhaustion of adenosine triphosphate (ATP) before complete axon injury, and promoting the progression of axon injury (Yang et al., 2015). A key pathway of inflammation in central nervous system injury is the NF-κB pathway, which regulates the expression of pro-inflammatory and pro-apoptosis genes in its active form (Howell and Bidwell, 2020). Neuroprotective effect may be achieved by regulating MAPK and NF-κB signaling pathways and thus regulating inflammatory responses .
By conducting conjoint analysis of MeRIP-seq & RNA-seq, we identified the transcripts which were differentially expressed as well as hypermethylated or hypomethylated. The results of the present study revealed that m6A modifications may provide a novel insight into TON.

CONCLUSION
In conclusion, our study indicated that the expression of methyltransferases and demethylases (METTL3, WTAP, FTO, and ALKBH5) are both up-regulated at 12 h after TON. Through MeRIP-seq and subsequent bioinformatics analysis, the potential functions of differentially expressed m6A modified transcripts were predicted. Furthermore, altered mRNAs with

DATA AVAILABILITY STATEMENT
The RNA-seq and MeRIP-seq data sets have been uploaded to the Gene Expression Omnibus repository (GEO accession number: GSE165520). All the records are now publicly-available.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Ethics Committee of Naval Medical University (Shanghai, China).

AUTHOR CONTRIBUTIONS
XQ: drafting/revision of the manuscript for content, including medical writing for content, major role in the acquisition of data, study concept or design, and analysis or interpretation of data. KZ: drafting/revision of the manuscript for content, including medical writing for content, major role in the acquisition of data, and analysis or interpretation of data. ZL: drafting/revision of the manuscript for content, including medical writing for content, and analysis or interpretation of data.
DZ: revision of the manuscript for content and including medical writing for content. LH: revision of the manuscript for content, including medical writing for content, and analysis or interpretation of data. All authors contributed to the article and approved the submitted version.