Decreased MEF2A Expression Regulated by Its Enhancer Methylation Inhibits Autophagy and May Play an Important Role in the Progression of Alzheimer’s Disease

Alzheimer’s disease (AD) is a neurodegenerative disease characterized by amyloid plaques and neurofibrillary tangles which significantly affects people’s life quality. Recently, AD has been found to be closely related to autophagy. The aim of this study was to identify autophagy-related genes associated with the pathogenesis of AD from multiple types of microarray and sequencing datasets using bioinformatics methods and to investigate their role in the pathogenesis of AD in order to identify novel strategies to prevent and treat AD. Our results showed that the autophagy-related genes were significantly downregulated in AD and correlated with the pathological progression. Furthermore, enrichment analysis showed that these autophagy-related genes were regulated by the transcription factor myocyte enhancer factor 2A (MEF2A), which had been confirmed using si-MEF2A. Moreover, the single-cell sequencing data suggested that MEF2A was highly expressed in microglia. Methylation microarray analysis showed that the methylation level of the enhancer region of MEF2A in AD was significantly increased. In conclusion, our results suggest that AD related to the increased methylation level of MEF2A enhancer reduces the expression of MEF2A and downregulates the expression of autophagy-related genes which are closely associated with AD pathogenesis, thereby inhibiting autophagy.


INTRODUCTION
which are neurotoxic and cause neuronal loss, synapse reduction, neurological degeneration, and brain atrophy. Pathological studies have demonstrated that NFTs and amyloid plaque deposition initially occurred in the cortical and hippocampal tissues of AD and subsequently spread to the whole brain (Braak and Braak, 1995). Inhibiting neurotoxicity by reducing amyloid plaque deposition and NFTs has been unsuccessfully attempted (Doody et al., 2013;Ostrowitzki et al., 2017).
Autophagy consists of a series of complex physiological processes in cells, which can eliminate misfolded proteins and damaged organelles, promote the synthesis of biofilms and transport of vesicles, and thus play a key role in reshaping the cell structure and regulating energy metabolism, resisting adverse external stimuli and stabilizing cell homeostasis. Autophagy disorders are closely related to AD (Nixon et al., 2005). Inhibiting lysosomal proteolysis produces similar neuropathological manifestations in wild-type mice and exacerbates amyloid plaque deposition and autophagy pathology in mouse models of AD (Nixon and Yang, 2011). Presenilin 1 mutations, associated with familial AD, result in decreased maturation of the lysosomal v-ATPase and, thus, directly increased lysosomal pH and impaired lysosome function, which would be predicted to reduce autophagosome clearance (Lee et al., 2010). Another genetic risk factor for AD is mutations in apolipoprotein E 4 (ApoE4). ApoE4 destabilizes lysosomal membranes in an allele-specific manner. Other factors include reactive oxygen species and amyloid plaque and oxidized lipids and lipoproteins, which also contribute to AD by impeding lysosomal proteolysis, damaging lysosomal membranes, and disrupting lysosomal integrity, thereby releasing proteases that can mediate neuronal cell death (Boya and Kroemer, 2008;Nixon, 2013). Moreover, a study found that autophagy inducers such as rapamycin in 3xTg-AD mice can effectively reduce the deposition of amyloid plaques in the brain and improve cognitive performance (Majumder et al., 2011). However, the relationships between the pathological progression of AD and autophagy-related genes remain unclear, and the mechanisms have not been elucidated.
Systems biology concepts and methods provide multivariate approaches to holistically analyze the larger interactive network of biological pathways and identify important players in AD onset and progression. In this study, multiple datasets were downloaded, and bioinformatics methods were used to analyze the correlation between autophagy-related genes and pathological progression of AD, assessing the reasons for their differential expression, thereby providing new clues for elucidating the mechanism of AD.
GSE80970 (Smith et al., 2018) contained prefrontal cortex and superior temporal gyrus tissue from 147 participants with varying levels of AD pathology. DNA modifications for these samples were quantified using the Illumina Infinium Human 450K Methylation Array.
Single-nucleus RNA-seq (snRNA-seq), GSE138852 (Grubman et al., 2019), on the entorhinal cortex from control and AD brains of 16 participants, yielding a total of 13,214 high-quality nuclei, were used to check the gene expression distributions across cells in Alzheimer's disease brains.
In GSE62420, brains of 4-, 12-, and 22-month-old C57Bl/6J mice were collected and dissected into four regions: cerebellum, cortex, hippocampus, and striatum. Microglia were extracted from each region using a magnetic bead-based approach. Total RNA was immediately isolated for purified microglia and stored (−80 • C) until performing microarray analysis of purified microglia and regional brain homogenates (n = 56).
The autophagy gene list with a total of 530 autophagyrelated genes was derived from the Gene Ontology with the term "autophagy" (GO: 0006914) in the Homo sapiens organism.

