Epitranscriptomic Analysis of N6-methyladenosine in Infant Rhesus Macaques after Multiple Sevoflurane Anesthesia

Clinical investigations to date have proposed the possibility that exposure to anesthetics is associated with neurodevelopmental deficits. Sevoflurane is the most commonly used general anesthetic in pediatric patients. Animal studies have demonstrated that multiple exposures to sevoflurane during the postnatal period resulted in neuropathological brain changes and long-term cognitive deficits.However, the underlying mechanisms remain to be clarified. In this study, methylated RNA immunoprecipitation sequencing (MeRIP-Seq) was performed to acquire genome-wide profiling of RNA N6-methyladenosine (m6A) in the prefrontal cortex of infant rhesus macaques. The macaques in the sevoflurane group had more m6A peaks than the macaques in the control group (p≤0.05). After sevoflurane treatment, the mRNA levels of YT521-B homology domain family 1 (YTHDF1) and YT521-B homology domain family 3 (YTHDF3) were decreased, and sevoflurane anesthesia dynamically regulated RNA m6A methylation. Gene ontology (GO) analysis revealed that after sevoflurane exposure, genes with increased methylation of m6A sites were enriched in some physiological processes relevant to neurodevelopment, mainly focused on synaptic plasticity. The female macaques had 18 hypermethylated genes. The males had 35 hypermethylated genes, and some physiological processes related to the regulation of synaptic structure were enriched. Rhesus macaques are genetically closer to human beings. Our findings can help in the study of the mechanism of sevoflurane-relevant neurodevelopmental deficits at the posttranscriptional level and can provide new insights into potential clinical preventions and interventions for the neurotoxicity of neonatal anesthesia exposure.


INTRODUCTION
In recent decades, the safety of anesthetic exposure in the developing brain during infancy and childhood has aroused significant public concern and controversy. A sibling-matched cohort study revealed that a single anesthesia exposure may not cause significant reductions in the intelligence quotient between siblings . Additionally, an hour of general anesthesia in early infancy did not alter neurodevelopmental outcomes (Davidson et al., 2016;McCann et al., 2019). However, increasingly more clinical investigations to date have proposed the possibility that exposure to surgery and anesthesia, especially multiple exposures, is associ-ated with neurodevelopmental deficits, including an increased risk of a language disability (Ing et al., 2017), cognitive impairment (Hu et al., 2017), neurobehavioral abnormality, and poorer performance in neuropsychological tests. (Warner et al., 2018;Zaccariello et al., 2019). In December 2016, the US Food and Drug Administration (FDA) issued a safety communication stating that the use of general anesthetic drugs in children less than 3 years old may negatively affect the neurodevelopment of children, with a specific focus on prolonged (defined as an anesthesia duration of > 3 h) or repeated exposures to general anesthesia. Subsequently, warnings were added to the labels for these medicines.
Sevoflurane, one of the most commonly used anesthetic agents in children, has been widely used to uncover the underlying mechanisms of anesthesiarelated neurotoxicity in animal models. Sevoflurane is thought to exert its function by activating inhibitory central nervous system (CNS) receptors, such as c-aminobutyric acid (GABA) receptors type A and glycine; as well as inhibiting N-methyl-d-aspartate (NMDA) receptors, a type of CNS excitatory receptor (Nishikawa and Harrison, 2003;Brosnan and 2012). Nonhuman primates, which are 98% genetically similar to humans, can more closely mimic the pediatric population and be highly translational to explore the effects of neonatal sevoflurane exposure . Researchers found that in developing monkey brains, a long duration of sevoflurane exposure (!8h) at clinical concentrations led to increased neuronal cell death and neuroinflammation (Liu et al., 2015a;Liu et al., 2020). In addition, as the forementioned safety statement of the FDA on repeated anesthesia warned, multiple sevoflurane exposures (mainly 3 times) for infant rhesus macaques could result in long-term neurobehavioral abnormalities, including visual recognition memory impairment (Alvarado et al., 2017) and increased anxiety-related behaviors (Raper et al., 2015;Raper et al., 2018). In our previous work, likewise, we found that the myelinationrelated gene was significantly downregulated in multiple sevoflurane-exposed macaques, which may be associated with myelination deficits and later fine motor impairment (Zhang et al., 2019). Thus, to explore the possible mechanism of anesthesia-induced neurotoxicity, assessing multiple sevoflurane exposures in NHP models are a reasonable approach.
To elucidate neonatal anesthesia-induced neurodevelopmental abnormalities and the mechanisms involved, many works have found that sevoflurane anesthesia leads to dramatic changes in gene expression in different brain regions, and whether posttranscriptional regulation of messenger RNA levels occurs upon sevoflurane anesthesia is unknown. m 6 A in mRNA is the most abundant reversible chemical modification and can mediate its function in a dynamic manner in response to various cellular and environmental signals (Jia et al., 2013;Chang et al., 2017). The reversible establishment of m 6 A modifications and the execution of their biological functions rely on a series of processes involving installation by methyltransferases (writers), removal by demethylases (erasers) and recognition by m 6 A binding proteins (readers). Studies have indicated that m 6 A exerts a critical function in nervous system development (Ma et al., 2018;Shafik et al., 2021), specifically in neurogenesis, axon regeneration, synaptic plasticity, learning and memory (Engel et al., 2018;Shi et al., 2018;Wang et al., 2018;Weng et al., 2018). The essential role of m 6 A in sevofluraneinduced neurodevelopmental deficits is largely unknown.
Here, we investigated the impacts of sevoflurane anesthesia on the dynamic regulation of m 6 A in the prefrontal cortex of young rhesus macaques. Furthermore, this study reveals that RNA m 6 A methylation plays an important role in the regulation of physiological processes pertinent to neurodevelopment after sevoflurane exposure.

