N6-methyladenosine (m6A) modification and its clinical relevance in cognitive dysfunctions

Collection of m6A regulators

Our study was designed and conducted according to the flow chart (Figure 1C). Briefly, twenty-six m6A regulators were selected accordingly to recent publications [15, 34–37], including ten writers (METTL3, METTL14, METTL16, RBM15, RBM15B, WTAP, KIAA1429, PCIF1, ZCCHC4, ZC3H13), two erasers (FTO, ALKBH5) and fourteen readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2, HNRNPC, FMR1, RBMX, LRPPRC, ELAVL1) (Figure 1B). Figure 1A summarized the landscape of related regulators, including classification, biological functions, and molecular mechanisms.

Acquisition of microarray datasets and preprocessing

The gene-expression dataset with full clinical annotation was obtained from Gene Expression Omnibus (GEO), a publicly sponsored genomic database operated by the National Institutes of Health (NIH). GEO provides open access to many gene expression data from biological and statically comparable samples [38]. In total, four eligible datasets regarding AD, VD, or MCI, including GSE122063, GSE63060, GSE63061, and GSE84422, were selected. There are 711 blood samples from GSE63060 (n = 382) and GSE63061 (n = 329), including 238 control, 189 MCI and 284 AD. GSE122063 consists of 136 brain samples from either frontal or temporal lobe, and patients are divided into control (n = 11), VD (n = 8), and AD (n = 12) groups. GSE84422 collected 1053 post-mortem brain samples from 125 subjects with a full spectrum of AD.

The microarray platforms provided by Illumina and Agilent were downloaded in the format of normalized matrix files. The dataset retrieved from Affymetrix was downloaded in raw “CEL” form. The R/Bioconductor algorithm “RMA” and package “SVA” were used to preprocess gene chips normalization among datasets and to remove batch effects and other latent variations [39]. The overall workflow was presented in Figure 1C.

Analysis of unsupervised clustering for m6A regulators

Unsupervised clustering analysis is an effective machine learning tool for exploring the patterns of datasets in a complex system, which has been applied to AD studies [40] and single-cell RNA sequencing applications [41]. In the current study, the m6A-related regulators were classified into several distinct endotypes by employing unsupervised clustering methods [42, 43] and m6A modification patterns based on the mRNA sequencing of 21 m6A regulators were hereafter determined for further research. Consensus Cluster Plus R package was conducted to perform 1000 times repetitions to guarantee the stability of classification [44]. The number of clusters was determined by the consensus clustering cumulative distribution function (CDF) result (Supplementary Figure 1). The purpose of the CDF plot is to find the k at which the distribution reaches an approximate maximum, which indicates maximum stability, and after which divisions are equivalent to random picks rather than the true cluster structure [44]. Besides, patients were classified into different groups for deeper analysis by adopting an unsupervised clustering method for analyzing the significant difference in different clusters by consensus clustering.

Identification of differential gene expression

After data normalization, differentially expressed genes (DEGs) in datasets of GSE122063, GSE63060, and GSE63061 were identified using the “Limma” package from R/Bioconductor software [45]. The significance of DEGs was set as the adjusted P-value <0.05 and threshold of |log2FC|≥1. Different expression levels of 26 m6A regulators among groups were further verified by unsupervised clustering analysis.

Exploring KEGG pathway enrichment

After identifying DEGs in frontal and temporal cortices, Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/, ver. 6.8) was further used to identify Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of DEGs in the above two brain regions. The cutoff criteria were set as P values of 0.05. The results of top10 KEGG pathways in both cortices were picked up and constructed in a bubble plot via R Studio.

Gene set enrichment analysis and functional annotation

Gene set enrichment analysis (GSEA) has successfully been applied to interpret the molecular pathway activated in different biological states [46]. In this study, software “GSEA” (https://www.gsea-msigdb.org/) was utilized to identify the gene up- or down-regulation after filtering for gene set size (min = 5, max = 500) and ranked by t-score [47]. The gene sets of “c2.cp.kegg.v7.1.symbols” (MSigDB database) were used for GSEA analysis. The FDR-corrected q-value <0.25 and P-value <0.05 were set for significance.

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) was used to extract highly correlated clinical traits and calculate module membership measures from the data sets [48]. “WGCNA” R package was applied to determine hub genes and clinical traits-related modules among microarrays [48]. Genes with variations in the top 25% were extracted from DEGs analysis. Biweight miscorrelation (corType = “Pearson”) was set to detect the outliers. The topological overlap matrix (TOM) was transformed to find the connectivity in the adjacent matrix. Genes were after that divided into multiple sensitive modules according to the TOM-based dissimilarity measurement. Other analysis setting regarding the identification of key modules included soft-threshold power = 7, scale free R2 = 0.9, height = 33, cut height = 0.2, and minimal module size = 10. Subsequently, genes from the highest correlated module were picked up to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. Hub genes were determined by defying gene significance (GS) >0.3 and module membership (MM) >0.8.

