MA was obtained from Must Bio-Technology Co., Ltd. The chemical structure of MA (Figure 1A) underwent confirmation through nuclear magnetic resonance (NMR) spectral analysis (Bruker) (Figures 1B, C), and the purity of MA was 98.979% (Figure 1D), as determined by high-performance liquid chromatography (HPLC) (Thermo Fisher Scientific Inc.).

The PharmMapper server was used to predict the potential target proteins of MA. First, the 2D structure of MA was obtained from the PubChem database and a mol2 file was generated using energy minimization principles. This file was then uploaded to the PharmMapper server. The UniProt database was utilized to standardize the UniProt IDs to gene symbols. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of the predicted targets were enriched using DAVID (National Institute of Allergy and Infectious Diseases, NIAID).

HT22 cell lines were originally obtained from the National Collection of Authenticated Cell Cultures. They were cultured in Dulbecco’s modified Eagle’s medium (Basal Media) containing 10% fetal bovine serum (Gibco) and then maintained in an incubator under 5% CO₂ at 37°C. MA was dissolved in dimethyl sulfoxide (DMSO; Sangon), and the final working concentrations of MA were 0, 400, 800, and 1,600 μM.

HT22 cells at a density of 4 × 10³ cells per well in the logarithmic growth phase were plated on a 96-well plate and then exposed to 0, 400, 800, and 1,600 μM of MA for 72 h. Afterward, each well was supplemented with 10% CCK-8 solution (Yeasen) and incubated at 37°C for 2 h. Subsequently, the OD value at 450 nm was detected using a microplate reader (BioTek).

LDH release in the cell supernatant was assessed following the manufacturer’s instructions after stimulation with 0, 400, 800, and 1,600 μM of MA for 72 h using an LDH assay kit (Beyotime). The OD value at 490 nm was determined using a microplate reader (BioTek).

The HT22 cells were plated in a 24-well plate at a density of 6 × 10⁴ cells per well, treated with 0, 400, 800, and 1,600 μM of MA for 72 h, and then incubated with 1 μM 2,7-dichlorofluorescein diacetate (DCFH-DA; Solarbio) in a dark environment at 37°C for 30 min. Reactive oxygen species (ROS) accumulation was assessed by measuring the dichlorofluorescein intensity using a fluorescence microscope (ZEISS).

Protein concentration was determined using the BCA method. The levels of malondialdehyde (MDA) content and superoxide dismutase (SOD) activity in HT22 cells were measured after treatment with 0, 400, 800, and 1,600 μM of MA for 72 h using an MDA assay kit (Njjcbio) and an SOD activity assay kit (Njjcbio), respectively. The OD values were measured at 530 and 560 nm using a microplate reader (BioTek), respectively.

Given the results from the cell viability experiments, we selected the group treated with 800 μM MA with HT22 cell viability at approximately 70% for MeRIP-seq. The HT22 cells were plated in 6 cm dishes at a density of 2 × 10⁶ cells per dish, followed by exposure to 0 and 800 μM MA for 72 h before being lysed using TRIzol Reagent (Invitrogen). The enrichment and sequencing of m6A-modified RNA were carried out by Seqhealth Technology Co., Ltd. In brief, polyadenylated RNA was enriched by VAHTS mRNA Capture Beads (Vazyme) from 10 µg of total RNA. Subsequently, mRNA was fragmented into 100–200 nt fragments, with 10% of the RNA fragments saved as input and the remaining portion subjected to m6A immunoprecipitation with an anti-m6A antibody (Synaptic Systems). Finally, the m6A-modified mRNA fragments were purified using an RNA Clean and Concentrator Kit (Zymo Research) and were then utilized for library construction and sequenced on a NovaSeq 6000 platform (Illumina). The RNA-seq and MeRIP-seq data have been uploaded to the GEO database (ID: GSE261320).