Anesthesia of rhesus macaques
The use of rhesus macaques in research at the Institute of Laboratory Animal Science was approved by the Institutional Animal Care and Use Committee (NO. XC17001). Efforts were made to minimize the number of animals in the study. The rhesus macaques received sevoflurane anesthesia as described in our previous studies (Zhang et al., 2019;Cheng et al., 2020). Briefly, for the induction (2-4 min) of general anesthesia, the rhesus macaques received 6-8% anesthetic sevoflurane with 100% oxygen. After reaching a state of tolerance, the animals were tracheally intubated and mechanically ventilated by receiving 2.5-3% sevoflurane and 100% oxygen for 4 h on postnatal day (P) 7 and then on P21 and P35. The heart rate and saturation of pulse oxygen (SpO2) were monitored, consistent with normal physiology. In addition, in previous work, under the same conditions, the vital parameters of exposed macaques were in the normal range, which excluded disturbances to experimental outcomes. (Raper et al., 2015;Alvarado et al., 2017;Zhang et al., 2019). The temperatures of the rhesus macaques were maintained by placing them in a warm box (37°C). The rhesus macaques in the control group received three maternal separations of the same duration (4 h). The rhesus macaques were sacrificed via decapitation under brief sevoflurane anesthesia (3% sevoflurane for 5 min) at the end of the third sevoflurane anesthesia on P35. Immediately afterward, the prefrontal cortex of the rhesus macaques was harvested. There were two male and four female rhesus macaques in the control group, marked as A, B, C, C1, C2, and C3, respectively; and two male and three female rhesus macaques in the sevoflurane group, marked as D, S1, F, S2, and S3, respectively (Fig. 1A). Furthermore, the male macaques in the control group were noted as Ctrl-M, and those in the sevoflurane group were noted as SEV-M. The females were named in the same way as Ctrl-F and SEV-F, respectively.