Weighted Gene Coexpression Network Analysis
We extracted the autophagy genes to perform WGCNA with expression data retrieved from GSE84422 (Wang et al., 2016) microarray data. The R package "WGCNA" was applied to find clinical trait-related modules and hub genes among them as previously described (Zhang et al., 2014). The adjacency matrix was transformed into topological overlap matrix. According to the topological overlap matrix-based dissimilarity measure, genes were divided into different gene modules. Herein, we set softthresholding power as 6 (scale free R 2 = 0.85), cut height as 0.25, and minimal module size as 10 to identify key modules. The module with the highest correlation with clinical traits (age, sex, race, PIM, pH, CDR, Braak, CERAD, PLQ_Mn, NPrSum, and NTrSum) was selected to explore its biological function through gene ontology (GO) analyses and to screen hub genes. Hub genes were defined as those with gene significance > 0.3 and module membership > 0.8.

Differential Expression Analysis
For the microarray differential expression analyses of GSE84422 (Wang et al., 2016), robust multichip average (RMA) was used for background correction of raw gene expression matrixes, then log 2 transformation of expression matrixes. The "affy" R package was used for quantile normalization and median polish algorithm summarization. Next, all gene probes were mapped into gene symbols by the affymetrix annotation files. The "limma" (linear models for microarray data) R package was performed for identifying DEGs between definite AD samples and healthy controls, and the results were visualized using the volcano plot and heat map. Cutoff criteria for screening DEGs were p < 0.05 and | log2fold change| ≥ 1.3. The screened seven genes were obtained by taking the intersection of the DEGs of GSE84422 and hub genes related to pathological progression of AD from WGCNA. The differential expression analyses were performed in GSE84422, GSE122063 (McKay et al., 2019), GSE118553 (Patel et al., 2019), and GSE132903 (Piras et al., 2019) datasets and log2fold changes in the screened seven genes and significantly differences (AD versus control expression levels) were represented with the corresponding bar plot.

Transcription Factors Enrichment
To shed further light on the functions of the candidate genes, DAVID online tools (Huang et al., 2009;Wishart et al., 2009) in the UCSC database were used for transcription factor annotations. The motif matrix profile MA0052.4 of MEF2A was downloaded from JASPAR Fornes et al. (2020), and Find Individual Motif Occurrences (FIMO, Grant et al., 2011) of motif-based sequence analysis tools (MEME suite 5.3.3, Bailey et al., 2009) was used to scan sequences of candidate genes for individual matches to the motif of MEF2A. The positions and sequences of the screened seven genes were inquired in UCSC, and promoters were defined as the 2,000-bp window centered on the transcript start site of genes.

Correlation Analysis
Correlations between each transcription factor and the expression of downstream genes were analyzed (Pearson's correlation) in GSE84422 (Wang et al., 2016), and result was represented with a heatmap. The correlation between MEF2A and the screened seven genes was validated (Pearson's correlation) in GSE118553 (Patel et al., 2019), and result was represented with scatter plots. MEF2A mRNA expression levels (signal intensity) in different groups in GSE84422 and GSE118553 were shown, and p-values were calculated using GraphPad. In GSE84422, three groups (control, definite AD, possible AD) were classified according to neuropathology category as measured by CERAD (to unify the results, we combined the diagnosis of possible AD and probe AD into possible AD). In GSE118553, participants in the control group were classified as showing no clinical sign of any form of dementia and no neuropathological evidence of neurodegeneration. Participants in the AsymAD group were defined as clinically dementia free at the time of death, but neuropathological assessment at autopsy revealed hallmark AD pathology. Participants in the AD group had both a clinical diagnosis of AD at death and received confirmation of this diagnosis through neuropathological evaluation at autopsy. A one-way analysis of variance (ANOVA) followed by a least significant difference (LSD) test was used for comparison among groups.

snRNA-Seq Data Analysis
Single-nucleus RNA-seq data were downloaded from the website 2 as described by Grubman et al. (2019). We then analyzed the expression of MEF2A in different cell types and different groups. The MEF2A mRNA expression levels (logCounts) in different cell types were represented and significant differences (compared with microglia) were calculated using GraphPad. The expression of MEF2A and screened seven genes (logCounts) in different groups were represented using the ggplot2 package in the R software and significant differences (AD versus control expression levels) were calculated using GraphPad. Following the standardization of GSE62420 (Grabert et al., 2016) microarray data, the signal intensity of MEF2A in microglia and other mixed cell types were shown and significant differences (compared with microglia) were calculated using GraphPad.