The raw MeRIP-seq data were analyzed using fastp, followed by alignment to the mm10 genome reference sequences using HISAT2 (Kim et al., 2019). RNA expression levels and differential expression were analyzed using StringTie (Pertea et al., 2016) and DESeq2 (Love et al., 2014), respectively. The m6A peak calling and differential methylation detection were performed using exomePeak2 (Meng et al., 2013; Tang et al., 2021). The Poisson generalized linear model was used in exomePeak2 for estimating the methylation level and identifying differentially methylated regions, while sequencing depth size factors were estimated using exomePeak on non-methylated background regions. The m6A motif sequence was identified using STREME (Bailey, 2021) by integrating a position weight matrix Markov model into its algorithm. MetaTX was used for visualizing the distribution of epitranscriptome profiles (Wang Y. et al., 2021). GO and KEGG annotations were conducted using DAVID (Jiao et al., 2012). The m6A conservation and disease association data were obtained from ConsRM (Song et al., 2021) and RMDisease (Song et al., 2023a), respectively. The substrates of m6A regulators were identified using CLIP-seq datasets (Wang et al., 2020; Liu et al., 2022; Yin et al., 2022). The m6A patterns of gene methylated sites were visualized using Integrative Genomics Viewer (IGV) software. The PPI analysis was conducted using STRING, and the results were visualized using Cytoscape software. The gene set enrichment analysis (GSEA) and heatmaps were generated at https://www.bioinformatics.com.cn (last accessed on 20 February 2023).

The mitochondrial membrane potential (MMP) was determined using the JC-10 probe (Solarbio). The JC-10 probe accumulates in the mitochondrial matrix, forming red fluorescent aggregates. Upon the decrease in MMP, the JC-10 probe diffuses out of the mitochondria, changes to monomeric form, and stains the cell with green fluorescence. The HT22 cells were plated in a 24-well plate at a density of 6 × 10⁴ cells per well and treated with 0, 400, 800, and 1,600 μM of MA for 72 h. Subsequently, the cells were incubated with the JC-10 probe at 37°C for 20 min. Finally, red and green fluorescence were analyzed using a microscope (ZEISS).

The opening of the mitochondrial permeability transition pore (mPTP) was assessed using an mPTP Fluorescence Assay Kit (Beyotime) following the manufacturer’s instructions. Calcein AM enters the cell by passive transport and accumulates in the mitochondria. In cells, the barely non-fluorescent calcein AM is hydrolyzed by intracellular esterases to produce calcein, which has no membrane permeability, thereby retaining calcein in the cell and causing the cytoplasm to fluoresce in strong green. When calcein binds to Co²⁺, the fluorescence signal is quenched. Under normal conditions, the mPTP is closed, preventing Co²⁺ from entering the mitochondria. However, when mitochondria are damaged, the mPTP opens, allowing Co²⁺ to enter and bind with calcein, resulting in the quenching of green fluorescence.

Briefly, HT22 cells were exposed to 0, 400, 800, and 1,600 μM of MA for 72 h, followed by incubation with calcein AM staining solution for 45 min at 37°C. Subsequently, the staining solution was replaced with a cell culture medium, and incubation continued for another 30 min at 37°C. Finally, green fluorescence was analyzed using a microscope (ZEISS).

The intracellular calcium ion concentration was measured using the Fluo-4 AM probe (Beyotime). Upon entering the cell, Fluo-4 AM was enzymatically cleaved by esterases to produce Fluo-4, which could bind to calcium ions and produce strong green fluorescence. Following the treatment of HT22 cells with 0, 400, 800, and 1,600 μM of MA for 72 h, the cells were exposed to a working solution of Fluo-4 AM at 37°C for 45 min. Following this, the cells underwent three washes with PBS and were subsequently incubated for an additional 20 min. The resulting green fluorescence was then analyzed using a microscope (ZEISS).

Mitochondrial activity was detected using the MitoTracker Red CMXRos probe (Beyotime). The MitoTracker Red CMXRos probe contains mildly thiol-reactive chloromethyl that specifically marks biologically active mitochondria and emits red fluorescence. After a 72-h treatment with 0, 400, 800, and 1,600 μM of MA, the cells were exposed to MitoTracker Red CMXRos at 37°C for 15 min. Following the removal of the solution, a fresh cell culture medium was added, and the cells were analyzed using a fluorescence microscope (ZEISS).

