Re-exploring the core genes and modules in the human frontal cortex during chronological aging: insights from network-based analysis of transcriptomic studies

Overlap of differentially expressed genes among studies within the compendium

To identify a common transcriptional signature reflecting age within the frontal cortex, we built a gene expression compendium using five independent studies (Table 1). Data were extracted and annotated, yielding a compilation of 7274 unique genes from 316 individuals, including 161 younger and 155 older individuals. An integrative analysis across datasets was then performed by computing the differentially expressed genes per dataset and assessing the overlap of the significant results. Eight genes (HTR7P1, MMD, OLFM1, ATP2B2, ANXA4, LPL, VCAN, and RHOBTB3) were significantly associated with age among the five datasets (Figure 1 and Table S1). Among them, LPL (lipoprotein lipase) is mainly expressed in adipose tissue and regulates the availability of polyunsaturated fatty acid within the central nervous system. HTR7P1, OLFM1 and ATP2B2 are highly expressed specifically in brain.

Figure 1. Significant age-associated genes in studies of the frontal cortex compendium. Shown is a Venn analysis performed to determine the overlap of significantly age-associated genes identified in five independent studies.

Meta-analysis of the frontal cortex throughout the normal aging

Five microarrays were analyzed using Integrative Meta-Analysis of Expression Data (INMEX), a web interface for integrative meta-analysis. The overall meta-analysis workflow in this study is shown in Figure 2A. By employing three meta-analysis methods, Fisher's method, Fixed effect model and Voting count, we identified 2367, 1856 and 1416 differentially expressed genes, respectively (Combined p < 0.001 or vote counts ≥ 2 were considered to be significant). Among those, 1260 genes were identified by all three methods (Figure 2B). A Venn diagram of these data was shown in Figure 2B. Among this group, 635 (50.4%) genes were downregulated and 625 (49.6%) were upregulated in older group as compared to the younger group. A heat map visualization of the top 20 genes across the different studies is displayed in Figure 2C. Among them, ANXA4 (annexin A4) was the top upregulated gene, and CALB1 (calbindin1) was the most significantly downregulated gene across the five microarray datasets. A complete list of the differentially expressed genes is provided in Table S2. The merged data of this meta-analysis is listed in Data S1.

Figure 2. Meta-analysis of the frontal cortex throughout normal aging. (A) Flowchart of the meta-analysis. (B) Venn diagram of differentially expressed genes identified from the meta-analysis using Fisher’s method, the vote counting method and a fixed effect model. (C) Heat map representation of the top 20 differentially expressed genes across different microarrays identified from the meta-analysis (row-wise comparison).

Functional enrichment of differentially expressed genes analysis

We next identified the enriched functional groupings that were significantly associated with age in the human frontal cortex. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using the set of up- and downregulated genes (False Discovery Rate (FDR) < 0.05) (Figure 3, Table S3). Among the enriched pathways were a number involving protein kinase signaling, including “PI3K-Akt signaling pathway” and “MAPK signaling pathway”. Enriched pathways were also related to chemical synapses (“Retrograde endocannabinoid signaling,” “Glutamatergic synapse” and “Dopaminergic synapse”) and electrical synapses (“Gap junction”). More specifically, MAPK and calcium signaling were shown to be involved in human longevity in an independent genome-wide association study (GWAS) of Han Chinese [19]. On the other hand, signal transduction systems that mediate “Long-term potentiation (LTP),” “Long-term depression (LTD)” and “Calcium signaling pathways” showed age-related downregulation. Using BiNGO in Cytoscape, we were able to obtain a global perspective of the changes in gene expression patterns (Figure 4). The gene sets showing upregulated expression were enriched for biological processes and molecular functions associated with “response to stimulus,” “cell migration,” “programmed cell death” and “apoptosis.” By contrast, the gene sets showing downregulated expression were enriched for biological processes and molecular functions associated with “synaptic transmission,” “phosphate metabolic process” and “learning and memory.” Given the functional enrichment of these differentially expressed genes, we delved further into the results in an integrative meta-analysis.