Comparing the expression levels according to genotypes of apolipoprotein E gene

The association between apolipoprotein E (APOE) gene ɛ4 allele and m6A methylation regulators in AD patients was examined by extracting data from GSE29652. Eighteen post-mortem brain samples of AD were categorized into APOE ɛ4+ or APOE ɛ4− subtype. The expression levels of m6A regulators were compared between APOE ɛ4 genotype groups after data normalization and DEGs extraction.

Statistical analysis

The Statistical Package for the Social Sciences (SPSS) version 24.0 was used for statistical analysis. Patients in GSE63060 and GSE63061 were sub-grouped by age according to their cognitive functions. Normal distributed continuous variables were described using mean ± SD; categorical variables were presented as percentages (%). Differences between groups were compared by t-test, one-way ANOVA, or Kruskal-Wallis test for continuous variables, and chi-square for categorical variables.

Correlation analyses were carried out to compute the strength of interrelationships between clinical traits and gene expression traits. Correlations between m6A regulators were computed by Spearman correlation analyses and visualized by the “corrplot” package in the R program. Univariate analysis examinations, filtering the meaningful independent variables, followed by multivariate logistic regression analysis, were conducted to estimate the association between m6A methylation levels and MCI and AD.

All statistical P values were two-tail, and p < 0.05 was regarded as statistically significant.

Data sharing statement

All relevant data supporting the key findings are available from the corresponding author upon reasonable request.

Overview of included datasets and m6A regulators

Seven hundred and seventy-one blood samples from GSE63060 and GSE63061 were stratified into three age categories (≤70, 71–79, ≥80). There are significant differences in age distribution among CTL, MCI, and AD groups (χ² = 26.2, P < 0.001) (Supplementary Table 1). 50% of patients over 80 years had AD, which is significantly higher than in younger age groups, supporting that AD is age-related.

Expression patterns of m6A-related regulators vary with cognitive dysfunctions

Firstly, expression levels of m6A RNA methylation regulators from different samples in the same groups were compared. As we can see, nearly half of the proteins have various expressive phenotypes in the control group. The expression profile of YTHDC2 mainly was enriched in the frontal lobe, while RBMX and FTO were in the temporal lobe, and IGF2BP3 was in the blood (Supplementary Figure 2A). The expression profiles were different in AD patients, and even more prominent expressive differences were found. The expression levels of METTL3 and KIAA1429 were highest in the frontal lobe, while PCIF1 had the highest expression in the temporal lobe, and YTHDC2 had the highest expression in the blood (Supplementary Figure 2B). The combined results indicated that the modification phenotypes of m6A-related regulators varied in different samples.

To further explore these m6A regulators’ expression patterns, our attention was fixed on comparing different levels of cognitive functions within the same tissues. Unsupervised clustering analysis was performed by classifying into several clusters, according to K-means. When classified according to brain regions, two distinct clusters were identified. We noted that m6A Cluster1 was relatively enriched in the frontal lobe (Supplementary Figure 2C), with writers METTL3, METTL14, readers YTHDF2/3, and YTHDC1/2 had the most significant enrichment (Supplementary Figure 2D–2K). No noticeable enrichment difference was found by visual inspection in m6A Cluster2. When classified by cognitive status, the exhibition was different (Figure 2). We discovered that m6A Cluster1 was distinctively enriched in the CTL group, especially in readers and eraser FTO (Figure 2A, 2C). As for m6A Cluster2, it was mainly enriched in readers of VD and AD groups, with IGF2BP1/2/3 and HNRNPA2B1 were the most obvious (Figure 2A, 2C). The most significant divergently expressed regulators between the three groups were seen in Figure 2E–2H (YTHDF2, F = 10.612; YTHDC2, F = 36.231; LRPPRC, F = 13.354; FTO, F = 22.789). Besides, our previous research has found the differences in epigenetics modifications in diverse brain regions (unpublished), so when further comparisons were conducted based on brain domains, the expression imbalance was similar to the combined results (Supplementary Figure 2C).

Figure 2.

