Transcriptome-Wide N6-Methyladenosine Methylome Alteration in the Rat Spinal Cord After Acute Traumatic Spinal Cord Injury

Recent studies showed that RNA N6-methyladenosine (m6A) plays an important role in neurological diseases. We used methylated RNA immunoprecipitation sequencing (MeRIP-Seq) technology to generate the m6A modification map after traumatic spinal cord injury (TSCI). A total of 2,609 differential m6A peaks were identified after TSCI. Our RNA sequencing results after TSCI showed 4,206 genes with significantly altered expression. Cross-link analysis of m6A sequencing results and RNA sequencing results showed that 141 hyper-methylated genes were upregulated, 53 hyper-methylated genes were downregulated, 57 hypo-methylated genes were upregulated, and 197 hypo-methylated genes were downregulated. Among these, the important inflammatory response factor Tlr4 and the important member of the neurotrophin family Ngf were both upregulated and hyper-methylated after TSCI. This study provides that in the future, the epigenetic modifications of the genes could be used as an indicator of TSCI.

a lack of effective pharmacotherapies due to our incomplete understanding of the physiological and pathological mechanisms on secondary injuries after TSCI. As currently known, the pathogenesis of TSCI can be divided into primary and secondary injuries (McDonald and Sadowsky, 2002;Dietz and Fouad, 2014;Ahuja et al., 2017). The primary injury immediately causes mechanical rupture and dislocation of the spine, resulting in compression or transection of the spinal cord. The secondary injury cascade then occurs within minutes of the primary injury and lasts for months or more than half a year (Pineau and Lacroix, 2007). The mechanism of secondary injury is complex and is not yet fully understood. Possible mechanisms include ischemia, hypoxia, inflammation, edema, excitatory toxicity, ion steady state imbalance, and apoptosis (Schanne et al., 1979;Kwon et al., 2004;Brockie et al., 2021). Moreover, current research on the treatment of spinal cord injury mainly focuses on treatments such as neuroprotective therapy, nerve regeneration therapy, cell transplantation, neuromodulation, and robotics (Hall and Braughler, 1982;Hurlbert et al., 2015). Unfortunately, there is still a lack of recognized clinically approved therapeutic drugs for the treatment of the subacute phase of TSCI (Kwon et al., 2004).
Post-transcriptional modification has become recognized as an important regulatory factor in various physiological and pathological processes and has attracted increasing attention in biological research (Fu et al., 2014;Meyer and Jaffrey, 2014). The N6-methyladenosine (m6A) modification is the most abundant mRNA modification (Geng et al., 2020). A methyltransferase (Writer), a demethylase (Eraser), and a reader protein (Reader) are all required for the regulation of the m6A of target genes (Zhuang et al., 2019;Qin et al., 2020;Zhang H. et al., 2020;Ye et al., 2021). The level of modification affects the stability and translation efficiency of a target gene (Ma et al., 2019). In recent years, research has focused on m6A modification in neurological diseases. It has been reported that m6A modification alleviates ischemic stroke by promoting Akt phosphorylation by disrupting the stability of Pten mRNA (Zhang Z. et al., 2020). During cerebral ischemia and reperfusion, ALKBH5 and FTO selectively demethylate Bcl2, preventing degradation of Bcl2 transcripts. These results highlight that m6A modification plays an important role in neurological diseases (Du et al., 2020;Xu H. et al., 2020;Xu K. et al., 2020). Downregulation of m6A mRNA methylation was found to be involved in Parkinson's disease . Oxygen glucose deprivation leads to an increase in the level of m6A modification of Lnc-D63785, which leads to the accumulation of miR-422a and induces neuronal cell apoptosis (Xu S. et al., 2020). In the early stage of acute ischemic stroke, the level of METTL3 increases, which promotes an increase in the methylation and expression level of miR-355, thereby promoting stress granule formation (Si et al., 2020).
However, the study of m6A modification after TSCI in mammals has not been reported. In this study, we obtained TSCI and sham spinal cord tissues for Methylated RNA Immunoprecipitation Sequencing (MeRIP-Seq). Moreover, we investigated the function of mRNAs with significantly changed peaks by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. Subsequently, a conjoint analysis was performed to reveal the significantly changed mRNAs with significantly altered enrichments that may play a vital role in the epigenetic modification of TSCI. Thus, we provide data for the exploration of TSCI pathological mechanisms and the development of effective treatment drugs.