Methylation Array Data Analysis
The ChAMP package in R software (Tian et al., 2017) was used to analyze the differentially methylated positions (DMPs), and the screening condition yielded a value of p < 0.05. The DMPs range around MEF2A (from the last gene to the next gene) was approximately 99,676,703-100,882,647 bp of chromosome 15, and the coMET package of the R software (Martin et al., 2015) was used to create a Manhattan plot with a threshold value of p = 0.05. These DMPs were annotated by UCSC, and the means of CpG methylation levels in different groups were represented in a violin plot in the enhancers.

Super Enhancers and Chromatin Interaction Analysis
The super enhancer (SE) analyses, as Nott et al. (2019) described, were performed in UCSC. The ATAC-seq, H3K27ac and H3K4me3 ChIP-seq, and PLAC-seq in all types of cells in the 99,950,000-100,200,000-bp section of chromosome 15 were queried, and the original images were downloaded.

Cell Culture
The mouse microglia cells (BV2), obtained from the Cell Resource Center, Peking Union Medical College (China), were cultured in Dulbecco's modified Eagle's medium (DMEM, Gibco, United States) with 10% fetal bovine serum (Gibco, United States) in 5% CO 2 at 37 • C. For in vitro transfection (n = 3), the target and control siRNA (GenePharma Co., Ltd., China) were transfected into BV2 cells using Lipofectamine 2000 (Invitrogen, United States) according to the manufacturer's guidelines. Cells were collected 48 h after transfection. The siRNA sequences were listed ( Table 2).

Quantitative Real-Time PCR
Total RNA was extracted using TRIzol R Reagent (Life Technologies, Grand Island, NY, United States) and reverse transcribed into cDNA using PrimeScript TM II 1 st Strand cDNA Synthesis Kit for qPCR (TaKaRa, Tokyo, Japan) according to the manufacturer's instructions. Quantitative real-time PCR was performed to detect the gene mRNA levels using 2 × Universal SYBR Green Fast qPCR Mix (ABclonal, Wuhan, China). Primers were synthesized by Sangon (Sangon Biotech Co., Ltd., Shanghai, China). The qPCR conditions were as follows: 95 • C for 3 min and 40 cycles of 95 • C for 5 s and 60 • C for 30 s. Melting curves were tested to assess the accuracy of the PCR analysis. The 2 − Ct was calculated to analyze the gene expression levels: Ct = Ct (target gene) -Ct (β-actin gene), Ct = Ct (treatment) -Ct (control). The primer sequences were listed (Table 3).

Statistical Analysis Software
RStudio (version 3.6.2) was used for data processing, GraphPad (version 8.0) was used for the calculation of significance and means of the differences, and Adobe Illustrator 2020 (version 24.0.1) was used for image processing.

WGCNA and Differential Expression Analysis Were Used to Determine Hub Genes Related to AD Clinical Phenotypes
To assess the relationship between gene expression and clinical phenotypes of AD, dataset GSE84422 (Wang et al., 2016) were downloaded and the WGCNA package, which can cluster genes and divide them into different modules and associate the genes of each module with the clinical phenotypes, was used to cluster autophagy genes and divide them into 11 hub gene modules with different module colors, according to the gene expression correlation patterns ( Figure 2B). The correlation between each module and clinical phenotypes was calculated (Figure 2A).
To investigate whether the expression levels of autophagyrelated genes were different between healthy controls and AD, we used the limma package to analyze differential expression of autophagy-related genes in definite AD and healthy controls in dataset GSE84422, and then screened 71 DEGs shown in a heatmap ( Figure 2C). Compared with healthy controls, the expression of autophagy genes in AD was mostly reduced (69 genes) and only two genes were increased ( Figure 2D). The autophagy-related genes of AD were inhibited, suggesting that their autophagy function was lower than that of healthy controls.
Seven genes were screened by taking the intersection of the autophagy-related genes with significant differences and the black module hub genes with strong correlation with pathological progression of AD screened by WGCNA (Table 4).
Validating the multiple differences of the screened seven genes (Figure 2E), we found that compared with the control group, all screened seven genes showed a significant decrease in the AD group in GSE118553 (Patel et al., 2019; p < 0.05); six genes showed a significant decrease in the AD group in GSE132903 (Piras et al., 2019; p < 0.05); four genes showed a significant decrease in the AD group in GSE122063 (McKay et al., 2019; p < 0.05).
FIGURE 2 | WGCNA coexpression analysis and differential expression analysis. GSE84422 (Wang et al., 2016) was used to perform WGCNA. (A) The correlation between each module and clinical phenotypes is shown, including 11 modules and 11 clinical phenotypes. In each unit, the numbers above show the correlation, and the numbers below show the p-value. Abbreviations: CDR, clinical dementia rating; Braak, Braak NFT score; CERAD, CERAD diagnoses and ratings of pathology; PLQ_Mn, average of neuritic plaque counts in five cardinal cortical regions; NPrSum, sum of CERAD semiquantitative rating scores for all cortical regions examined neuropathologically; NTrSum, sum of semiquantitative NFT density ratings for all cortical regions examined. (B) Cluster dendrogram of the coexpression network modules was produced based on the autophagy genes, including 11 modules. (C) DEGs in healthy control and AD in GSE84422 were shown in the heatmap. Cut-off criteria for screening DEGs were p < 0.05 and | log2fold change| ≥ 1.3. (D) Autophagy genes with significantly different expression in GSE84422 were shown in the volcano plot. Red spots indicate upregulated genes, and blue spots indicate downregulated genes. Cut-off criteria for screening DEGs were p < 0.05 and | log2fold change| ≥ 1.3. (E) Validation of the gene expression levels of BNIP3, CDK5R1, HERC1, ITPR1, OPTN, UBQLN2, and USP33 between healthy control and AD was shown in three other datasets, GSE122063 (McKay et al., 2019), GSE118553 (Patel et al., 2019), and GSE132903 (Piras et al., 2019). Y-axes in the left indicate log2fold change of the screened seven genes in each dataset. Significant differences (AD versus control expression levels) were performed, *p < 0.05 compared with the control group.