HT22 cells were plated in 6 cm dishes at a density of 2 × 10⁶ cells per dish and treated with 0, 400, 800, and 1,600 μM MA for 72 h. Afterward, the HT22 cells were prefixed in a mixture of 1.5% paraformaldehyde and 3% glutaraldehyde for 2 h at 4°C and in a mixture of 1% osmic acid and 1.5% potassium ferrocyanide for 1 h at 4°C. The samples were dehydrated in an ethanol gradient solution and embedded in epoxy resin for staining and sectioning. Ultrathin sections were observed using an EM 208 transmission electron microscope (Philips, Amsterdam, Netherlands).

Quantitative real-time polymerase chain reaction (qRT-PCR) was utilized to detect the expression of genes. After treating HT22 cells with 0, 400, 800, and 1,600 μM of MA for 72 h and treating zebrafish larvae with 0, 30, and 60 μM of MA at 1 day post-fertilization (dpf) for 2 days, TRIzol Reagent (Invitrogen) was used to extract total RNA. Subsequently, the extracted RNA was reverse-transcribed into cDNA. qRT-PCR was performed using a MonAmp ChemoHS qPCR Mix Kit (Mona) according to standard protocols. Actb served as an internal reference gene. All PCR primer sequences are shown in Table 1.

HT22 cells were plated in a 6-well plate at a density of 5 × 10⁵ cells per well and treated with 0, 400, 800, and 1,600 μM MA for 72 h. Total protein was extracted by RIPA lysate containing phosphatase inhibitors and proteinase inhibitors (Beyotime). The proteins were separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) and then transferred onto PVDF membranes and blocked with 5% skimmed milk. Then, the membranes were incubated with primary antibodies for PINK1, LC-3B, p62, and β-actin (Proteintech) overnight at 4°C. Next, the membranes underwent three washes with TBST buffer and were co-incubated with secondary antibodies for 1 h. Finally, the membrane proteins were exposed to ECL ultra-high-sensitivity luminescent liquid, and images were obtained using an Amersham Imager 680 System (Cytiva).

The interaction between MA and m6A-modified differentially expressed regulators was investigated through molecular docking using SYBYL-X 2.0 software. The 3D structures of proteins were obtained from AlphaFold (Jumper et al., 2021) and pre-treated using SYBYL-X 2.0 to add hydrogen atoms, remove heteroatoms and water molecules, and repair side chains (Wu et al., 2022b). The 2D structure of MA was obtained from the PubChem compound database and subsequently transformed into 3D structures based on energy minimization principles (Zhao J. et al., 2020). The total docking score, which comprehensively evaluated solvation, entropy, hydrophobic complementarity, and polar complementarity, was considered indicative of a stable interaction between the proteins and molecules when it exceeded 5 (Zhu et al., 2019).

To simulate the protein–ligand interaction, we utilized the Simulation Package tOward Next GEneration (Huang et al., 2022). The FF14SB united-atom force field (Maier et al., 2015) was used consistently throughout the simulation study. The complex systems were solved using the SPC/E water model (Berendsen et al., 1987) in cubic boxes, ensuring a minimum distance of 1.2 nm from the cube edge. Potassium and chloride ions were added to neutralize the systems. For each system, molecular dynamics (MD) runs lasting 50 ns were executed. To evaluate binding stability in a dynamic environment, various MD trajectory analyses, including root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), H-bond occupancy, and total binding free energy, were applied to analyze the simulation results.

HT22 cells were treated with 0 and 800 μM MA in 6-cm dishes for 72 h. After washing the cells once with PBS, 350 μL of PBS was added to collect and re-suspend the cells. The cell suspension was divided into aliquots, with each tube containing 50 μL. The metal-bath temperature was set at 37, 41, 45, 49, 53, 57, and 61°C for heating the individual tubes for a duration of 3 min. Following heat treatment, the cells underwent three cycles of freeze–thawing using liquid nitrogen and were subsequently centrifuged at 12,000 ×g for 10 min. Following centrifugation, protein samples were prepared using the supernatant, and the interaction between METTL14 and MA was confirmed by Western blot.