Animals
All processes in this study followed the guidelines of the National Institutes of Health Laboratory Animal Care and Use Guidelines, and the study was approved by the Institutional Animal Care and Use Committee of Wuhan University (IACUC: ZN2021119). The male Sprague-Dawley rats (7 weeks) used in the study were purchased from Vital River Laboratory Animal Technology Co. Ltd. (Vital River, Beijing, China) and housed at the Animal Experiment Center of Zhongnan Hospital at Wuhan University. All the rats were placed in a regular environment with a 12 h light/dark cycle for 1 week before experiments with a weight ranging from 250 to 300 g.

Animal Model of Traumatic Spinal Cord Injury
Rats were anesthetized with an intraperitoneal injection of 5% pentobarbital (50 mg/kg) and placed in the prone position. Next, after skin preparation and positioning, we performed laminectomy to expose the T10 spinal cord (Lee et al., 2014). An impact device (6800 II; RWD Life Science Corp, Shenzhen, China) was used to create an TSCI with an energy of 25 gcm. Wagging tail reflection, body trembling, and lower limb retraction were observed after impact. After surgery, the rats were injected with antibiotics for 3 days and were helped to empty their bladder until bladder function was restored.

Hematoxylin and Eosin Staining
The T7-T12 spinal cords of rats were dissected 14 days after trauma. All the spinal cords were fixed in a 4% paraformaldehyde fix solution and then embedded in paraffin. Next, longitudinal serial sections (5 µm) from paraffin blocks were taken, and histopathological examination with Hematoxylin and Eosin (HE) staining was performed under a light microscope (Olympus BX53).

Basso-Beattie-Bresnahanlocomotor Rating Scale Measurements
Basso-Beattie-Bresnahan (BBB) motor function scoring (Basso et al., 1995) was performed 1, 7, and 14 days post surgery to evaluate hind limb motor function. Each test was double-blinded by three professionally trained research graduate students in the laboratory. Before the experiment, the motor function of each rat was determined to be normal. For each test, the rats were placed in an open basin, and a researcher gently tapped the wall of the basin, making it crawl, and observing the hip, knee, ankle joint walking, trunk movement, and coordination. The total score was 21 points, and the higher the score, the more perfect the motor function.

Inclined Plane Test
An inclined plane experiment (Kazanci et al., 2017) was carried out at 1, 7, and 14 days after the operation to evaluate the balance ability of each rat. A rat was placed on a plane equipped with an indicator indicating the inclination of the plane. Then, the plane was tilted at a speed of 2 • per second, and the angle at which the rat fell from the plane was recorded as the falling degree. Each test was repeated three times, and the average value was recorded.

Quantitative Reverse Transcription Polymerase Chain Reaction
The T7-T12 spinal cords of rats were dissected 14 days after trauma and immediately frozen in liquid nitrogen and stored at −80 • C. Total RNA was extracted and purified using the RNAeasy TM Animal RNA Isolation Kit with Spin Columns (R0024). Next, we analyzed the quality of RNA using enzymelabeled instruments, and samples with A260/A280 ratios between 1.9 and 2.2 were used for further experiments. The primers used in quantitative reverse transcription polymerase chain reaction (qRT-PCR) were designed using the Primer3 website 1 and listed in the additional file (Supplementary Table 1). Following the Biomarker One Step SYBR Green RT-qPCR Kit (RK02012, Biomarker Technologies) to construct the reaction system. After centrifugation and blending (D1008, SCILOGEX), qPCR reaction was performed using QuantStudio 1 real-time fluorescence quantitative PCR system (QuantStudio 1, Applied Biosystems). The original CT data were exported and input by Excel, and the data were calculated and analyzed by Delta CT method.