Mutation frequency distribution of m6A regulators across different groups in brain and blood samples. (A) Unsupervised clustering of 26 m6A regulators in GSE122063, annotated according to cognitive statues. Red represented high expression of regulators and blue represented low expression. (B) Unsupervised clustering of 26 m6A regulators in GSE63060 and GSE636061, annotated according to cognitive statues. Yellow represented high expression of regulators and blue represented low expression. (C–D) The expression profiles of 26 m6A regulators in brain and blood samples. (E–L) Box plots showing representative expression differences between CTL, VD and AD from brain sample, or between CTL, MCI and AD from blood sample. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001, ^****p < 0.0001, ns: no significance.

There were three m6A Clusters, and 13 differently expressed regulators between CTL, MCI, and AD in blood samples (Figure 2B, 2D). Cluster 2 was mainly enriched in the VD group, Cluster 3 was mainly enriched in the CTL group, with WTAP (F = 22.354), RBM15 (F = 12.609), PCIF1 (F = 13.629), HNRNPC (F = 12.133) had the most distinctive mutations among groups (all Ps < 0.0001, Figure 2I–2L). Interestingly, we noticed that the expression tendencies of MCI were always different from those of CTL or AD groups (Figure 2D), indicating that expressional alterations of m6A regulators played a crucial role in mediating the progression of cognitive dysfunctions.

Interactions between m6A-related regulators

Considerable evidence has proved that cross-talks of m6A regulators were ubiquitous [35, 49]. Therefore, we further investigated collaboration among writers, erasers, and readers by constructing protein-protein interaction (PPI) networks. In general, correlated expression patterns and genetic alterations were not only within the same biological regulators but also among writers-erasers-readers (Supplementary Figure 3). Except for recently reported proteins PCIF1, ZCCHC4, and LRPPRC, the other regulators' expressive patterns were significantly associated with each other. Mostly evident interactions were observed between proteins within the ten writers. Close interrelationships were also found between writers with readers or erasers. Relatively few interactions were existed between the 14 readers, while complicated interactions between erasers FTO and ALKBH5 were apparent in the PPI network (Supplementary Figure 2).

Correlations between m6A regulators in different samples and two brain cortices were computed by Spearman correlation analyses and visualized by the “corrplot” package in the R program. The lines linking regulators showed their interactions, and thickness showed the correlation/interference strength between regulators (Figure 3, Supplementary Figure 3). As we can see, the result of brain samples showed strong co-occurrences between eraser FTO with readers RBMX, YTHDF2, and YTHDC1. At the same time, negative correlations were exhibited with writers METTL16, IGF2BP3 and HNRNPA2, and PCIF1. The closest interference was found between METTL16 and IGF2BP3 (Figure 3A, 3C; R = 0.792, p-value < 0.0001). According to brain regions, we found more inseparable expression profiles (Supplementary Figure 4). For example, the absolute values of the correlation coefficient between IGF2BP3 with FTO and YTHDF2 were both more than 0.80 in the temporal lobe (Supplementary Figure 4B, 4D; R = −0.818 and −0.813, respectively, both Ps < 0.0001). While in the frontal lobe, the interaction strength between regulators was somewhat weaker, but the absolute values of the correlation coefficient between IGF2BP3 and YTHDF2, FTO and YTHDC1 were also close to 0.80 (Supplementary Figure 4A, 4C; R = 0.773 and −0.753, respectively, both Ps < 0.0001). Besides, the interactions between METTL16, TYHDF2, YTHDC1, and IGF2BP3 with other regulators were always the strongest, indicating intimidate connections among m6A-related regulators.

Figure 3.

Interaction among m6A RNA methylation regulators. (A, B) Spearman correlation analysis of the 26 m6A methylation regulators in brain (A) and blood (B) samples. Positive correlation was marked with red and negative correlation with blue. (C, D) The interaction between the 26 m6A methylation regulators by constructing PPI network in brain (C) and blood (D) samples, respectively. The circle size was determined by the degree value. The lines linking regulators showed their interactions, and thickness showed the correlation strength between regulators. Positive correlation was marked with red and negative correlation with blue.

As for blood samples, we found that the overall correlations were less close to those in the brain. The writer RBM15 was positively correlated with writer WTAP and readers YTHDC1/2 and YTHDF3, while negatively correlated with writer PCIF1 and readers IGF2BP1/2, with the closest relationship existed within the same functional protein class between writers RBM15 and PCIF1 (Figure 3B, 3D; R = −0.638, p-value < 0.0001).

Identification of m6A-related biological functions