Quantitative real-time PCR
The total RNA of the prefrontal cortex of the rhesus macaques after treatment was harvested with TRIzol (Invitrogen, USA) according to the manufacturer's instructions. Total RNA was reverse transcribed into cDNA using the Prime-Script RT Reagent Kit (Takara, Japan). The standard SYBR-Green method was used to detect the relative expression of FTO, ALKBH5, METTL3, METTL14, YTHDF1, YTHDF2, YTHDF3 and YTHDC1 with an ABI7500 Real-Time PCR System (Applied Biosystems, USA), and target gene expression was normalized to the level of GAPDH. The relative expression of target genes is presented as 2 ÀDDCt . Additionally, the relative expression of genes with significantly different methylation levels, including FKBP2, ELAVL3, SAFB2, UBALD2, ELFN2, NCOR2, PDE2A, STXBP1, NUP155, TFCP2, and EEF1AKNMT, between groups (Ctrl vs. SEV, Ctrl-M vs. SEV-M, and Ctrl-F vs. SEV-F) was evaluated as described above. Furthermore, in the PCR test, a new sample was added to Ctrl-M and SEV-M, respectively. The primer sequences are listed in Table S1.

RNA extraction and m 6 A MeRIP
Total RNA was extracted using TRIzol reagent (Invitrogen), and the concentration and integrity of total RNA were measured by a Qubit RNA HS Assay and an Agilent 2100 Bioanalyzer (Agilent Technology), respectively. For the MeRIP experiment, approximately 20 mg of total RNA from each sample was fragmented using the 10X RNA Fragmentation Buffer (Invitrogen) by incubating in a preheated thermal cycler for 10 min at 70°C. Fragmented RNA was pelleted by ethanol precipitation. Protein A and protein G magnetic beads were washed twice with IP buffer (150 mM NaCl, 10 mM Tris-HCl pH 7.5, 0.1% IGEPAL CA-630 in nuclease-free H2O) before incubation with 5 mg m 6 A antibody (Synaptic Systems) at 4°C for 2 h. After two washes with IP buffer, the antibody-bead complexes were resuspended in 500 ml of the IP reaction mixture including fragmented total RNA and incubated for 4 h at 4°C. The m 6 A RNA immunoprecipitated with protein A/ G magnetic beads was then washed 3 times using IP buffer for 10 min each at 4°C. Finally, the bead complexes were resuspended in 100 ml of m 6 A competitive elution buffer with continuous shaking for 1 h at 4°C. The supernatant containing the eluted m 6 A RNA was collected in a new tube and purified using phenol:chloroform:isoamyl alcohol (125:24:1).

Library construction and sequencing
The MeRIP libraries using eluted RNA were constructed using the SMARTer Stranded Total RNA-Seq Kit version 2 (Takara/Clontech) according to the manufacturer's protocol. Briefly, the eluted m 6 A RNA and input RNA were directly subjected to first-strand cDNA synthesis without fragmentation. Then, all following steps were based on the manufacturer of the SMARTer Stranded Total RNA-Seq Kit version 2. Libraries for IP RNA were PCR amplified for less than 16 cycles, and input libraries were amplified for less than 12 cycles. All libraries were analyzed by an Agilent 2100 Bioanalyzer (Agilent Technologies) and quantified by real-time PCR before sequencing. Then, different libraries were pooled according to the requirements of the effective concentration and target off-board data quantity and then sequenced on the NovaSeq platform. The sequencing strategy was PE150.