MEF2A Expression Is Related to the Screened Seven Genes and AD Neuropathological Category
The online annotation website DAVID was used to perform transcription factor analysis for 71 DEGs and 16 genes in the black module (total 80 genes) in the UCSC database. The result showed that 18 transcription factors were identified by enrichment analysis of 80 genes (p < 0.05), including 73 genes enriched in organic cation transporter 1 (OCT1), 66 genes enriched in Ecotropic Virus Integration Site 1 (EVI1), and 65 genes enriched in MEF2 (Supplementary Figure 1). This data suggested that these 18 transcription factors were likely to play important roles in the regulation of autophagy genes in AD. To further investigate the autophagy genes with differences and related to pathological progression of AD, the screened seven genes were annotated by DAVID, revealing that the screened seven genes were regulated by transcription factors MEF2A and CUX1 ( Figure 3A). Transcription factors can specifically bind to target motifs to regulate the expression of downstream genes; therefore, the expression of transcription factors was consistent with the expression of downstream genes. We analyzed the correlation between transcription factors and the expression of the screened seven genes obtained from the query of the UCSC database in GSE84422 (Wang et al., 2016; Figure 3C). Among them, the expression of MEF2A was most correlated to the screened seven genes. The motifs of the screened seven genes matched to MEF2A (p < 0.001) were scanned by MEME tools (described in Methods section "Transcription Factors Enrichment") and were shown in UCSC ( Figure 3B). The correlation between transcription factors and the expression of the screened seven genes were validated in GSE118553 (Patel et al., 2019 ; Supplementary Figure 2). These results showed that MEF2A had the strongest correlation with the screened seven genes, so that we selected MEF2A for further study. We downloaded dataset GSE84422 and extracted the mRNA expression data of MEF2A, and then compared different neuropathological categories of the sample ( Figure 3D) and found no significant difference in the possible AD group compared with the control group; in addition, the definite AD group showed a significant difference (p < 0.05). The expression of MEF2A in the control group was the highest, followed by the possible AD group, and finally the definite AD group. To validate our results, another dataset GSE118553 was used ( Figure 3E). The MEF2A expression levels were significantly lower in the AD group than in the control group (p < 0.05) in different brain regions (entorhinal cortex, frontal cortex, and temporal cortex), except for the cerebellum, while AsymAD revealed no significant difference (p > 0.05) compared with control.
In our experiment, to verify the relationship between MEF2A and the screened seven genes, three siRNA were transfected into BV2 cells to knock down the expression of MEF2A, and the mRNA levels of the seven autophagy-related genes were detected using qPCR (Figure 4). These results showed that most of the screened seven genes were significantly decreased following siRNA transfection (p < 0.05), suggesting that the expression of the screened seven genes was influenced by MEF2A.

MEF2A Expression Is Cell Type Specific and Mainly Concentrated in Microglia
The brain contains different cell types, which are responsible for different physiological processes. Therefore, transcriptome sequencing of different subgroups of cells can reflect the functions of different types of cells. Grubman et al. (2019) performed singlecell sequencing on the entorhinal cortex of participants with AD and healthy controls, and the cells with different types of markers were clearly divided into eight subgroups ( Figure 5A). The expression levels of MEF2A detected in microglia were significantly higher than those in other subgroups ( Figure 5B). The single-cell sequencing dataset, GSE138852 (Grubman et al., 2019), showed that the expression of MEF2A in microglia was significantly (p < 0.05) higher than that in other subgroups (Figure 5C), validating the above results. Furthermore, the expression levels of MEF2A in AD were significantly (p < 0.05) higher than that in healthy controls in microglia, doublet, oligodendrocyte, and oligodendrocyte progenitor cells (OPC, Figure 5E). Besides, there were no significant difference between AD and healthy controls of the screened seven genes in microglia (p > 0.05, Supplementary Figure 3).
In another dataset, GSE62420 (Grabert et al., 2016), microglia from mouse brain tissue were purified according to specific markers, and gene expression was measured. We downloaded the expression matrix and screened the expression levels of MEF2A in each tissue and found that MEF2A was significantly (p < 0.05) upregulated in microglia compared with mixed brain cells across different brain regions ( Figure 5D). These results suggested that the expression of MEF2A in brain tissue was cell type specific, and the expression of MEF2A in microglia was significantly higher than that in other cells.