Zebrafish (AB line) were raised at 28.5°C under a 14-h light/10-h dark cycle in a circulating system, and all procedures were conducted in compliance with the regulations of the Institutional Animal Care and Use Committee (IACUC) of Fujian Medical University. The Fujian Medical University IACUC reference number was FJMU2023-Y-0824. The zebrafish larvae were treated with 30 and 60 μM MA at 1 dpf for 2 days. Siblings were treated with 0.1% DMSO as a control. The survival rate and hatching rate were calculated at 2 days post-treatment (dpt), and the body length was measured at 1 dpt and 2 dpt, respectively, to evaluate the influence of MA treatment on zebrafish development. Zebrafish were anesthetized using a 0.02% tricaine solution before imaging.

Wild-type zebrafish larvae were exposed to 5 μM MA or 0.1% DMSO at 4 dpf, and locomotor behavior was assessed at 9 h post-treatment (hpt) using Noldus EthoVision^® XT video tracking software (version 10.1.856; Noldus Information Technology, Netherlands) and the Noldus DanioVision^® zebrafish tracking hardware system. Zebrafish larvae were placed in 12-well plates (one larva per well) at 28°C and acclimated for 10 min. The behavior of zebrafish larvae was then recorded for 5 min, followed by the analysis of the locomotor distance, duration, speed, and meander.

The experimental data were analyzed using SPSS and presented as the mean ± standard deviation. A one-way analysis of variance was used for group comparison, with a significance level set at p < 0.05.

We first used the PharmMapper server to predict the potential target proteins associated with MA and found 302 potential target proteins. In order to gain further insights into their functional characteristics, GO and KEGG enrichment analysis were conducted using DAVID. This analysis resulted in the identification of 471 significant GO entries and 142 significant KEGG entries, with a statistical level set at p < 0.05. The GO enrichment terms mainly included mitochondrion, lysosome, and oxidoreductase activity (Figure 1E). The KEGG enrichment terms mainly included drug metabolism-cytochrome P450, adherens junction, and autophagy-animal (Figure 1F).

We first used HT22 cells to examine the neurotoxic effects of MA in vitro. HT22 cells were treated with different concentrations of MA for 72 h. In the control group, HT22 cells exhibited a polygonal shape with well-defined boundaries (Figure 2A). The cell body demonstrated extensive branching growth and interconnectedness, forming a complex network structure. Following exposure to MA, the cellular morphology became rounder, more oval-shaped, and heteromorphic. Additionally, the interconnection of cell processes was reduced, and the network structure became less apparent. The CCK-8 results showed that the viability of HT22 cells decreased to 70.67% (p < 0.05) and 45.11% (p < 0.001) after treatment with 800 and 1,600 μM MA for 72 h, respectively (Figure 2B). The release of LDH was detected following MA treatment, revealing a dose-dependent induction of LDH release from MA-treated cells (p < 0.05), thereby indicating the cytotoxicity of MA toward HT22 cells (Figure 2C).

ROS are byproducts generated during normal cellular metabolism. Under physiological conditions, intracellular antioxidant systems predominantly eliminate ROS to maintain their low levels. However, when the intracellular antioxidant system is inhibited or there is an imbalance between the antioxidant system and the production of ROS, oxidative stress occurs. DCFH-DA staining results showed significant ROS accumulation after treatment with 400, 800, and 1,600 μM MA (p < 0.001), indicating that MA induces oxidative stress in HT22 cells (Figures 2D, E). Furthermore, the MA-induced oxidative damage was confirmed by an increase in the lipid peroxidation marker MDA, accompanied by a significant decrease in the activity of the antioxidant enzyme SOD (Figures 2F, G).