Figure 3. KEGG pathway analysis. KEGG pathway analysis was performed using the set of genes showing upregulated (red) or downregulated (green) expression (FDR< 0.05). The x-axis represents the fold enrichment (indicating the magnitude of enrichment in our dataset against the population background based on analysis using DAVID bioinformatics tools). The y-axis shows the pathway categories.

Figure 4. Gene Ontology analysis of significantly age-associated genes in the frontal cortex. Shown is a BiNGO (the Biological Network Gene Ontology tool) analysis depicting upregulated (A) and downregulated (B) genes in over-represented categories in the ontology of biological processes. The size of the nodes is proportional to the number of genes in the test set that are annotated to that node. Colored nodes are significantly over-represented along a color scale ranging from yellow (p = 0.05) to dark orange (p = 5.00E-7).

Analysis of PPI networks of selected genes

A novel integrative PPI network-based approach was applied to explore all the corresponding genes found in our meta-analysis search. This method has been validated as a way of identifying age-associated bio- data [13]. Considering the betweenness (BC) and degree of centrality (DC), we ranked the 1260 genes that showed significant differences between young and old in our integrative meta-analysis. We screened 33 hub genes that had at least 10 node degrees and ranked the top 10 percent of the total genes. Among them, 14 genes were upregulated and 19 genes were downregulated in the older group as compared to the younger group (Figure 5, Table S4). CDC42 was the most significantly downregulated hub gene (BC = 12193; DC = 30), followed by MAPK1 (BC = 9241.11; DC = 29) and HRAS (BC =7136.24; DC = 27). HDAC1 was the centermost upregulated gene (BC = 6305.53; DC = 25) followed by YES1 (BC = 6574.41; DC = 22).

Figure 5. Network analysis of significantly age-associated genes in the frontal cortex. (A) Zero-order interaction network of upregulated genes (red). (B) Zero-order interaction network of downregulated genes (green).

Neural aging is regulated through extensive interaction among genetic networks, signaling pathways and cellular metabolic responses, not through any single gene. However age-related genes are not spread throughout the interactome; instead, they cluster into tightly connected modules [20]. The joint expression of modules plays critical roles in human aging. After mapping the 1260 significant genes onto the PPI network, the zero-order network contained five subnetworks, including one big subnetwork (“continent”) (nodes 373, edges 743) and four smaller ones (“islands”) (p < 0.05, gene numbers ≥ 5). The big subnetwork was divided into 18 significant modules (p < 0.05, gene numbers ≥ 5) (Figure S1). In the following steps, we applied Gene Ontology (GO) analysis to these subnetworks and modules. The enriched terms characterizing the four satisfactory PPI islands and eighteen significant modules are shown in Table S5. The genes in module 1 were enriched for “apoptotic process.” Those in modules 2 and 3 were enriched for “regulation of transcription from RNA polymerase II promoter” and “regulation of small GTPase mediated signal transduction,” respectively. These enriched terms also included “semaphorin-plexin signaling pathway,” which is brain-specific expression and controls presynaptic neurotransmitter release and homeostatic plasticity [21]. In addition, three subnetworks were related to the mitochondria. Mitochondrial dysfunction is a critical characteristic of the aged brain and neurodegenerative diseases. This includes reduced respiration, dynamic structural modifications, loss of membrane potential and the changes of proteomic profile [22]. Structurally or functionally damaged mitochondria are more proficient at producing reactive oxygen species (ROS) but less efficient at ATP production, which accelerates pathological processes such as AD [23]. The remaining subnetwork was enriched for “mRNA splicing via spliceosome.” Of all the human organs, the brain has the largest amount of alternative splicing [24], and the genes that are alternatively spliced have been associated with neurodegenerative disorders such as AD [25].