Methylated RNA Immunoprecipitation-QPCR (MeRIP-QPCR)
Take out the Invitrogen Dynabeads Antibody Coupling Kit (14311D, LC BIO) and prepare magnetic beads as instructed after balancing at room temperature for half an hour. Blend magnetic beads with Synaptic System M6A-Antibody (202003, LC BIO). Add 2 g RNA, 100 ml m6A binding buffer and add ddH 2 O to supplement to 500 ml. The incorporation of these fragment RNA with m6A-immunomagnetic beads at room temperature for 1 h with vertical mixing machine for 7 cycles/min. Place EP tube on magnetic holder for at least 5 min. Then, the preheated 125 µl elution buffer was added, slowly blown and mixed, and incubated at 42 • C for 5 min. After 5 min, blow again slowly and drain the supernatant to a new 1.5 ml EP tube. Repeat 3 times to elute the RNA fragments on the magnetic beads. The previous Elution buffer is mixed together and the EP tube is placed on ice. After 4 rounds of elution, the solution in the final EP tube was 500 ml. Use the RNeasy MinElute Cleanup Kit (74204, Qiagen) to purify RNA. Finally, a qPCR quantitative experiment was carried out with the Biomarker One Step SYBR Green RT-qPCR Kit (RK02012, Biomarker Technologies). Using Ct of control input, Ct of control IP, Ct of treatment group input, and Ct of control IP to calculate the changing trend of the peak.

Methylated RNA Immunoprecipitation Sequencing
Using a mixture of three rats as a biological replicate, approximately 10 µg of double poly(A) selected RNA was generated from 400 µg of total RNA. Then, the poly(A) RNA was fragmented into 100 nucleotide long pieces using a thermocycler. m6A-methylated RNA fragments were then enriched by immunoprecipitation with an anti-m6A antibody (Synaptic Systems, Cat. No 202 003). Then, 100 ng of RNA (100 ng of input and 100 ng of post-m6A-IP positive fraction) was used to construct libraries utilizing the Illumina TrueSeq Stranded mRNA library prep kit. Using the KAPA library quantification kit (KAPABIOSYSTEMS Cat.NO KK4824) to quantify the library, we then submitted these libraries for highthroughput sequencing (Illumina Hiseq X10; Meng et al., 2014;Lan et al., 2019).

Data Analysis
In order to obtain high-quality clean reads that could be used for subsequent analysis, Trimmomatic software (v0.36) was used to remove adapters from FASTQ data. We then used the default parameters of HISAT2 software (v2.1.0) to compare clean reads with the reference genome (vRnor_6.0) of Sprague-Dawley rats and retained only comparison reads for subsequent analysis. The Guitar R package (v1.1.1.18) was used to check the quality of the genome comparison result data. We then used MeTDiff software (v1.1.0), and an input sample was used as a control for peak detection. We also used MeTDiff software to detect significantly enriched peaks. Then, GO and KEGG enrichment analyses were performed using hypergeometric distribution testing. Motifs for peak sequences were generated using MEME-ChiP (v5.0.5). The input of MeRIP-seq was used as the RNA-seq data. DESeq2 software (V1.14.1) was used to standardize the count number of each sample gene (Basemean value was used to estimate the expression level), the multiple of difference was calculated, and NB (negative binomial distribution test) was used to test the significance of reads number. Finally, the differential protein coding genes were screened according to the difference multiple and difference significance test results.