The process and specific mechanism of neurogenesis in zebrafish closely resemble those observed in mammals. Consequently, zebrafish serves as a valuable vertebrate model for investigating and assessing the neurotoxicity of drugs. In order to elucidate the in vivo neurotoxic effects of MA, we exposed 1-dpf zebrafish embryos to varying concentrations of MA. We observed a significant dose-dependent decrease in both survival and hatching rates at 2 dpt with MA (Figures 3A, B). Furthermore, pericardial edema and cardiac congestion were observed following MA treatment, consistent with previous reports of cardiotoxicity. Notably, surviving zebrafish exhibited a small-head phenotype and significantly reduced body length (Figures 3C, D). After treating zebrafish larvae with 0, 30, and 60 μM of MA at 1 dpf for 2 days, the mRNA expression levels related to neural development were measured. As shown in Figure 3E, MA treatment significantly decreased the mRNA expression levels of genes associated with neurodevelopment, such as elavl3, gap43, gfap, and mbpa. In order to investigate the potential interference of MA on the neuromotor system of zebrafish, we conducted behavioral analyses following MA treatment. The administration of MA resulted in a significant reduction in locomotor distance, movement duration, and speed (Figures 3F–I). Furthermore, larvae exposed to MA exhibited twitching movements characterized by enhanced meandering behavior, indicative of neuronal damage-related motor impairments (Figure 3J). The above results indicated that MA treatment inhibited neural development in zebrafish.

To explore the neurotoxicological mechanism of MA, we performed RNA-seq analysis on HT22 cells with or without MA treatment. Given the results from the cell viability experiments, we selected the group treated with 800 μM MA with HT22 cell viability at approximately 70% for mRNA-seq. Principal component analysis (PCA) showed that the cells treated with MA were segregated into distinct groups compared to the untreated cells, indicating a significant alteration in the cellular state following MA treatment (Figure 4A). A total of 2,291 differentially expressed genes (DEGs) were identified, comprising 1,126 upregulated and 1,165 downregulated genes after MA treatment (Figure 4B).

DAVID was utilized to perform KEGG and GO pathway analysis in order to predict the associated signaling pathways and biological functions related to the DEGs. The KEGG pathway analysis showed that the upregulated genes were mainly enriched in the lysosome, the mTOR signaling pathway, and autophagy (Figure 4C). GO enrichment analysis was categorized into three groups: biological process (BP), cellular component (CC), and molecular function (MF). The most enriched BP terms were associated with lysosome organization, autophagy, and phagosome–lysosome fusion (Figure 4D). CC terms were involved in the mitochondrion, proton-transporting v-type ATPase complex, and autophagosomes. The MF terms exhibited significant enrichment in proton-transporting ATPase activity, RNA binding, and oxidoreductase activity. In addition, the downregulated genes after MA treatment were found to be involved in several KEGG pathways, including the TNF signaling pathway, focal adhesion, PI3K-Akt signaling pathway, and cellular senescence (Figure 4E). The most enriched BP terms were associated with cell migration, regulation of cell proliferation, motor neuron axon guidance, and response to calcium ions. CC terms were enriched in the extracellular matrix, cell junction, and cytoskeleton. MF terms were involved in calcium ion binding, metal ion binding, and actin filament binding (Figure 4F).

Next, we performed GSEA to further elucidate the signaling pathways associated with MA treatment. Consistently, phagosome acidification and phagosome maturation were upregulated after MA treatment compared to control cells (Figures 4G, H). Additionally, RNA-seq data revealed that the expression level of autophagy-associated genes was upregulated in MA-treated HT22 cells, suggesting that MA induced autophagy in mouse neurons (Figure 4I).

In order to investigate the effects of MA treatment on protein expression, we carried out a quantitative proteomic analysis to elucidate comprehensive protein changes in HT22 cells with or without MA treatment. The PCA revealed a distinct segregation of MA-treated cells compared to the control cells, aligning with our RNA-seq findings (Figure 5A). Volcano plots from the proteomic analysis displayed a total of 113 differentially expressed proteins (DEPs) in the MA-treated group, with 70 upregulated and 43 downregulated proteins relative to the untreated control (Figure 5B). GO and KEGG pathway enrichment analysis showed that DEPs induced by MA treatment were mainly enriched in biological processes, including the mitochondrion, autophagosome assembly, lysosome organization, autolysosome, and autophagy (Figure 5C). The heatmap illustrates the expression patterns of DEPs associated with the aforementioned terms (Figure 5D).