Enhancer Region Methylation Regulates MEF2A Expression
We then assessed whether epigenetic regulation could alter the expression of MEF2A because both the pathological progression of AD and tissue specificity could change it. SEs are important gene control elements composed of a series of enhancers. We analyzed datasets involved in the interactions between SEs and promoters of different cell types in brain tissue previously The motifs of the screened seven genes matched to MEF2A (p < 0.001) were scanned by MEME tools (described in Methods section "Transcription Factors Enrichment") and were shown in UCSC. The positions and sequences of the screened seven genes were inquired in UCSC, and promoters (shown as green bar) were defined as the 2,000 bp window centered on the transcript start site of genes. The motif of MEF2A downloaded from JASPAR 2020 is shown on the bottom right. (C) Correlations between transcription factors and the expression levels of BNIP3, CDK5R1, HERC1, ITPR1, OPTN, UBQLN2, and USP33 were analyzed (Pearson's correlation) in the dataset GSE84422 (Wang et al., 2016) and was shown as the heatmap. (D) The dataset GSE84422 (Wang et al., 2016) was downloaded, and the mRNA expression (signal intensity) of MEF2A in different neuropathological category of sample (possible AD group, control group, and the definite AD group) were shown using GraphPad. ANOVA followed by LSD test was used for comparison among groups. *p < 0.05 compared with the control group. (E) The dataset GSE118553 (Patel et al., 2019) was downloaded and the mRNA expression (signal intensity) of MEF2A in different groups (AD, AsymAD, and control group) and brain regions of sample (cerebellum, entorhinal cortex, the frontal cortex and temporal cortex) were shown using GraphPad. An ANOVA followed by an LSD test was used for comparison among groups. *p < 0.05 compared with the control group. published in Nott et al. (2019). We found that ATAC-seq, H3K27ac, and H3K4me3 ChIP-seq showed significant peaks in the MEF2A promoter, with no cell type specificity. The H3K27ac ChIP-seq (a characteristic marker of enhancers and promoters) in microglia showed significant peaks in the MEF2A enhancers, but no significant peaks were observed in neurons or astrocytes. The PLAC-seq revealed a series of regions, mainly concentrated in microglia, which had a high degree of interaction with the promoter of MEF2A, suggesting only microglia had SEs ( Figure 6A). Therefore, this suggested that the specific expression of MEF2A in microglia was due to the SEs, which promoted the expression of MEF2A.
We analyzed the DMPs in the dataset GSE80970 (Smith et al., 2018), and the results were shown in a Manhattan plot (Figure 6B). In the SEs region between MEF2A and LRRC28, there were eight DMPs which had methylation differences (p < 0.05) in AD compared with healthy controls. The methylation levels of SEs in participants with AD were significantly higher than that in healthy controls (Figure 6C, p < 0.05). This suggested that AD may lead to an increase in the methylation of the SE region of MEF2A, thereby reducing the mutual binding with the MEF2A promoter, leading to a decrease in the expression of MEF2A.

DISCUSSION
In this study, we used bioinformatics methods to investigate the relationship between autophagy-related genes and pathological progression of AD. We downloaded the datasets and identified the seven autophagy-related genes to be correlated with pathological progression of AD using WGCNA. A correlation analysis showed that MEF2A, a transcription factor enriched in UCSC database, was closely related to the screened seven genes. Furthermore, MEF2A was highly expressed in microglia due to the existence of SEs, and the decreased expression of MEF2A in AD was caused by the increased methylation of the SEs.
As a chronic neurodegenerative disease, AD has been found to be closely related to autophagy, given that autophagyrelated pathology including accumulated autophagic vacuoles (AVs) were found in a number of dystrophic neurites in AD (Nixon and Yang, 2011). Through live-imaging studies of cortical neurons, this study showed that the inhibition of lysosomal proteolysis could selectively disrupt the axonal transport of autophagy-related compartments, causing an AD-like axonal dystrophy. Our analysis of the array dataset showed that many genes related to autophagy were significantly decreased in AD and closely related to the pathological progression of AD. Among them, gene expression in MEblack module was most strongly associated with CDR, CERAD, PLQ_Mn, NPrSum, and NTrSum, which can reflect the pathology of AD. Thus, genes in MEblack module, including BNIP3, OPTN, CDK5R1, UBQLN2, ITPR1, USP33, and HERC1, were selected for further analysis.
Currently, many studies showed that age-related declines in cognitive fitness were associated with a reduction in autophagy (Hou et al., 2019). However, according to our results, age was positively related to autophagy gene expression in MEgrey module. This may be related to the fact that participants were all over the age of 60. This was consistent with the The mRNA expression of MEF2A in different cell types was performed using GraphPad. An ANOVA followed by an LSD test were used for comparison among groups. *p < 0.05 compared with the microglia group. (D) The mRNA expression of MEF2A in microglia and in mixed brain cell from different brain regions in dataset GSE62420 (Grabert et al., 2016) were performed using GraphPad. An ANOVA followed by an LSD test was used for comparison among groups. *p < 0.05 compared with the microglia group. (E) The mRNA expression levels of MEF2A in AD and in healthy controls are shown. t-Test was performed for comparison between groups. *p < 0.05 compared with the control group.
report of Glatigny et al. (2019) which stated that autophagy was increased in many long-lived model organisms and contributed significantly to their longevity. The relationship between autophagy and longevity warrants further study.
For the genes in MEblack module, the functions of BNIP3, OPTN, and UBQLN2 are related to LC3 on the lysosomal membrane and participate in the fusion process between lysosomes and autophagosomes. BNIP3, a BH3-domain protein FIGURE 6 | Epigenetic regulation analysis of MEF2A expression. (A) UCSC browser of the MEF2A locus showed ATAC-seq, H3K27ac and H3K4me3 ChIP-seq, and PLAC-seq in brain cell types, including neurons, microglia, astrocytes, and oligodendrocytes (OLs). Shared active promoter region is indicated in yellow; microglia-specific enhancer region is indicated in blue. (B) The DMPs range around MEF2A, approximately 99,676,703-100,882,647 bp of chromosome 15, was used to create a Manhattan plot by coMET package of the R software with a threshold p-value of 0.05 (shown as red dotted line). The -log10(p-value) of DMPs for MEF2A in dataset GSE80970 (Smith et al., 2018) were shown in the Manhattan plot. The SEs region is shown as a green box. (C) The methylation levels for eight DMPs of the enhancer region in AD and healthy people were shown using GraphPad, and a Student's t-test was used for comparisons across groups. *p < 0.05 compared with the control group.
of the Bcl-2 family located mainly on the outer membrane of mitochondria, is an important mitochondrial autophagy receptor that can specifically bind to LC3 on the lysosomal membrane to promote the fusion of autophagosomes containing mitochondria with lysosomes to induce autophagy (Chourasia and Macleod, 2015;Tang et al., 2019). OPTN is a ubiquitinbound autophagy receptor involved in pathogen autophagy and mitochondrial autophagy (Bussi et al., 2018). In aged APP-PSEN1-SREBF2 mice, chronic cholesterol accumulation results in an age-dependent impairment of OPTN translocation to mitochondria, inhibiting mitophagosome formation (Roca-Agujetas et al., 2021). APP/PS1 mice had enhanced Aβ clearance, improved cognition and mobility when treated with miR-331-3p and miR-9-5p, two microRNAs targeting autophagy receptors SQSTM1 and OPTN, respectively, and antagomirs at a late stage (Chen et al., 2021). The function of UBQLN2 in the ubiquitin protease system was to direct misfolded or redundant proteins to the proteasome for degradation. Several studies have shown that UBQLN2 was involved in the process of autophagy and directly combined with LC3 to promote the fusion of lysosomes and autophagosomes (Osaka et al., 2015;Hjerpe et al., 2016;Chen et al., 2018). UBQLN2(P497H) transgenic mice, causes amyotrophic lateral sclerosis (ALS) and frontotemporal type of dementia, had the feature of a dendritic spinopathy with protein aggregation in the dendritic spines and an associated decrease in dendritic spine density and synaptic dysfunction, related to impaired protein degradation (Gorrie et al., 2014).
The other four genes were also closely related to autophagy. CDK5, a serine/threonine kinase, is essential for neuronal migration and synaptic plasticity. CDK5 activation depends on the protein p35 encoded by the CDK5R1 gene, which forms a complex with CDK5 to perform its biological functions (Roufayel and Murshid, 2019). CDK5 genetically interacts with Acinus (Acn), a primarily nuclear protein, which promotes starvation-independent, basal autophagy. Downregulation of CDK5 influences pathologic processes of AD, including the formation of amyloid plaques and tau hyperphosphorylation (Nandi et al., 2017;Nandi and Krämer, 2018). CDK5R1 determines the risk for AD, with a 12.5-fold decrease in AD risk associated with both homozygosity for CDK5R1 (3 -UTR, rs735555) A allele and homozygosity for GSK-3β (-50, rs334558) C allele (Mateo et al., 2009). Moncini et al. (2017) demonstrated that two microRNAs, miR-103 and miR-107, regulate CDK5R1 expression and affect the levels of p35. As the autophagy receptor of cells, ITPR1, a member of the IP3 receptor family, encodes the endoplasmic reticulum (ER) receptor and mediates the release of ER calcium to induce autophagy (Messai et al., 2015;Ren et al., 2017;Xu et al., 2019). Recently, Seo et al. (2020) found ITPR1 displayed associations with the neuroimaging features of AD pathologies through a targeted sequencing analysis of the coding and UTR regions of 132 AD susceptibility genes including 557 participants. USP33 is a deubiquitination enzyme associated with the regulation of lysosomal activity and cell membrane surface receptors (Kommaddi et al., 2015). HERC1, a giant protein belonging to the HERC family, is involved in regulating the ubiquitination of intracellular proteins and can interact with mTOR to regulate the autophagy process (Mashimo et al., 2009;Ruiz et al., 2016;Bachiller et al., 2018).
MEF2A, a member of the MEF2 family, belongs to the MADS-box superfamily and is involved in the transcription of many important genes in the cell life cycle in the form of dimers, including the growth, differentiation and apoptosis of neurons (Zhu et al., 2018). Vargas et al. (2018) used the transcription regulatory network and master regulator analyses on transcriptomic data of human hippocampus from GEO to identify transcription factors that can potentially act as master regulators in AD, and then 34 master regulator candidates were identified including MEF2A. Our results showed that MEF2A, a transcription factor for the seven autophagy-related genes, was significantly decreased in AD. González-Velasco et al. (2020) also reported that 158 genes were regulated by transcription factors MEF2A among the transcriptional changes in the cerebral cortex and hippocampus caused by aging. Genome-wide association studies (GWAS) were performed to assess the significance of the overlap between genome-wide significant AD risk variants and sites of open chromatin from data sets representing diverse tissue types (Tansey et al., 2018). AD risk variants of MEF2A were significantly enriched both in macrophage and microglia. González et al. (2007) stated that variation in the MEF2A gene could be involved in the risk of developing late-onset AD. Besides, MEF2C, another member of the MEF2 transcription factor family, identified by GWAS as also having effects on AD risk, was inferred as a MEF2A target. Moreover, using rapid time-lapse two-photon calcium imaging of network activity and single-neuron growth within the unanesthetized developing brain, Chen et al. (2012) found that MEF2A was the major regulator of neuronal response to plasticity-inducing visual stimulation directing both structural and functional changes.
Microglia, accounting for 10-15% of the total number of brain cells, are an innate immune cell in the central nervous system that uses phagocytosis to engulf apoptotic cells and cellular debris (Plaza-Zabala et al., 2017). Microglia participates in the regulation of tissue repair, synaptic plasticity and synaptogenesis, resist invasion of foreign pathogens, and maintaining the stability of central nervous system tissues. In our study, we found that MEF2A, as a transcription factor to regulate BNIP3, OPTN, and UBQLN2, was highly expressed in microglia. Phagocytosis is very similar to autophagy in vacuole formation and lysosomal digestion. However, unlike autophagy, which is present in all cells, phagocytosis is a unique function of innate immune cells such as microglia. The dysregulation of microglia is closely related to the pathological process of AD, in particular in the context of its role in the phagocytosis of amyloid plaques (Condello et al., 2015;Sacks et al., 2018). The potential regulatory effect of autophagy on phagocytosis may occur in different phases of the phagocytosis cascade, including phagocytosis of substrates, maturation of phagocytes, and fusion with lysosomes, thereby affecting the degradation of phagocytes. For examples, in the LC3-related phagocytosis process, autophagy is partially transferred to phagocytosis to promote the effective intracellular degradation of phagocytic extracellular substances (Martinez et al., 2016). Heckmann et al. (2019) found that the clearance of Aβ in microglia cells correlated with LC3-related phagocytosis, which can promote the phagocytosis efficiency of cells and reuse the membrane receptors related to Aβ phagocytosis, such as CD36, TERM2, and TLR4, to improve cognitive levels.
In the past, microglial phenotypes were characterized by cell surface molecules and were classified as M0, M1 like (exhibiting proinflammatory signaling and neurotoxicity), or M2 like (participating in the resolution of inflammation). However, with the help of newly developed technologies, including single cell RNA sequencing, quantitative proteomics, and epigenetic studies, the characterization of microglial diversity in health and disease has therefore been redefined (Ransohoff, 2016). Recent genome−wide transcriptomic analyses of microglial cells under different disease conditions have uncovered a new subpopulation named disease−associated microglia (MGnD, Verheijen and Sleegers, 2018;García-Revilla et al., 2019). Krasemann et al. (2017) found that MEF2A was significantly decreased in microglia of EAE (multiple sclerosis models), SOD1 (ALS models), and APP/PS1 (AD models) mice. The aggregation of amyloid β (Aβ) changed the M0-homeostatic microglial phenotype to the neurodegenerative phenotype MGnD identified by two major gene clusters, after which the expression of MEF2A was significantly decreased. TREM2 induced APOE signaling and targeting the TREM2-APOE pathway restored the homeostatic signature of microglia in ALS and AD mouse models and prevented neuronal loss in an acute model of neurodegeneration. Our results may explain that the changes in microglia homeostasis in AD were related to autophagy regulated by MEF2A. A genome-wide analysis of gene expression in microglia from different brain regions across the adult lifespan of the mouse was performed, revealing that there were region-specific transcriptional profiles and agedependent regional variability in gene expression (Grabert et al., 2016). The presence of microglial heterogeneity may underly the different expression patterns of MEF2A in the different brain regions.
The relationship between AD and methylation has been well investigated. Wu et al. (2008) assessed AD-related gene methylation in peripheral blood leucocytes of diagnosed AD, which revealed decreased DNA methylation at the amyloid precursor protein (APP) promoter regions accompanied by upregulated APP transcripts. The methylation of other genes, including BACE1, PSEN1, SORL1, and NEP, involved in amyloidogenic pathway, was also shown to be related to AD (Poon et al., 2020). The epigenetic mechanisms, including the methylation of GSK3β, BDNF, ANK1, BIN1, and RELN, have been consistently reported to play critical roles in neurochemical processes including long-term potentiation (LTP) and synaptic plasticity (Jager et al., 2014;Poon et al., 2020). We found that the MEF2A in AD had significantly higher levels of methylation in the SE region than in healthy controls. SEs are composed of a series of enhancers and are gene control elements with tissue specificity. The significant peaks of H3K27ac ChIP-seq (a characteristic marker of enhancers and promoters) in MEF2A enhancers and promoters in microglia illustrated the reason for which MEF2A is highly expressed in microglia. Changing the methylation levels of SEs can regulate their interaction with gene promoters and thus regulate the expression of related genes (Flam et al., 2019). This suggests that SEs play a key regulatory role in the expression of the MEF2A.
Currently, there is no direct evidence that the screened seven genes are significantly reduced in microglia from the snRNAseq dataset of GSE138852 (Grubman et al., 2019). However, this does not mean that the screened seven genes are not altered in microglia in AD. We speculate that the expression of these genes may be too low to be accurately detected in microglia, especially OPTN, CDK5R1, and BNIP3, based on the report of Li et al. (2019). Similarly, the low expression levels of these genes do not mean that they are dispensable to microglia. Indeed, these genes are crucial to the function of microglia. Bussi et al. (2018) found that autophagy induced by exogenous fibrillar in microglia correlated with lysosomal damage and was characterized by the recruitment of the selective autophagy-associated proteins TANK-binding kinase 1 (TBK1) and OPTN to ubiquitinated lysosomes. Ma et al. (2013) found that enhanced CDK5 activity by increasing p35-to-p25 conversion promoted Aβ phagocytosis in microglia, whereas the inhibition of CDK5 reduced Aβ internalization. Gong et al. (2020) reported that pinocembrin protected microglial cells against intermittent hypoxia (IH)induced cytotoxicity by activating BNIP3-dependent mitophagy through the JNK-ERK signaling pathway. In addition, another reason may be linked to the fact that the MGnD was not captured accurately; it was found mainly concentrated around the amyloid plaques and was not evenly distributed throughout AD brain tissue (Krasemann et al., 2017). Besides, whether the reduced expression of these genes in microglia was influenced by epigenetics mechanisms warrants further investigation.

CONCLUSION
In summary, our results indicated that AD is associated with the increased methylation levels of MEF2A enhancer, reducing the expression of MEF2A and downregulating the expression of autophagy-related genes which were closely related to AD pathogenesis, thereby inhibiting autophagy. Although further conformational studies are warranted, our findings provide further insights into the role of MEF2A in the prevention and treatment of AD. The association between the reduction of MEF2A expression and autophagy-related genes in AD warrants further investigations.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
YJ contributed to the perception and finally approved the version to be published. HL participated the whole work, drafted the article, and did the data analysis. FW and XG collected and downloaded the data. All authors contributed to the article and approved the submitted version. Supplementary Figure 3 | The mRNA expression levels of the screened seven genes in microglia. In GSE138852, the mRNA expression levels (logCounts) of BNIP3, CDK5R1, HERC1, ITPR1, OPTN, UBQLN2, and USP33 in AD and in healthy controls were shown. A t-test was performed for comparison between groups. * p < 0.05 compared with the control group.