Basic Characteristics of N6-Methyladenosine Modification After Traumatic Spinal Cord Injury
In order to evaluate the success of TSCI modeling, on the 14th day after surgery, samples were taken for HE and Nissl staining ( Figure 1A). Compared with the control group, a large number of spinal cord nerve cells had died and tissues had atrophied after TSCI. Moreover, we performed BBB scoring and inclined plane experiments on 1, 7, and 14 days after surgery to verify the success of this model (Figures 1B,C). The BBB score and inclination angle of TSCI group were lower than those of sham group, indicating a significant decline in posttraumatic exercise ability. We used the tissues (T7-T12) obtained 14 days after TSCI to conduct transcriptome-wide m6A sequencing (m6A-seq) and RNA sequencing (RNA-seq). A total of 90.17 G of clean bases were obtained from six samples (each sample was mixed from two spinal cord tissues). The effective data volume of each sample was distributed from 5.37 G to 5.6 G, Q30 scores ranged from 90.32 to 95.3%, and the average GC content was 52.89% (Supplementary Table 2). By comparing the reads to the reference genome, the genome comparison of each sample was obtained, and the mapping rate was 87.23-93.32% (Supplementary Table 3). And reads are distributed throughout the genome (Figure 2A). A total of 22,701 peaks were identified in the TSCI group, and 24,108 peaks were identified in the sham group. There was only one peak for most genes, most of them were downregulated, and there were at most seven peaks for one gene (Figure 2B and Supplementary Table 4). Moreover, 43.35% of peaks were enriched in the 3 UTR, 5% of peaks were enriched in the 5 UTR, and 43.99% of peaks were enriched in the exon in the Sham group. Furthermore, 48.8% of peaks were enriched in the 3 UTR, 10.13% of peaks were enriched in the 5 UTR, and 41.07% of peaks were enriched in the exon in the TSCI group ( Figure 2C). Our study identified the top five motifs to be "RGAAR, " "DCTGTR, " "AWATAH, " "GCTGGGRW, " and "CACRK" (Figure 2D). Figure 2E shows the peak of Tlr4 (3 UTR) and Ngf (3 UTR) in the sham and SCI groups.

Peak Alterations After Traumatic Spinal Cord Injury
To investigate m6A methylation alteration after TSCI, differential m6A peaks were further analyzed. We identified an average of Frontiers in Neuroscience | www.frontiersin.org 22,701 peaks in the TSCI group and an average of 24,108 peaks in the sham group. The average length of peaks in the TSCI group was 4,696.27 bp, and that in the sham group was 4,476.93 bp. By comparing the peaks in the TSCI and sham groups, 1,948 hypomethylated m6A peaks and 661 hyper-methylated m6A peaks were subsequently identified (p < 0.05, FC > 1.5; Figure 3A). A volcano chart showed the p-value and multiple changes of these differential peaks ( Figure 3B). Moreover, the top 20 altered m6A peaks are listed in Table 1. The information on all the different genes is listed in Supplementary Table 5. We further investigated the location of significantly changed peaks and found that 1,537 differential peaks were located in the 3 UTR and exon regions, 438 differential peaks were located in the intron and exon regions, 255 differential peaks were located only in introns, and 192 differential peaks were located only in exons. A total of 122 differential peaks were located in the 5 UTR, intron, and exon regions, 59 differential peaks were located in the 5 UTR and exon regions, and 6 differential peaks were located in the 5 UTR, intron, 3 UTR, and exon regions ( Figure 3C). In order to verify the accuracy of the sequencing results, we selected the hyper-methylated (Ngf, Tlr4) and hypo-methylated (Sp3, Chil1) genes after TSCI for MeRIP-qPCR, and the results showed good consistency ( Figure 3D).

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analyses
Gene ontology and KEGG analyses were carried out using the genes with differential peaks to determine biological functions or pathways that these genes with differential peaks were mainly associated with. GO analysis showed that the differential peaks were related to the biological processes such as sterol biosynthetic process, RNA splicing, regulation of stress fiber assembly, and regulation of resting membrane potential. They were related to cellular components such as spliceosomal complex, nucleus, nucleoplasm, and nucleolus. Molecular functions such as ubiquitin-protein ligase activity, RNA binding, and protein kinase A binding were related to these significantly changed peaks ( Figure 4A). KEGG analysis indicated that the differential peaks were concentrated mainly in cellular processes such as tight junctions, signaling pathway regulating pluripotency of stem cells, regulation of actin cytoskeleton, and p53 signaling pathway, and these were related to environmental information processing such as Wnt signaling pathway, TGF-beta signaling pathway, and Rap1 signaling pathway. There were also relations to genetic information processing such as ubiquitin-mediated proteolysis, spliceosome, RNA transport, and RNA degradation. These genes were also related to diseases such as transcriptional misregulation in cancers and to primary immunodeficiency and metabolism such as ubiquinone and other terpenoid-quinone biosynthesis and terpenoid backbone biosynthesis. In addition, these genes were related to organic systems such as Toll-like receptor signaling pathway, T-cell receptor signaling pathway, and retrograde endocannabinoid signaling pathway ( Figure 4B).