Potential biological processes of m6A-related regulators were further investigated. 502 and 674 DEGs were identified in frontal and temporal cortices, respectively, and they were subsequently performed by the Kyoto Encyclopedia of Genes and Genomes (KEEG) pathway analysis. Bubble plots showed that DEGs in two lobes were both mainly enriched in neuro-modulatory activities, including neuroactive ligand-receptor interaction, serotonergic synapse, GABAergic synapse, glutamatergic synapse, and cholinergic synapse. Other pathways, such as the retrograde endocannabinoid signaling pathway and Calcium signaling pathway, were also mainly involved (Figure 4A, 4B). In addition, the Ras signaling pathway, MAPK signaling pathway were typically enriched in the frontal lobe (Figure 4A), while Rap1 signaling pathway, phagosome, and Gap junction were in the temporal lobe (Figure 4B), demonstrating that biological functions in various brain regions were similar in general but also slightly different.

Figure 4.

Identification of differentially regulated molecular pathways and m6A-related biological functions. (A, B) Functional annotation of the genes with different expression in frontal (A) and temporal (B) lobes using KEGG pathway. (C–J) Representative differentially regulated pathways are shown by analyzing blood samples. Pathways with increased activation included pathways indicative of Alzheimer's Disease (C) and VEGF signaling pathway (E), pathways with reduced activation included Huntington's Disease (D) and taste transduction (F) in AD group compared with CTL group; (G–J) Pathways with dysregulated activation in MCI group compared with CTL group.

Estimation of differentially regulated molecular pathways

We used the Gene Set Enrichment Analysis (GSEA) algorithm to identify molecular pathways differentially regulated according to cognitive dysfunctions. Of the 65 pathways in the AD group compared to the CTL group, 41 pathways were activated, while 24 pathways were down-regulated (Additional file 1). Generally, GSEA-based analysis highlighted a broad dysregulation of genes related to neurodegenerative disorders. Among the latter are genesets associated with AD, Parkinson’s Disease and Huntington’s Disease, taste transduction, sphingolipid metabolism, and calcium signaling pathway (Figure 4C–4J). It is interesting to find that taste transduction was the only pathway mediated between AD and MCI, and it was up-regulated in AD vs. CTL. However, down-regulated both in MCI vs. CTL and AD vs. MCI (Figure 4J, 4F, Additional file 1), indicating that taste transduction malfunction might play a role in cognitive impairment progress. Besides, Pathways such as the VEGF signaling pathway, complement and coagulation cascades, JAK-STAT signaling pathway, and MAPK signaling pathways were activated in the AD group compared to the CTL group (Additional file 1). The identified dysregulated pathways and related genes, e.g., VEGF, JAK, MAPK, and complements, could be therapeutic targets of AD.

Construction of gene co-expression network and module of interest identification

In this study, weighted gene co-expression network analysis (WGCNA) [50] was performed to identify the key modules most associated with AD and MCI clinical traits. Ages were stratified into three subgroups: ≤70, 71–79, and ≥80 years old. After setting the soft threshold = 7, 6 outliers were removed (height = 33), and the left 705 genes were, after that, theoretically classified according to the expression pattern (Supplementary Figure 5). Firstly, we identified 11 modules of highly co-expressed genes by considering clinical traits, including age/age stratification, gender, cognitive status (Figure 5A). Eigengene adjacency heatmap was further conducted and revealed four main branches among the genes, which verified the above interconnections (Figure 5B). Unique color identifiers were assigned to each module, with gray represented the remaining poorly connected genes. Then, the co-expressed gene network was constructed, and the topological overlap matrix (TOM) heatmap plot was employed to show the network landscape (Figure 5C). The rows and columns in the TOM plot corresponded to various genes. The color intensity represented values of Pearson correlation coefficients, which meant that the higher color intensity indicates higher co-expression similarity between genes included in the network. Herein, genes appeared highly interconnected between module green with modules turquoise and blue, module purple with module brown by visual inspection (Figure 5C).

Figure 5.

Identification of key modules correlated with clinical traits by WGCNA. (A) Cluster dendrograms of all genes, with dissimilarity based on topological overlap, and then various module colors were assigned. (B) The upper panel displays the hierarchical clustering dendrogram of hub genes that summarize the analyzed modules and branches of the dendrogram group with eigengenes are closely correlated. The lower panel shows the eigengene adjacency heatmap, with the trait weight included. The darker red color represents higher adjacency, while darker blue color represents low adjacency. (C) Heatmap plot of Topological Overlap Matrix (TOM) among selected genes. Each module corresponds to a branch in the hierarchical clustering dendrogram. Modules demonstrate more saturated yellow or even red colors indicate higher co-expression interconnection. Genes locate at the tip of each branch indicate highest interconnection with the rest of the genes in the module. (D) Heatmap of the associations between module eigengenes and clinical traits. Each row and column correspond to a module eigengene or a clinical trait. The plot is colored by corresponding correlation according to the legend, and each cell contains the corresponding P-value. The red color represents positive correlation, while blue color represents negative correlation.