The correlation among biological replicates within the same group was high, with the absolute values of R minima higher than 0.7 and p < 0.05, which validated the reliability of our model and the reproducibility of the sample processing approach (Figure 5E). Subsequently, we examined the regulatory relationship among these genes using STRING to construct a PPI network. MAP1LC3B, LAMP2, CTSD, LAMP1, and GABARAPL2 emerged as the top five proteins, exhibiting degrees of 13, 13, 12, 10, and 9, respectively (Figure 5F). GSEA results indicated an upregulation in the regulation of autophagy, autophagosomes, and lysosomes following MA treatment, aligning with our RNA-seq findings (Figures 5G–I).

Furthermore, we examined the concordance in the directions of change between mRNA and protein expression levels (Figure 5J). The correlation between the mRNA and protein (r = 0.1151 and p < 0.05) was positive and highly significant. We identified genes that exhibited concordant increases in both mRNA and protein (n = 17), discordant patterns with mRNA decreasing and protein increasing (n = 16), discordant patterns with mRNA increasing and protein decreasing (n = 3), and concordant decreases in both mRNA and protein (n = 32). A total of 19 mRNAs and proteins showed inconsistent trends, which may be attributed to RNA post-transcriptional modification and post-translational processing and modification. Genes that overlap between mRNA-seq and proteomic analysis are listed in Supplementary Table S1.

The results of RNA-seq and proteomics showed that MA exposure can affect the mitochondria. Therefore, we examined the functional status of the mitochondria. In our assessment of mitochondrial activities post-MA treatment, we observed that aggregates of JC-10, indicative of healthy and polarized mitochondria, decreased significantly, while monomeric JC-10, which signified membrane depolarization and dysfunction, increased (Figure 6A). Concurrently, the fluorescence of accumulated calcein in the mitochondria decreased, suggesting that MA treatment induced damage to the mPTP of HT22 cells (Figure 6B). The elevation of cytoplasmic calcium concentrations occurs through either the influx of extracellular calcium ions or the release of calcium ions from cellular organelles. This increase in cytoplasmic calcium levels facilitates the increased binding of intracellular Fluo-4 AM to calcium ions, resulting in a subsequent augmentation of green fluorescence intensity within the cell. Our findings indicate a dose-dependent increase in Fluo-4 AM fluorescence intensity post-MA treatment, reflecting an augmentation of mitochondrial calcium efflux (Figure 6C). In line with these observations, the fluorescence intensities of MitoTracker Red CMXRos showed a significant decrease, suggesting compromised mitochondrial activities in HT22 cells after MA treatment (Figure 6D).

Mitochondrial damage serves as a pivotal trigger for the induction of autophagy, which is a widely presented degradation/recycling system in eukaryotic cells and plays a crucial role in maintaining neuronal homeostasis. When stimulated, misfolded proteins or abnormal cell components are transported to the lysosome for degradation to maintain cellular self-renewal. In the TEM examination, the control group exhibited normal mitochondrial ultrastructure in HT22 cells. However, exposure to MA resulted in the dissolution of the mitochondrial membrane, the disappearance of mitochondrial cristae, and abnormalities in mitochondrial morphology across all MA-exposed groups. Autophagosomes, indicated by yellow arrows, are double-membrane vesicles formed during autophagy. They phagocytose a variety of intracellular substances and transport them to lysosomes. Subsequently, autophagosome–lysosome fusion occurs to generate autolysosomes for the degradation of these substances. Furthermore, there was an increase in the formation of autophagosomes and autolysosomes following the exposure of HT22 cells to 400, 800, and 1,600 μM MA (Figure 7A).

Through qRT-PCR analysis, we observed elevated expression levels of multiple genes involved in the autophagy pathways after MA treatment (Figure 7B). Consistently, Western blot analysis demonstrated a significant increase in the protein levels of PINK1, LC3-II, and p62 in HT22 cells exposed to 800 and 1,600 μM MA (p < 0.05) (Figure 7C). The above results demonstrated that MA exposure leads to mitochondrial damage and triggers autophagy in HT22 cells, playing a crucial role in the neurotoxicity of MA.