RNA Expression Profiles After Traumatic Spinal Cord Injury
By using the input library, this study generated a differential gene expression map of TSCI vs. sham, with a total of 1,720 downregulated genes and 2,486 upregulated genes (p < 0.05, FC > 1.5; Figure 5A). In the results of the principal component analysis, confidence ellipses separated the TSCI and sham groups, indicating that there was comparability between these two groups ( Figure 5B). The heat map shows the relative expression levels of each sample in the TSCI and sham groups; samples in the  same group showed similar patterns, and the two groups showed obvious different patterns ( Figure 5C). The volcano graph shows the overall P value and foldchange of each difference peak ( Figure 5D). Upregulated genes such as Mmp12, Ifitm1, and Cd68 and downregulated genes such as Fut4, Lrrc17, and Dsc2 are marked in this graph. Moreover, the top 20 differentially expressed genes are listed in Table 2. All genes with significant changes in expression after TSCI were listed in Supplementary  to conduct the validation experiment. The qRT-PCR results showed good consistency with the sequencing result ( Figure 5E).

Crosslinking Analysis Results
Cross-linking analysis of differential peaks and differential expression of given genes was performed (Supplementary Table 7), and the four-quadrant graph shows the FC of differential peaks and differential RNAs (Figure 6A). A total of 141 hyper-methylated genes were upregulated, 53 hypermethylated genes were downregulated, 57 hypo-methylated genes were upregulated, and 197 hypo-methylated genes were downregulated in the two groups ( Figure 6B). Using the significantly changed genes with differentially methylated peaks, a protein interaction network diagram was created, wherein the different colors represent different gene expression and methylation trends. As shown in Figure 6C, the hypermethylated genes Ptprc, C3ar1, Myd88, and Cdc20 were upregulated after TSCI. The hypo-methylated genes Grm3, Gria4, and Gria2 were downregulated after TSCI. The hyper-methylated genes Brnp2K, Pld1, and Gpr65 were downregulated after TSCI. The hypo-methylated genes Gria1, Agap2, and Kcnj9 were upregulated after TSCI.