To further strengthen the study of crucial module identification, we defined a measurement to discern the statistical significance between modules with clinical traits. As depicted in the Heatmap, the associations between module eigengenes and clinical traits were colored by corresponding plot, with a darker color (red or blue) indicating strong correlations. Therefore, modules turquoise and yellow were closely intimated with cognitive statues (Figure 5D), indicating that genes in these two modules had high relationships with different cognitive performance levels. The investigation of gene significance (GS) was further performed to find the modules most biologically connected to clinical traits. (Figure 6A–6C). What is more, the existence of significant correlations between GS and module membership (MM) implied that genes within the module turquoise tended to strongly interrelate to cognition (correlation coefficient = 0.35, P = 3.4e-49; Figure 6B).

Figure 6.

Identification of hub genes and functional annotation of the WGCNA module highly correlated with clinical traits. (A–C) Scatter plots of eigengenes in the representative modules yellow (A), turquoise (B), and brown (C), which have highly significant correlation between Gene Significance (GS) and Module Membership (MM). (D) PPI network of the top20 hub genes in the module turquoise. The circles with darker red represent higher gene rank. (E) Biological functional annotation of the top20 hub genes in the module turquoise by Geno Ontology (GO) enrichment analysis.

A total of 1666 genes in module turquoise and the top20 hub genes, including SNRPG, SNRPD2, RPL26, ERH, SNRPB2, and SSB, were selected, setting GS >0.3 and module membership (MM) >0.8. We can tell that these genes were inseparably associated with each other (Figure 6D). Herein, GO enrichment analysis was performed to identify potential biological functions of module turquoise-related genes (Figure 6E). The result revealed that genes within the module turquoise were most significantly enriched in translation, peptide biosynthetic process, nuclear-transcribed mRNA catabolic process, mRNA/RNA catabolic process.

Prediction of associations between m6A-related methylations and cognitive dysfunctions

In order to assess the effect of m6A-related regulators on the prevalence of cognitive dysfunctions, unconditional logistic regressions were used. As expected, age predicts AD (OR = 1.055, 95%CI 1.025–1.087) and MCI (OR = 1.056, 95%CI 1.022–1.091), while gender has no prediction values. WTAP, ZCCH4, and HNRNPC from blood samples were found to have protective roles in MCI prevalence (Supplementary Table 2-1). While AD, METTL3, YTHDF3, and ALKBH5 were predicted to significantly increase the odds ratios, with YTHDF3 having the most apparent effect (OR = 8.033, 95%CI 2.047–31.523). On the contrary, METTL14, WTAP, YTHDC1, IGF2BP2, HNRNPC, and FTO had a protective effect on AD (Supplementary Table 2-2). When we further analyzed the relationships between MCI and AD, we surprisingly found that males with MCI had prone to deteriorating into AD compared to females (OR = 1.490, 95%CI 1.003–2.213), while age failed to predict the outcome. Besides, METTL3, WTAP, and RBM15 increased the possibility of occurrence, with METTL3 had the most apparent effect (OR = 6.984, 95%CI 1.778–27.432), while YTHDF1, YTHDC1, and LRPPRC were just the opposite (Supplementary Table 2-3).

Further analysis from brain regions of frontal and temporal lobes demonstrated age (OR = 1.067, 95%CI 1.050–1.085) and IGF2BP2 (OR = 2.633, 95%CI 1.248–5.552) were promotive to AD, while METTL16 and LRPPRC were protective to the prevalence. The model’s overall predictive value, sensitivity, and specificity were 79.4%, 70.8%, and 80.0%, respectively.

Exploration of correlations between APOE ɛ4 and m6A-related methylations

The ɛ4 allele of the human apolipoprotein E gene (APOE) is a well-proven genetic risk factor for the late-onset form of AD [51]. Studies have proved that the ɛ4 allele of APOE was differentially methylated in AD [52, 53]. However, whether RNA methylation changes the presence of the APOE ɛ4 allele is still unknown. Thus, data from GSE29652, which contains 18 postmortem brain samples, were extracted and analyzed according to APOE ɛ4+/−. As a result, five DEGs were dug out (Figure 7). Except for YTHDC2 was downregulated in the ɛ4+ group, the other four DEGs (METTL3, METTL16, RBMX, and LRPPRC) were upregulated epigenetic alterations of RNA methylation might be related to APOE ɛ4 dysfunction in AD (Figure 7A–7E).

Figure 7.

Expression differences on m6A-related regulators between APOE ɛ4+ and APOE ɛ4− groups. (A–E) Five representative DEGs between two groups. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.