Co-expressed PPI modules are enriched in GenAge human aging genes

To validate the candidate modules and nodes, we determined whether GenAge, a database providing a comprehensive overview of aging-related genes in humans and model systems, contained these aging-related genes. Human aging genes were enriched within module 1 (Odds ratio (OR) = 11.3, 95% Confidence Interval (CI): 4.73-27.16, p < 0.001) and module 16 (OR = 13.1, 95% CI: 2.35-73.53, p = 0.019) (Table S6). We also found that human aging genes are enriched within the big subnetwork (OR = 4.57, 95% CI 2.38-8.76, p < 0.001) (Table S5). Notably, the genes contained in GenAge had higher centralities among all 1260 genes, and almost half had at least 10 degrees of centrality. The higher centralities indicate that the function of these genes may be more variable and essential in aging.

Verification of functional roles using an established aged mouse model

To further verify the candidate genes screened with this network-based approach, we employed an established aged mouse model. Wild-type C57/BL6 mice maintained in a specific pathogen free (SPF) environment were used. Six-month-old mice comprised the young mature group, while 20-month-old mice comprised the old group. The Morris water maze was employed to evaluate the spatial learning and memory of two groups, which were assessed for 5 days. The escape latency significantly differed between the young and old mice (Day 3, p = 0.0313; Day 4, p = 0.0003; Day 5, p = 0.0015; Figure 6A), as did the numbers of crossings to the correct platform (Figure 6B) and the times in the correct quadrant on the fifth day (Figure 6C). qPCR performed to examine the mRNA expression in this model revealed that levels of Hdac1 and Yes1 mRNA were increased, while Cdc42 mRNA was decreased in the frontal cortex of the old group as compared to the young group, which was consistent with the results from the network analysis. By screening all age information in the GSE71620 dataset, we found that HDAC1, CDC42 and YES1 are associated with age. Moreover, by analyzing the GSE48350 dataset we found that HDAC1 and YES1 mRNAs were upregulated and CDC42 mRNA was downregulated in the frontal cortex of AD patients as compared to healthy subjects (Figure 6D-F). We then tested whether the three genes have an intrinsic relationship in aging or AD. By analyzing the STRING database, we found that there is only one significant connection between YES1 and CDC42 (Figure 6G). In both neural aging and AD (GSE71620 and GSE48350 datasets), CDC42 negatively correlates with YES1 (Figure 6F). We also verified 30 other genes using an established aging mouse model, and 18 were consistent with the results from the network analysis (Figure S2).

Figure 6. Changes in CDC42, HDAC1 and YES1 mRNA expression in the frontal cortex of an aging mouse model. (A) Escape latency on each day revealed a significant difference between younger and older mice on day 3 (p = 0.0313), day 4 (p = 0.0003), and day 5 (p = 0.0015) (n = 6). The average search errors of the younger and older mice did not differ in the first training trial (p = 0.3409) (n = 6). (B) (C) The numbers of times mice crossed to the correct platform (p = 0.0027) and the times in the correct quadrant (p = 0.0014) differed on the fifth day (n = 6). (D) Expression of CDC42 mRNA was associated with age through screening the GSE71620 dataset (ρ = -0.2537; P < 0.0001). Levels of Cdc42 mRNA were lower in older than younger mice (n = 6, p = 0.0255). Analysis of the GSE48350 dataset revealed levels of CDC42 mRNA are lower in the frontal cortex of people with AD than healthy subjects (p = 0.0406). (E) HDAC1 mRNA expression was associated with age through screening the GSE71620 dataset (ρ = 0.2347, p < 0.0001). Levels of Hdac1 mRNA were higher in older than younger mice (n=6, p = 0.0034). Analyzing the GSE48350 dataset revealed levels of HDAC1 mRNA are higher in the frontal cortices of people with AD than in those of healthy subjects (p = 0.0310). (F) Expression of YES1 mRNA was associated with age through screening the GSE71620 dataset (ρ = 0.1183; p = 0.0153). Levels of Yes1 mRNA were higher in older than younger mice (n = 6, p = 0.0357). Analyzing the GSE48350 dataset revealed levels of YES1 mRNA were higher in the frontal cortices of people with AD than in those of healthy subjects (p = 0.0011). (G) Functional link between YES1 and CDC42 determined using the STRING database to detect protein interactions. (H) YES1 mRNA levels associate with CDC42 in both neural aging and AD, as indicated by screening the GSE71620 (ρ = -0.1293; p = 0.0080) and GSE48350 (ρ = -0.3520; p < 0.0001) datasets.