Data analysis
In brief, low-quality bases were trimmed by using Trim Galore (version: 0.5.0). Adapter sequences, which are Illumina TruSeq adapter sequences, were trimmed by using cutadapt implemented in Trim Galore. Clean reads were mapped to the genome (rheMac8) using hisat2 (version: 2.1.0) with default parameters. The R package exomePeak (version: 2.1.2) was used to call peaks and detect differential peaks from filtered alignment files with parameters that ''WINDOW_ WIDTH=200 SLIDING_STEP=30 FRAGMENT_ LENGTH=150 DIFF_PEAK_ABS_FOLD_CHANGE=2 FOLD_ENRICHMENT=2".
Only peaks with FDR < 0.05 and fold change >=2 were identified as significantly differential peaks. Gene expression calculation was performed using the StringTie software using default parameters. Gene expression profiling was based on the number of reads. TPM (transcripts per million mapped reads) values were used to estimate the expressed values and transcript levels. DESeq2 was selected to identify differentially expressed genes (DEGs). Genes with an adjusted P value (padj) < 0.05 and abs (log2(fold change)) > 1 were characterized as DEGs.
Gene ontology (GO) and pathway enrichment analyses were performed using the Database for QRT-PCR indicated that FKBP2 and UBALD2 mRNA expression increased in the macaques prefrontal cortex after multiple sevoflurane exposure (in the Ctrl group, n = 6, in the SEV group, n = 5; F = 1.069, p = 0.0325; F = 2.286, p = 0.0150, respectively, one-way ANOVA). Relative mRNA level of ELAVL3 and SAFB2 didn't change (in the Ctrl group, n = 6, in the SEV group, n = 5; F = 3.527, p = 0.6226; F = 1.081, p = 0.9574, respectively, one-way ANOVA). (E) Volcano plot of differentially expressed genes (DEGs). DEGs with log2FoldChange (log2FC) > 1 were shown in red; DEGs with log2FoldChange (log2FC) < -1 were in green (P < 0.05). (F) Unsupervised hierarchical cluster image shows the differential gene expression profiles in two groups (SEV vs. Ctrl). The Heat maps display a color scale: red indicates upregulation, whereas blue represents downregulation. (SEV: sevoflurane; Ctrl: Control). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Annotation, Visualization and Integrated Discovery. The ontology consists of three parts: the molecular function (MF), cellular component (CC) and biological process (BP). The p value indicates the significance of the GO term enrichment of the genes. Pathway enrichment analysis is a functional analysis that maps genes to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The p value denotes the significance of the pathway correlated with the conditions.

General features of m 6 A methylation in rhesus macaques
To explore whether and how the transcriptome-wide profiles of m 6 A marks change after sevoflurane treatment, we used the exomePeak R package software to identify the significant m 6 A peaks characteristic of different groups (FDR < 0.05, fold_enrichment ! 2). In the control group, there were 601 m 6 A peaks in sample A, 535 in sample B, 483 in sample C, 452 in sample C1, 703 in sample C2, and 1210 in sample C3. Furthermore, in the SEV group, there were 1146 m 6 A peaks in sample D, 1178 in sample F, 852 in sample S1, 1673 in sample S2, and 1236 in sample S3. After sevoflurane exposure, we noted a dramatic increase in the number of m 6 A marks in the SEV group compared to the control group (p 0.05, Fig. 1B). We performed a de novo motif analysis on the detected m 6 A sites and found them to be enriched in the conserved RRACH motif (R is A/G and H is A/C/U) that was reported previously in other species (Luo et al., 2014;Haussmann et al., 2016;Wang et al., 2021), proving the quality of our data (Fig. 1C). In the examination of the m 6 A distribution, peaks were categorized into three transcript segments: the 5 0 UTR, the coding sequence (CDS), and 3 0 UTR. We found in all of the samples that the m 6 A peak positions in transcripts were most often harbored in coding sequence (CDS) regions (Fig. 1D). Furthermore, for most of the samples, the coding sequence segment was the most highly enriched segment type with m 6 A after each segment was normalized. We displayed the representative enrichment distribution of sample C1 (Fig. 1E). To confirm whether the m 6 A eraser, writer, and reader systems may be altered after sevoflurane treatment, we examined the mRNA levels of FTO, ALKBH5, METTL3, METTL14, YTHDF1, YTHDF2, YTHDF3, and YTHDC1 in the prefrontal cortex of macaques after multiple sevoflurane anesthesia by qPCR.
We found that among numerous m 6 A-related proteins, only the mRNA levels of YTHDF1 and YTHDF3 were decreased after sevoflurane treatment (Fig. 1F). These data suggest that sevoflurane-induced m 6 A regulation likely occurs at the m 6 A reader level instead of the writer or eraser level.

