Regulators relate to AD
To evaluate whether m6A regulators are related to AD, we first performed gene differential expression analysis between AD and control cases in different brain regions (Fig. 2A, p values showed in Supplementary table 1). Separate regulators’ results vary in different brain regions. Only 2 writers differed in the PG region of the postcentral gyrus: VIRMA and RBM15 were downregulated in AD groups. ZCCHC4 and HNRNPNPA2B1 were significantly up regulated in the superior frontal gyrus. Several regulators of different types are differentially expressed in the AD patients' hippocampus. Two writers (RBM15 and METTL3) were downregulated in AD samples. As for readers, IGFBP1 and ELAVL1 were up-regulated, and YTHDC1 and YTHDF2 were down-regulated in AD samples' hippocampus. Eraser FTO was downregulated in AD samples. In the entorhinal cortex, 3 readers (IGFBP3, YTHDC1, and YTHDC2) were up regulated to varying degrees. Only YTHDC1 was differentially expressed across multiple brain regions, and the propensity of changes was varied. In addition to the hippocampus, the other three brain regions have an obvious tendency to certain regulator types. In the PG region, m6A writers tended to be downregulated in AD group. Writers and readers tended to be up regulated in AD patients’ SFG region. And m6A readers tended to be up regulated in the EC region.
In addition to DEGs exploration, we also assessed which regulators might influence disease progression. We investigated the relationship between cognition level, neuropathology stages, and regulators (Fig. 2B, p values showed in Supplementary table 2). Multiple regulators were considered to affect the cognitive level or pathological progression in multiple brain regions. Among erasers, the level of FTO transcript in the PG region is associated with both MMSE score and Braak stage, negatively correlated with disease severity. The level of FTO transcript in the HP region is also associated with neuropathology, and negatively correlated with the Braak stages. For readers, ELAVL1 was positively correlated with the Braak stage in both the SFG and HP regions. The levels of YTHDF2 in the HP and EC regions were correlated with pathological progression and tended to have opposite effects and were also positively correlated with cognitive levels in the PG region. Among writers, the levels of HAKAI in HP and EC regions were negatively correlated with Braak score and MMSE score, respectively. There are many genes with altered expression in AD patients compared to controls and related to pathology and cognition at the same time, which may be indicating critical regulatory roles in AD, such as down-regulation of FTO and YTHDF2 and up-regulation of ELAVL1 also happening with more severe neuropathological progression in AD patients' HP region. Many genes correlating with cognition level and neuropathology were not previously found in DEGs screening. Despite no expression change in erasers between groups in the PG region, FTO and ALKBH5 were found correlate with the Braak stages and the MMSE scores. Similar results were found in other brain regions, suggesting that there may be potential contributions of regulators to the pathological deterioration of AD that don’t have significant differences in expression level between control and disease groups.
Regulators' interactions with each other are very common in many circumstances, so we also measured the correlation between regulators in different brain regions (p < = 0.005, Fig. 3A-H, Supplementary table 3). The correlations of regulators existed in all brain regions and were mostly positive correlations. However, the correlation status varied between brain regions. There were extensive intra- and inter-group correlations between writers and readers in the PG region, while erasers are only related to readers. In the SFG region, writers and readers are still commonly correlated. There was an increase in the number of correlations among writers compared to the PG region, while the correlation between readers decreased. Erasers were not only related to readers. FTO was also associated with HAKAI, an m6A writer. The HP region was similar to the SFG region, and a strong correlation between writers and writers was detected. There is a weak correlation between readers in the HP region, but ALKBH5 (eraser) was discovered to be related to writers. The number of correlations between regulators in the EC region was the smallest, but the erasers had a stronger correlation with writers and readers, which might indicate a different pathogenic mechanism in this region. In conclusion, the correlation of regulators varies in different brain regions. Considering the different functions of each brain region, this may be related to different pathogenic mechanisms of m6A and AD.
Biological process influenced by m6A regulators
To investigate which biological processes the above regulators may affect, the GSVA method was used to measure the activation of KEGG and REACTOME biological pathways and to assess the changes of these pathways and terms in AD and their correlation with regulators (Fig. 3, Supplementary Fig. 1A-D, Supplementary table 4). Several pathways with altered activation levels in AD were found in different brain regions (p <= 0.1). 164 in PG, 420 in SFG, 686 in HP, and 713 in EC. Among them, 23 pathways' scores differed in all 4 brain regions, covering several metabolism pathways such as N-glycan biosynthesis, lipid metabolism, amino acid metabolism, and protein export pathways. Several signaling pathways were involved, such as G protein activation and beta-gamma signaling, chemokine receptor binding, NTRK2 and NTRK3 activation via RAS, etc. We also evaluated the pathways potentially regulated by m6A regulators. For the regulators related to AD in the corresponding brain regions (DEGs and genes associated with Braak stage and MMSE score), the Pearson correlations between these regulators and the pathway GSVA scores were calculated. Among the 23 pathways, 13 were also associated with at least one disease-related regulator (p <= 0.05, Fig. 4, Supplementary Fig. 1H-I, Supplementary table 5) in all brain regions, mainly REACTOME pathways (11 out of 13). Those pathways included most of the metabolic pathways described previously (N glycan, amino acid, and protein metabolism) and are associated with most disease-related regulators in the hippocampus.
Regulators closely related to most of these pathways (associated with more than half of the 13 pathways) in 4 regions had a significant difference. YTHDF2, YTHDF3, and FTO in the PG region are associated with more than half of these pathways. HNRPA2B1 in the SFG region assumed this duty. YTHDF2 and FTO in the hippocampus still played an important role, with RBM15 associated with multiple pathways. In the EC region, only HNRNPC is associated with most pathways. Most of these regulators, other than FTO and RBM15, are readers. YTHDF2 and FTO played important roles in 2 of all 4 brain regions, indicating an important relationship with AD.
Different m6A modification patterns in AD patients
To investigate whether there are distinct m6A modification patterns in AD patients' transcriptome profiles in different regions, we used the regulators' expression profile to perform consensus clustering. The clustering results among 4 brain regions were divided to 2 patterns. The PG and SFG regions were divided into 3 clusters, while the HP and EC regions were divided into 2 clusters (Fig. 5, Supplementary table 6). To explore the consistency of the above groupings, we examined the similarities and the differences of differentially expressed regulators in different clusters of each brain region (Supplementary table 7). For the 3 clusters in each the PG region and the SFG region, 8 regulators (HNRNPA2B1, HNRNPC, IGFBP2, IGFBP3, YTHDC1, YTHDC2, YTHDF1, and ZC3H13) had similar expression patterns (Fig. 5B and D). For the 2 clusters in each the HP and EC regions, HNRNPA2B1 and IGFBP3 had similar expression patterns (Fig. 5F and H). The above regulators are all readers except ZC3H13.
To explore the biological differences of the different modification patterns above, we conducted a difference analysis on the GSVA scores of biological pathways between different clusters in each brain region. We found several KEGG and REACTOME pathways have different activation state in all 4 brain regions, some with the same trend among multiple regions (p < 0.1). 98 KEGG pathways and 797 REACTOME pathways had different activation state among clusters in the PG region. In the SFG region, 136 KEGG and 1095 REACTOME pathways were differentially activated. Of these two regions, 23 KEGG pathways and 189 REACTOMEs had similar trends among the three clusters, such as the Alzheimer pathway, the autophagy pathway, the oxidative phosphorylation pathway, and the TCA cycle pathway (Fig. 6A-B). There were differences in 81 KEGG pathways and 703 REACTOME pathways between 2 different modification patterns in the HP region, also including the Alzheimer pathway, autophagy pathway, and TCA cycle pathway (Fig. 6C). There are few differentially activated pathways in the EC region. Only 8 KEGG pathways and 72 REACTOME pathways have differences, presumably because of the lack of numbers of AD samples in this region. In summary, we found different m6A modification patterns can be distinguished in AD patients’ transcriptome in multiple brain regions, commonly affecting several pathways, including autophagy, oxidative phosphorylation, the TCA cycle, etc.
Discovery of key regulators based on WGCNA
To study the relationship between regulators and AD and assess their importance in the regulatory network, the WGCNA method was used to construct a weighted co-expression network for each brain region. The importance of regulators was evaluated based on the relationship between regulators and phenotype modules (Fig. 7). Regulators with GS > 0.2 for a phenotype and MM > 0.65 for a module were screened, and their GS values with all phenotypes and MM values with all modules are shown in Supplementary Fig. 2. Some important regulators related to AD/Control phenotype were found. FTO was considered essential in the PG region (Supplementary Fig. 2A-B). RBM15B and FTO had critical roles in the SFG region (Supplementary Fig. 2C-D). Multiple regulators greatly influenced the HP region, including RBM15, FTO, YTHDF2, YTHDC1, and IGFBP1 (Supplementary Fig. 2E-F). In the EC region, YTHDC1 and YTHDC2 played important roles in the regulatory network (Supplementary Fig. 2G-H). To assess the disease-related biological processes that the above regulators may modulate, disease-related (p < 0.1) and modules related to regulators (MM > 0.8 for regulators of corresponding modules) were selected in each brain region. The final selection was as follows: the PG region did not obtain disease-related modules corresponding to essential regulators. In the SFG region, RBM15B corresponds to the grey60 module, and FTO corresponds to the blue module. In the HP region, RBM15, YTHDC1, and IGFBP1 all correspond to the turquoise module, while FTO and YTHDF2 correspond to the blue module. In the EC region, YTHDC1 and YTHDC2 correspond to light cyan and dark red modules, respectively. To confirm the key genes in these modules regulated by m6A regulators, the hub genes in each corresponding module were selected for the following analysis (MM > 0.8 and GS > 0.2). We assessed whether these genes were correlated with the corresponding regulators above (p < 0.05) and obtained potential relationships between regulators and hub genes. We also evaluated the reliability of the possible regulatory relation based on the Me-RIP seq result from an article and potential ties between regulator and target from the m6A2targets database (Supplementary Fig. 3-5). In conclusion, we obtained key regulators and regulatory relationships that may have a decisive impact on AD.
Diagnose Model construction using key genes
For the key regulators and possibly regulated genes obtained in the previous step, use the above genes to fit a diagnostic model to evaluate the relation between these genes and AD in each brain region. Specifically, using the above regulators and hub genes (GS > 0.2, MM > 0.8) associated with regulators (p < 0.005), which were disease-related and highly correlated with disease-related modules. We performed Lasso penalized logistic regression in each brain region. Then the diagnostic model is validated using the validation dataset.
We obtained different diagnostic models in the regions above. No corresponding diagnostic model was constructed in the PG region due to lacking key regulators and target modules found. For the SFG region, we used RBM15B, FTO, and a total of 55 related genes to fit the model, with λ=-3.9. The model consisted of 13 genes, and the AUC was 0.958 (Fig. 8A-B). The frontal lobe data of the validation dataset GSE118553 had an AUC of 0.98 (Fig. 8C-D). The HP region training data set used RBM15, FTO, YTHDF2, YTHDC1, IGFBP1, and related 4524 genes for penalized regression, with λ=-6.1, the final model consists of 21 genes, and the AUC is 1. For the hippocampal data of the validation datasets GSE28146, the AUC was also 1. The EC region training dataset used YTHDC1, YTHDC2, and a total of 33 related genes for model fitting, with λ=-4.9, the model consists of 8 genes, and the AUC is 1 (Fig. 8E-F). The entorhinal cortex data of the validation dataset GSE118553 had an AUC of 0.731 (Fig. 3F). In conclusion, highly reliable diagnostic models for SFG, HP, and EC regions were obtained.