Building the frontal cortex gene expression compendium

Gene expression data from microarray studies were downloaded from the Gene Expression Omnibus (GEO) using the terms (Frontal Lobes) OR (Frontal Cortex) OR (Anterior Central Gyrus) OR (Superior Frontal Gyrus) OR (Prefrontal Cortex) AND (aging) AND (transcriptional profiling). Microarray studies using RNA samples from human frontal cortices were included in our study. The basic characteristics used to identify the studies included first author, year of publication, dataset, platform, region, sample number and type. Because it has been shown that the changes in gene expression across cortical regions occur mainly when people are in their 20’s and 60’s [17,18], statistical group comparisons were made between subjects classified as younger (20-40 years) or older (60-90 years). Only microarrays that contained three or more samples in both the younger and older group were included in our study. Using these criteria, five studies were ultimately selected, all of which exhibited a common transcriptional signature in the prefrontal cortex throughout the aging. Postmortem human brain materials were obtained from subjects without a medical history of neuropathology, drug use, alcohol abuse, or psychiatric illness.

Microarray meta-analysis

We conducted a microarray meta-analysis using INMEX (a web-based tool for integrative meta-analysis of expression data) [45]. INMEX is in accordance with the Preferred Reporting Items for Systematic Reviews and Meta-Analyses guidelines for meta-analysis [46]. All gene probes were converted to a common Entrez ID using the gene/probe conversion tool in INMEX. After changing to Entrez ID, all datasets were preprocessed through a log2 transformation and Variance Stabilizing Normalization (VSN), followed by quantile normalization. Each individual dataset was visualized in box plots to ensure identical distribution among the samples. Differential expression analysis was performed independently for each dataset using INMEX with an FDR of 0.05 and a significance of p < 0.05. The moderated t test was based on the Limma algorithm. In INMEX, the results from individual microarray datasets are only for reference comparison and are not required for the subsequent steps of the meta-analysis [45]. For the meta-analysis, we used the Fisher’s method, a fixed effect model, and Vote counting (a significance level of p < 0.001, p < 0.001, votes number < 2). Fisher’s method (-2*∑Log (p)) is a statistical approach widely used in meta-analysis for combining P values from different studies independently of the sample size. This method is generally more sensitive than a combined analysis (i.e. it detects more DE genes). Fixed effect models combine effect sizes, and the estimated effect size in each study is assumed to come from an underlying true effect size plus measurement error. A fixed effect model can be selected based on statistical heterogeneity estimated using Cochran’s Q tests. Estimated Q values that approximate a chi-squared distribution suggest the fixed effect model assumption is appropriate. This method is usually conservative (fewer DE genes are detected, but with higher confidence). Vote counting is a meta-analysis method in which differentially expressed genes are first selected based on a threshold (a significance level of p < 0.05) to obtain a list of DE genes for each study. The vote for each gene can then be calculated as the total number of times it occurs in all DE lists. The final DE genes can be selected based on a minimal number of votes set by the user.

Pathway enrichment analysis