Altered m 6 A methylation and DEGs between groups (Ctrl vs. SEV)
To investigate the potential regulatory roles of m 6 A marks after sevoflurane exposure, we examined the differences in the m 6 A-containing transcripts and the differences in the methylation levels of those transcripts compared with the control group. A volcano plot showed that there were 22 hypermethylated genes and 1 hypomethylated gene in the SEV group compared with the control group ( Fig. 2A). Some well-known genes, such as RAB13, which is related to the regulation of neurite outgrowth (Sakane et al., 2010), and RIN1, which is related to the modulation of neuronal plasticity (Szı´ber et al., 2017), were found to be highly methylated ( Fig. 2A).
To uncover the functions of m 6 A in sevofluranetreated rhesus macaques, protein-coding genes containing differentially methylated m 6 A sites compared with the control group were selected for GO enrichment analysis and KEGG pathway analysis. GO analysis revealed that genes with hypermethylated m 6 A sites were enriched in the intracellular organelle part, intracellular part, RNA binding, organic cyclic compound binding, regulation of mRNA processing, mRNA metabolic process, response to endogenous stimulus, response to hormone, and regulation of synaptic plasticity (Fig. 2B). This finding is of particular interest because synaptic plasticity is known to form the basis for learning and memory and has been shown to be involved in cortical development (Greenhill et al., 2015;O'Reilly et al., 2019). Furthermore, genes with hypomethylated m 6 A sites were found to be involved in the TGF-beta signaling pathway, cell cycle, and cellular senescence (Fig. 2C).
To explore the relative expression of genes with hypermethylated m 6 A peaks, we examined the mRNA levels of FKBP2, ELAVL3, SAFB2, and UBALD2, which were among the top 10 hypermethylated genes between groups, in the prefrontal cortex of rhesus macaques after multiple sevoflurane anesthesia by qPCR, as depicted in Fig. 2A. We found that the mRNA levels of FKBP2 and UBALD2 were increased after sevoflurane treatment. The mRNA levels of ELAVL3 and SAFB2 were not altered after sevoflurane treatment (Fig. 2D). DEG analysis revealed 158 upregulated and 166 downregulated genes in the SEV group compared with the control group (Fig. 2E). The top 30 differentially expressed genes are shown in Table S2. Furthermore, a significantly different gene expression pattern was discovered after sevoflurane treatment, and all samples could be easily distinguished through the DEG cluster diagram (Fig. 2F).

Altered m 6 A methylation and DEGs of different genders after sevoflurane anesthesia
To evaluate the effect of sevoflurane treatment on macaques of different genders, we analyzed the possible functions of m 6 A marks in males and females, respectively, compared with the control group. The volcano plot showed that males had 35 hypermethylated genes and 1 hypomethylated gene in the SEV group compared with the control group (Fig. 3B), while females had 18 hypermethylated genes and 2 hypomethylated genes (Fig. 4A). For the CC category, GO enrichment analysis revealed the enrichment of genes with hypermethylated m 6 A sites in sevofluraneinduced male macaques for synaptic development, such as the synapse, neuron part, and synaptic membrane (Fig. 3A). In females, we observed that m 6 A-highly modified genes in the SEV group compared with the control group were enriched in mRNA processing, regulation of mRNA metabolic process, mRNA polyadenylation, and RNA binding (Fig. 4B). Furthermore, KEGG pathway analysis showed enrichment of genes with hypomethylated m 6 A sites in sevoflurane-induced female macaques for the regulation of methylation, including transferase activity, transferring one-carbon groups, and methyltransferase activity (Fig. 4C). As described above, to further identify the relative expression of genes containing differentially methylated m 6 A sites, we examined the mRNA levels of ELFN2, NCOR2, PDE2A, and STXBP1, which were among the top 10 known hypermethylated genes between groups, in the prefrontal cortex of male macaques after multiple sevoflurane anesthesia by qPCR and discovered that only the mRNA level of STXBP1 was decreased after sevoflurane treatment. The mRNA levels of ELFN2, NCOR2, and PDE2A were not changed after sevoflurane exposure (Fig. 3C).
Furthermore, in females, the mRNA level of the hypomethylated gene EEF1AKNMT was increased after sevoflurane anesthesia. The mRNA levels of NUP155, ELAVL3, and TFCP2, which were among the top 10 hypermethylated genes, were not changed (Fig. 4D), as depicted in Fig. 4A. Next, we performed DEG analysis. In particular, in males, there were 311 upregulated and 384 downregulated genes in the SEV group compared with the control group (Fig. 3D); and in females, 80 upregulated and 207 downregulated genes were observed (Fig. 4E). Furthermore, DEG cluster analysis revealed different gene expression patterns for both males and females (Figs. 3E, 4F).

DISCUSSION
N6-Methyladenosine (m 6 A) has been regarded as a pivotal epigenetic marker that regulates various aspects of RNA function, including RNA stability, translation efficiency, and RNA-protein interactions (Wang et al., 2014;Liu et al., 2015b;Wang et al., 2015). In the nervous system, m 6 A is highly abundant and sensitive to neural activity, which plays a critical role in various developmental and physiological processes, such as neurogenesis, axon regeneration, oligodendrocyte maturation, synaptic transmission, circadian rhythm, cognitive function and stress response (Yoon et al., 2017;Anders et al., 2018;Fustin et al., 2018;Koranda et al., 2018;Weng et al., 2018;Xu et al., 2020;Widagdo et al., 2021). Recent studies have revealed that m 6 A imbalance can contribute to the pathogenesis of brain cancer, neurodegenerative diseases and neuropsychiatric disorders (Choudhry et al., 2013;Li et al., 2018;Su et al., 2018;Huang et al., 2020). In our study, to investigate the possible impact of m 6 A posttranscriptional regulation on sevofluranerelated neurotoxicity in the developing brain, we revealed changes in the m 6 A methylation profile of the prefrontal cortex of infant rhesus macaques after sevoflurane anesthesia. There were eleven macaques in the research, which are genetically closer to humans than rodents and widely used as models of healthy humans in neuroscience research. (Poirier and Bateson, 2017). Thus, regarding neonatal sevoflurane exposure-related neurodevelopmental effects, the findings of our work are plausible to some extent.
Among a series of m 6 A-related proteins, we found that YTHDF1 decreased after multiple sevoflurane exposures, consistent with our previous work . The m 6 A-binding protein YTHDF1, one of the m 6 A readers, directly interacts with m 6 A sites and promotes the translation efficiency of m 6 A-containing mRNAs (Wang et al., 2015;Shi et al., 2017). YTHDF1 has been found to perform critical functions in the central nervous system, such as expediting new protein synthesis required for long-lasting changes in synapse plasticity and thus facilitating learning and memory . In addition, YTHDF1 can attenuate sevoflurane-induced fine motor and cognitive impairments and decrease the expression of synaptophysin, a presynaptic marker . Therefore, we speculate that the reduction of YTHDF1 has an important role in sevoflurane-induced neurotoxicity. YTHDF3, another m 6 A reader, was also found to be decreased in our study. Interestingly, YTHDF3 is highly expressed in brain tissues compared with YTHDF1 and YTHDF2 (Lein et al., 2007). Previous studies have found that YTHDF3 can promote protein synthesis in synergy with YTHDF1 (Li et al., 2017) and affect methylated mRNA decay mediated through YTHDF2 (Shi et al., 2017). However, its biological functions in the central nervous system are underexplored, and our study found that it was decreased after sevoflurane treatment. Thus, its specific role in sevofluraneassociated neurotoxicity merits further research in our future work.
In the current study, the results of GO analysis of hypermethylated genes between groups were shown to be enriched in many important biological processes, which were mainly related to mRNA metabolism, the regulation of phosphorylation activity, and the regulation of synaptic plasticity or synaptic structure. Increasing evidence suggests that the synaptic structure is highly dynamic and that synaptic plasticity, the cellular basis of learning and memory, plays a key role in functional network alterations in neurodevelopmental and memory disorders (Holtmaat and Svoboda, 2009;Diering and Huganir, 2018;Mansvelder et al., 2019;Sun et al., 2020). Additionally, the impairment of synaptic plasticity is often correlated with cognitive deficits (Salazar et al., 2017;Mayne and Burne, 2019), and a large body of animal studies have demonstrated that neonatal sevoflurane-induced memory impairment co-occurs with the inhibition of synaptic plasticity (Liang et al., 2017;Yu et al., 2018;Fan et al., 2021). However, the mechanism is obscure. Our study revealed that m 6 A exerted an important influence on sevoflurane-induced alteration of synaptic plasticity via posttranscriptional regulation. In male macaques, in addition to synaptic regulation, we also found that the hypermethylated genes in the regulation of phosphorylation activity were enriched after sevoflurane treatment. Sevoflurane can induce tau phosphorylation and cognitive impairment (Tao et al., 2014;Yu et al., 2020). We speculate that the hypermethylation of phosphorylationrelated genes may be involved in the process of tau phosphorylation and therefore influence the neurodevelopment of young rhesus macaques. However, the findings need further investigation.
After sevoflurane anesthesia, we found that the relative expression of the hypermethylated gene FKBP2, which plays a role in immunoregulation by binding the immunosuppressive compound rapamycin that inhibits mTOR signaling, was increased (Hendrickson et al., 1993). This can provide new insights into why rapamycin may act as a potential therapeutic agent for sevofluraneinduced changes in the developing brain , possibly because of the upregulation of FKBP2 and the corresponding increased binding to rapamycin. For male macaques, the relative expression of the hypermethylated gene STXBP1 was decreased after sevoflurane anesthesia. STXBP1 is an essential protein for presynaptic vesicle release (Stamberger et al., 2017), and mutationrelated loss of functions in STXBP1 has been associated with a series of neurodevelopmental disorders because of reduced synaptic transmission (Kovacevic et al., 2018;Lanoue et al., 2019;Lammertse et al., 2020), which can help us understand the underlying mechanism of the effect of sevoflurane on synaptic transmission. In addition to the differential expression of m 6 A highly modified genes, we displayed the top 30 differentially expressed genes in Table S2. After multiple sevoflurane exposures, brain-derived neurotrophic factor (BDNF) was downregulated, which agreed with our previous study . BDNF not only regulates neuronal survival and differentiation but also plays an essential role in modulating synaptic transmission and plasticity. Imbalances in BDNF levels can result in deficits in cognition, learning and memory (Su et al., 2014;Mizui et al., 2016;Yu et al., 2016). Additionally, we observed the downregulation of immediate early genes (IGEs), including NPAS4, EGR1, and ARC, which are critical regulators of synaptic plasticity and long-term memory formation (Plath et al., 2006;Ramamoorthi et al., 2011). Studies have revealed that deficiency in IGE expression could contribute to the impairment of learning ability and memory formation (Sun and Lin, 2016;Zhang et al., 2018). Altogether, the altered expression of these genes, along with the aforementioned results of GO analysis, reveal the possible effect of sevoflurane on synapses, which may shed light on the underlying mechanism of sevoflurane-related learning and memory disorders.
There are several limitations in this study. First, the rhesus macaque experiment had a limited sample size (Ctrl = 6 and SEV = 5). A larger sample size may be required to draw stronger conclusions. Second, based on the analysis of the sequencing data, it is difficult to determine whether the observed alterations of the m 6 A methylation profile in the developing brain are transient or lasting in our current study. But we have preliminarily explored the impacts of multiple sevoflurane exposures on neurodevelopment at a posttranscriptional regulation level and displayed the DEGs, which may provide clues for future studies. Finally, after sevoflurane treatment, we found that some GO terms and DEGs could affect synaptic plasticity. Although we could not verify the corresponding long-term behavioral changes in this research, Baxter et al. applied the same experimental design as ours and found that infant rhesus monkeys exposed to sevoflurane at P7, P21, and P35 showed memory impairment after the first year of life (Alvarado et al., 2017). Thus, our study may elucidate the possible role of m 6 A methylation and DEGs in synaptic plasticity and the later influences on memory.
In summary, sevoflurane anesthesia can dynamically regulate RNA m 6 A methylation in the prefrontal cortex of young rhesus macaques, and some physiological processes relevant to neurodevelopment, mainly focused on synaptic plasticity, were altered after sevoflurane treatment.

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/geo/, GSE174742.

ETHICS STATEMENT
The use of rhesus macaques in research at the Institute of Laboratory Animal Science was approved by the Institutional Animal Care and Use Committee (protocol number #XC17001).