RNA epigenetic modifications, particularly m6A, have been reported to play an essential role in the regulation of gene expression. To investigate whether m6A modification is involved in neurotoxicity following MA treatment, we performed MeRIP-seq on both MA-treated and untreated control cells. After stringent filtration, a total of 9,137 m6A peaks containing transcripts from 4,712 genes were identified. Specifically, 6,484 peaks corresponding to 3,665 genes were detected in control groups, and 8,323 peaks were identified after MA treatment, representing transcripts of 4,383 genes (Figures 8A, B).

Moreover, transcripts with m6A modifications were categorized according to the number of modification peaks present in each transcript. Statistical results revealed that each group had over 3,000 genes containing 1–2 m6A modification peaks, with a relatively small number of transcripts containing more than 4 m6A modification peaks (Figure 8C). Additionally, a chi-square test was conducted to examine the composition ratio of the number of m6A-modified peaks between the control and MA groups. The results indicated a significant difference in the composition ratio of the number of m6A-modified peaks between the two groups (Figure 8D). The conserved m6A sequence of RRACH (R represents purine and H represents a non-guanine base) was determined using STREME analysis, revealing the presence of classical m6A motifs in both control and MA-treated groups (Figure 8E). To determine the distribution pattern of m6A modification across the entire transcriptome, the density of m6A modification peaks in both the control and MA-treated groups was examined. The result suggested that m6A modification was enriched in the 5′ untranslated region (5′UTR), coding sequence (CDS), and 3′UTR (Figure 8F). Notably, the density of the m6A modification increased significantly in the 5′UTR and CDS regions, reaching its peak at the stop codon before rapidly decreasing in the 3′UTR.

Subsequently, we examined the differentially m6A-modified genes after MA treatment. Specifically, 616 genes exhibited increased m6A modifications, while 153 genes displayed weakened m6A modifications (Figure 8G). The KEGG pathway analysis revealed that these differentially modified genes were predominantly enriched in the cytochrome P450, Notch signaling pathway, and lysosome (Figure 8H). The most enriched BP terms were associated with the regulation of autophagy, brain development, cell cycle, and negative regulation of cell migration. CC terms were enriched in the endoplasmic reticulum, axon, mitochondrion, and autophagosome. MF terms were found to be involved in ATP binding, metal ion binding, and ubiquitin–protein transferase activity (Figure 8I). These annotation results were consistent with our transcriptomic and proteomic data, indicating that MA induces neurotoxicity by regulating m6A modification of genes involved in autophagy and cell proliferation and migration signaling pathways.

To further investigate the genes whose expression changes were regulated by m6A modification after MA treatment, we analyzed the overlapped genes with differential m6A modification and gene expression between control and MA-treated cells through MeRIP-seq and mRNA-seq data. The results showed that 92 genes and 17 genes were upregulated and downregulated after MA treatment, respectively, accompanied by enhanced m6A modification (Figure 9A). Additionally, 10 genes were upregulated and 13 genes were downregulated, accompanied by a decrease in m6A modification. The level of m6A modification on mRNA transcripts was visualized using IGV. Notably, autophagy-related genes Tsc1, Ctsb, and Atg2a exhibited elevated m6A methylation after MA treatment, suggesting that the activation of autophagy genes might be caused by enhanced m6A modification (Figure 9B). In order to investigate m6A conservation and disease association, we performed ConsRM and RMDisease analyses, respectively. The ConsRM analysis result revealed that 99.5% of m6A-modified regions were non-conservative in differentially modified m6A and DEGs, suggesting that the dysregulation of non-conserved m6A sites was more likely to be associated with disease pathogenesis (Figure 9C). Furthermore, RMDisease analysis showed that 8.9% of m6A genes and 8.7% of m6A peaks were associated with diseases, suggesting a potential correlation between MA and diseases (Figure 9D).

To identify potential regulators of m6A methylation in MA-based neurotoxicity, we analyzed the expression levels of 22 previously reported m6A methylation writers, readers, and erasers. The results revealed that 6 m6A methylation regulators were significantly differentially expressed (p < 0.05), as shown in Table 2. Specifically, after MA treatment, writers VIRMA and METTL14 and reader YTHDF3 were upregulated, while writer RBM15B, reader IGF2BP2, and eraser ALKBH5 were downregulated. Meanwhile, MA treatment induced alterations in both gene expression levels and m6A methylation levels of genes related to the lysosome (Fuca2, Galc, Mcoln1, Naglu, Slc11a2, and Wdr7), autophagy (C9orf72, Irs1, and Tsc1), and adherens junction (Tcf7l2) (Figure 9A).