DISCUSSION
An increasing number of studies have found that RNA m6A methylation plays an important role in central nervous system (CNS) injury disease and neurodegenerative diseases (Weng et al., 2018;Meng et al., 2019;Wang et al., 2021). However, there are currently no studies on how m6A is related to acute TSCI, partly because the methylation modification profile in TSCI is not yet clear. Using MeRIP-Seq, we have obtained a panoramic view of m6A methylation in the spinal cord after TSCI, which provides material for further in-depth study of m6A methylation in TSCI. A total of 2,605 significantly changed peaks were identified in this study, among which 1,948 peaks were downregulated and 661 were upregulated. A summary of the peak sites indicated that the m6A peaks with different methylation in traumatic brain injury (TBI) were mainly distributed at the 3 UTR near the stop codon and at the 5 UTR near the start codon, which was consistent with the characteristics found in other sequencing studies. Peaks identified in a MeRIP-Seq conducted on the hippocampus of mice following TBI were enriched mainly in the coding sequencing near the stop codon . Moreover, research conducted on the cerebral cortex of rats following TBI showed a similar peak distribution pattern as that observed in our study .
We also conducted GO and KEGG analyses, and the most enriched pathway and the biological process had some similarities with other NCS injury disease datasets. The MeRIP-Seq conducted on mice hippocampus indicated that the significantly changed peaks were enriched in the regulation of actin cytoskeleton, rap1 signaling pathway, and MAPK signaling pathway, which were also the most enriched pathways in our research . A study conducted on C56BL/6J mice with focal ischemia revealed that significantly changed peaks were also enriched mostly in the p53, NF-KB, and MAPK signaling pathways, as was found in our study (Chen et al., 2018;Chokkalla et al., 2019;Ye et al., 2020). However, the m6A sequencing conducted on the neurodegenerative diseases showed different characteristics. In Alzheimer's disease (AD), the significantly changed peaks were mostly enriched in synaptic growth at the neuromuscular junction, snRNA transcription, and the smoothened signaling pathway involved in dorsal/ventral neural tube patterning (Han et al., 2020). This phenomenon might be due to the fact that our sample was taken in the subacute phase of TSCI, while neuro degeneration mostly occurs during the chronic phase. Our results from cross-linking analysis showed that some TSCI-related genes were changed significantly and hypo-or hyper-methylated after TSCI. As is well known, all main CNS cell types express toll-like receptor 4 (TLR4). Moreover, our study found that TLR4 was upregulated and hypermethylated after TSCI. We found that TLR4 signaling is essential for maximal sparing and replacement of acute oligodendrocyte lineage cells after spinal cord injury (Church et al., 2016). Studies have revealed that FGF10 inhibits the activation and proliferation of microglia/macrophages by regulating the TLR4/NF-B pathway and inhibits the release of pro-inflammatory cytokines after TSCI (Chen et al., 2017). All these results indicate that TLR4 plays important mediating roles in the pathological process of TSCI, and this has been well demonstrated. However, RNA methylation has been previously underrecognized, and its role remains to be clarified. As a member of the neurotrophin family, nerve growth factor (Ngf ) has been demonstrated to modulate the survival of damaged neurons and facilitate axonal sprouting and regeneration (Iulita and Cuello, 2016;Indo, 2018;Hu et al., 2020). Moreover, our study found that Ngf was upregulated and hyper-methylated after TSCI. A study conducted in 2019 using a new delivery system to deliver Ngf to the spinal cord trauma site and the pathological and behavioral test results have confirmed that exogenous Ngf enables neural regeneration, tissue remodeling, and functional recovery in mice with TSCI. Also, as is well known, m6A can mediate mRNA expression by directly influencing its stability or by attracting readers to attach to mRNA. Some drugs [such as FB23-2 ] have been developed to affect the enzymes related to m6A, and whether these drugs can be used to treat TSCI by influencing the Ngf level is an interesting question.
The study of m6A modification in CNS injury and neurodegenerative diseases has also unveiled some potential intervention targets (Mathiyalagan et al., 2019;Pan et al., 2021). FTO and METTL14, as important m6A mediation enzymes, have been found to be downregulated in rat cerebral cortex following TBI, and FB23-2 (an FTO inhibitor) contributes to the recovery of neural function in TBI. A study conducted in C56BL/6J mice revealed that m6A levels were elevated globally after transient focal ischemia, and the FTO level was also downregulated; further measurements of the m6A profile after focal ischemia with microarrays indicated that 95% of the transcripts were significantly hypermethylated (Chokkalla et al., 2019). Researchers have revealed that METTL3 is upregulated and FTO is downregulated in AD. More importantly, although some research has focused on the m6A modification in CNS injury disease, the downstream mechanism of these significantly changed genes with peaks has not been illustrated, and it awaits further investigation. There is currently no research that has focused on creating an m6A methylation atlas of TSCI in mammals. Our research fills this gap and provides potential genes that can be considered as key targets for TSCI treatment.

CONCLUSION
This study demonstrates that the m6A modification and RNA expression levels change after TSCI, as determined by MeRIP-Seq. Some genes were identified in this study such as Tlr4 and Ngf showed alteration in both expression and methylation level, providing substantial indicator in the epigenomic study of TSCI.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA791459.

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee of Wuhan University (IACUC: ZN2021119).

AUTHOR CONTRIBUTIONS
YZ and XJ conceived and designed the research plans, supervised, and fixed the writing. ZZ, XZ, RL, and PW established the TSCI model and isolated RNA samples from cerebral cortex tissues. JY and HC completed the data processing, normalization, and bioinformatics analyses. JY and HC wrote the article with contributions from all the authors. All authors have read and approved the manuscript.

FUNDING
This study was supported by the Research Fund from Medical Sci-Tech Innovation Platform of Zhongnan Hospital, Wuhan University (PTXM2020001). The funder had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.