The pathways of the identified proteins were classified using the DAVID program (http://david.abcc.ncifcrf.gov) for KEGG annotation. To determine GO categories, we used Cytoscape - BinGO. A custom annotation file was created using the built-in annotation file for GO biological processes.

Network-based meta-analysis and extracting co-expressed PPI modules

A network-based meta-analysis was performed using NetworkAnalyst and STRING (Search Tool for the Retrieval of Interacting Genes/Proteins). Microarray datasets were processed as described above. Network construction was restricted to contain only the original seed proteins. Protein-protein interactions were predicted using the STRING database v10.5 (http://www.string-db.org/). Proteins that linked to each other were detected based on experimental determination, curated databases, gene neighborhood, gene fusions, gene co-occurrence, co-expression, text mining and protein homology. An extended network was constructed using a minimum required interaction score (> 0.9) as the selection parameter, which implies that only interactions with a high level of confidence were extracted from the database and considered as valid links for PPI networks.

Morris water maze

Wild-type C57/BL6 mice were purchased from the Experimental Animal Center at Tongji Hospital (Wuhan, China). All animals were housed at the animal care facility at Tongji hospital. These animal studies were approved by the Institutional Animal Research Committee of Tongji Hospital. The experimental protocols were approved by the Institutional Animal Care and Use Committee (Approved number: TJ-A20160503). Six-month-old mice comprised the younger group, while 20-month-old mice comprised the older group. The maze consisted of a circular tank (diameter, 150 cm) filled with water (27 °C) to which powdered milk was added. An escape platform (diameter, 10 cm) was located 1 cm beneath the water surface. Acclimation to the water maze was performed on day 0, after which learning trials were conducted on days 1-5. The open field test was performed on day 5. Mice participated in 3 trials per day for 5 consecutive days using a 60-s inter-trial interval.

Quantitative real-time polymerase chain reaction

The frontal cortex was collected from younger and older mice. Total RNA was extracted using a HI Pure RNA extract kit (Magen, China) and converted to cDNA using a ReverTra Ace qPCR RT kit (TOYOBO, Japan). Relative mRNA expression was detected using RT-PCR performed on an ABI Step One Plus (Applied Biosystems, USA) with SYBR green PCR master mix (TOYOBO, Japan). The primer sequences are as follows: GAPDH — forward, 5′ - AATGGTGAAGG TCGGTGT - 3′; reverse, 5′- GTGGAGTCATACTG GAACATGTAG - 3′; CDC42 — forward, 5′ - TGCT CTGCCCTCACACAGAAAG - 3′; reverse: 5′ - GCGG CTCTTCTTCGGTTCTG - 3′; and HDAC1 — forward, 5′ - CTCACCGAATCCGCATGACT - 3′; reverse, 5′ - GGCTTTGTGAGGACGGTAGA - 3′, YES1 — forward, 5′ - AATGAGGACCAGA GGGTAGGG - 3′; reverse, 5′- CATTATCAAATCCGCTCGCTCC - 3′. The sequences of other primers used in this study are shown in Table S7. The fold change in relative mRNA expression was calculated using the 2^-ΔΔCt method.

Statistical analysis

Network-based microarray meta-analysis was performed using INMEX and NetworkAnalyst. In the Morris water maze experiment, measures of performance during acquisition trials (i.e., escape latency) were averaged within each day for each animal. Data are reported as the mean ± standard error. To evaluate differences among the days, data were analyzed using two-way repeated measures ANOVA, with day as the within-subjects’ factor and treatment as the between-subjects’ factor. Differences between the two groups were analyzed using one-way ANOVA, after which Fisher’s LSD was used for post hoc comparisons. For RT-PCR experiments and the analysis of the GSE48350 dataset, statistical significance was determined using Student’s t-test (for independent or dependent samples, as appropriate). Values of p < 0.05 (two-tailed) were considered significant. Spearman correlation analysis was used to determine statistical significance for CDC42, HDAC1 and YES1 adjusting for age in the GSE71620 dataset. Other statistical analyses were performed using SPSS 19 software (SPSS Inc., USA).