The CLIP-seq datasets were used to identify potential substrates of m6A methylation regulators, and the results suggested that METTL14 and ALKBH5 were broadly involved in regulating most genes in these four terms (Figure 10A). The readers YTHDF3 and IGF2BP2 specifically recognized m6A methylation on Naglu and C9orf72, respectively. To further explore the regulatory relationships between differentially expressed m6A regulators and these target genes, we used STRING to construct a PPI network. The top five genes, namely, Irs1, Fuca2, Mcoln1, Naglu, and Wdr7, exhibited degrees of 4, 3, 3, 3, and 3, respectively, suggesting the network regulation of m6A modification during neuronal damage caused by MA (Figure 10B).

In order to explore the potential targets of MA, we assessed the molecular interactions between MA and differentially expressed m6A modification regulators, including ALKBH5, IGF2BP2, METTL14, and YTHDF3. All total scores were above 5, indicating stabilized interactions between MA and these m6A modification regulators (Table 3). Theoretically, MA bids ALKBH5 through the formation of a hydrogen bond at Cys201 and Arg399, as well as 12 hydrophobic contacts with Pro271, Phe355, Val203, lle238, Leu14, lle340, Pro129, Arg131, Ala128, Leu341, Asp126, and Arg127 (Figure 11A). MA binds IGF2BP2 through the formation of a hydrogen bond at Gln113, Lys490, and Gly483, as well as nine hydrophobic contacts with Val108, Asp99, Trp95, Val111, Phe558, Arg484, Gly487, Phe559, and Glu112 (Figure 11B). MA bound to METTL14 through the formation of a hydrogen bond at Glu317, as well as 11 hydrophobic contacts with Glu317, Pro319, Arg11, Glu320, Asn323, Arg15, Val328, Lys326, Leu139, Pro327, and Thr316 (Figure 11C). MA bids YTHDF3 through the formation of two hydrogen bonds at Gln351 and one hydrogen bond at Leu350 and Asn468, as well as seven hydrophobic contacts with Ser425, Asn352, Gln494, Tyr424, Gln349, Arg353, and Gly469 (Figure 11D). These results indicated that ALKBH5, IGF2BP2, METTL14, and YTHDF3 were potential targets of MA.

The substrate prediction results showed that the m6A regulator METTL14 was the substrate of most genes (Figure 10A), and the molecular docking results demonstrated a stable interaction between MA and METT14 (Table 3). Therefore, based on these findings, we conducted molecular dynamics research on METTL14. In order to evaluate the binding stability and conformational variability of the MA-METTL14 complex during the MD simulation, we performed RMSD, Rg, and RMSF analyses. The average RMSD values of the complex ranged from 0.5 to 2.7 Å, while the Rg values exhibited fluctuations within 1 Å (Figures 11E, F). The RMSF analysis quantifies residue-specific fluctuations. During the simulation, none of the amino acid residues displayed conformational variations exceeding 3.2 Å relative to their average structure, indicating remarkable stability among backbone atoms in this complex system (Figure 11G). In the 456 amino acid residues, aspartic acid at position 136 demonstrated the highest hydrogen bond occupancy at 36% (Figure 11H). During the MD simulation, the binding free energy of the MA-METTL14 complex gradually decreased, suggesting a more stable binding state between MA and the METTL14 protein (Figure 11I). Molecular docking and MD simulations suggested that METTL14 could be a potential target of MA, which we then verified using the cellular thermal shift assay (CETSA). After treatment with 800 μM MA, the degradation rate of METTL14 in HT22 cells exhibited a deceleration trend with increasing temperature, suggesting that MA treatment enhanced the thermal stability of METTL14 (Figure 11J). The CETSA result suggested that MA binds to METTL14, in accordance with the predictions made by molecular docking and MD simulations.