Integrated analysis of differential gene expression profiles in hippocampi to identify candidate genes involved in Alzheimer's disease.

Alzheimer's disease (AD) is a complex neurodegenerative disorder with largely unknown genetic mechanisms. Identifying altered neuronal gene expression in AD may provide diagnostic or therapeutic targets for AD. The present study aimed to identify differentially expressed genes (DEGs) and their further association with other biological processes that regulate causative factors for AD. The present study performed an integrated analysis of publicly available gene expression omnibus datasets of AD hippocampi. Gene ontology (GO) enrichment analyses, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Protein‑Protein interaction (PPI) network analysis were performed. The present study detected 295 DEGs (109 upregulated and 186 downregulated genes) in hippocampi between AD and control samples by integrating four datasets of gene expression profiles of hippocampi of patients with AD. Respiratory electron transport chain (GO: 0022904; P=1.64x10‑11) was the most significantly enriched GO term among biological processes, while for molecular functions, the most significantly enriched GO term was that of protein binding (GO: 0005515; P=3.03x10‑29), and for cellular components, the most significantly enriched GO term was that of the cytoplasm (GO: 0005737; P=8.67x10‑33). The most significant pathway in the KEGG analysis was oxidative phosphorylation (P=1.61x10‑13). PPI network analysis showed that the significant hub proteins contained β-actin (degree, 268), hepatoma-derived growth factor (degree, 218) and WD repeat‑containing protein 82 (degree, 87). The integrated analysis performed in the present study serves as a basis for identifying novel drug targets to develop improved therapies and interventions for common and devastating neurological diseases such as AD.

Introduction

Alzheimer's disease (AD) is one of the most common and complex neurodegenerative disorders and is characterized by a progressive decline of memory and cognition (1). The disease is defined by specific neuropathological changes of neurofibrillary tangles (NFT) and amyloid plaques that accumulate in vulnerable brain regions (2,3). Neurodegeneration in the development of AD varies substantially across cell types and regions. Of note, it has been demonstrated that hippocampal CA1 pyramidal neurons are particularly vulnerable to neurodegeneration and bear NFTs during the early stages of AD (4,5); however, the underlying mechanisms of their degeneration have remained elusive.

AD is thought to be caused by the dysregulation of a large number of genes and the consequent alteration of their complex interactions, which finally contributes to the broad spectrum of disease phenotypes (6–9). Microarray technology, which provides researchers with a tool to assess the expression levels of thousands of genes simultaneously, offers the possibility of gaining insight into gene networks disturbed in intricate human disease such as AD, and to obtain possible molecular clues regarding the underlying mechanisms of the pathophysiology of AD. Previous studies have used this technique to more comprehensively enhance the knowledge of the cellular and molecular changes underlying AD (10–14). Although these studies have yielded significant novel insights, inconsistencies are present across these studies due to limitations based on small sample sizes and various results obtained by different groups with different laboratory protocols, microarray platforms and microarray data interpretations (15). In view of this, the present study integrated hippocampus gene expression datasets from multiple AD microarray studies to overcome these limitations of individual studies, resolve inconsistencies and provide significant novel insight into the complex biological processes involved in AD.

Materials and methods

Identification of eligible gene expression profiles of hippocampi of patients with AD

Hippocampal gene expression profiling studies in patients with AD were identified by searching the Gene Expression Omnibus database (GEO; http://www.ncbi.nlm.nih.gov/geo) (16). The following key words and their combinations were used: 'Alzheimer's disease', 'hippocampus', 'gene expression' and 'microarray.' Only experimental studies that had performed hippocampal gene expression profiling in patients with AD as well as normal control (NC) subjects were used. Non-human studies, review articles and integrated analyses of expression profiles were excluded.

Data preprocessing

Normalization is crucial for comparing different microarray datasets. The heterogeneity caused by different microarray platforms, gene nomenclature and clinical samples may make it difficult to compare the expression data directly. However, inappropriate normalization may contribute to the skewing of results and reduce their statistical significance. Consequently, a global normalization approach to minimize any inconsistencies should be included. For this propose, MATLAB Bioinformatics Toolbox was used in the present study to pre-process the raw microarray data of each study by Quantile normalization and log2 transformation to obtain intensity values.

Statistical analysis

MATLAB software, version 2013a (MathWorks, Natick, MA) was used to identify the differently expressed probe sets in the hippocampal tissues of patients with AD compared to those of NC subjects. A gene-specific t-test was performed, followed by calculation of the P-value and the effect size of the individual microarray study. Fisher's combined probability method was used to combine P-values from multiple studies, and the random effects model was used to combine effect sizes from multiple studies. Genes with an effect size >0.8 and a P-value <0.01 were selected as the significantly differentially expressed genes (DEGs).

Functional annotation of DEGs

To gain insight into the biological functions of DEGs, gene ontology (GO) classification was performed. GO provides a common descriptive framework as well as functional annotation and classification for analyzing the gene expression datasets. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) pathway enrichment analysis was performed to map the potential pathways of the DEGs. The KEGG pathway database is a recognized and comprehensive database, which includes an extensive variety of biochemical pathways (17). The online-based software GENECODIS, version 3 was utilized in the present analysis (18).

Protein-protein interactions (PPIs) network construction

PPI analysis allows for the assessment of protein functions at the molecular level, which are divided into the categories of cellular growth, development, metabolism, differentiation and apoptosis (19). The detection of key protein-interacting ions in the PPI networks of AD is important for the interpretation of cellular regulatory mechanisms in the development of the disease (20). The present study adopted the Search Tool for the Retrieval of Interacting Genes/Proteins (http://www.string-db.org/), a database of known and predicted protein interactions, to construct the PPI network and then visualized the distribution characteristics of the top 10 up- and downregulated DEGs in the network with Cytoscape software, version 3.2.0 (21).

Results

Identification of DEGs in hippocampi of patients with AD

The present study collected a total of four datasets of gene expression profiles in hippocampi of patients with AD according to the inclusion criteria; in total, data on the gene expression in 73 samples from patients with AD and 61 samples from control subjects were analyzed. The studies containing the individual hippocampal expression profiles in patients with AD are listed in Table I (12,22–24). A total of 11,494 genes from four expression profiling studies were assessed. For the purpose of global normalization, the raw microarray data were pre-processed by Quantile normalization and log2 transformation to obtain intensity values for each probe, which were used in the gene expression profiling. Subsequently, MATLAB software was utilized to identify DEGs in hippocampi between patients with AD and control subjects. Finally, a total of 295 DEGs were regarded as significantly differentially expressed between samples of patients with AD and NC subjects (109 upregulated and 186 downregulated genes) when the threshold was set as P<0.01 and effect size >0.8. A list of the top 10 most significantly up- or downregulated genes is presented in Table II. The pattern of expressional changes of the top 50 most significantly DEGs is displayed in a heat map in Fig. 1.

Table I

Characteristics of the individual studies.

GEO ID	Platform	Samples (n) (cases:controls)	Country	Year	Author
GSE29378	GPL6947 Illumina HumanHT-12V3.0 expression beadchip	31:32	USA	2013	Miller JA (22)
GSE36980	GPL6244 [HuGene-1_0-st]Affymetrix Human Gene 1.0 ST Array	7:10	Japan	2013	Hokama M (23)
GSE5281	GPL570 [HG-U133_Plus_2] AffymetrixHuman Genome U133 Plus 2.0 Array	13:10	USA	2007	Liang WS (24)
GSE1297	GPL96 [HG-U133A] AffymetrixHuman Genome U133A Array	22:9	USA	2004	Blalock EM (12)

GEO, Gene Expression Omnibus.

Table II

Top 10 most significantly up- or down-regulated differentially expressed genes.

Gene ID	Gene symbol	Official full name	P-value	Effect size
Upregulated genes
51663	ZFR	Zinc finger RNA binding protein	1.12×10−6	0.94853
2669	GEM	GTP binding protein overexpressed in skeletal muscle	1.21×10−6	1.1503
6277	S100A6	S100 calcium binding protein A6	2.74×10−6	1.0918
80335	WDR82	WD repeat domain 82	3.05×10−6	0.80325
5209	PFKFB3	6-Phosphofructo-2-kinase/fructose-2, 6-Biphosphatase 3	3.24×10−6	1.1908
3895	KTN1	Kinectin 1 (kinesin receptor)	3.36×10−6	1.1643
3068	HDGF	Hepatoma-derived growth factor (high-mobility group protein 1-like)	4.57×10−6	1.1580
2077	ERF	Ets2 repressor factor	5.16×10−6	1.0677
7049	TGFBR3	Transforming growth factor, beta receptor III	5.34×10−6	1.1818
5042	PABPC3	Poly(A) binding protein, cytoplasmic 3	7.09×10−6	1.0283
Downregulated genes
22820	COPG1	Coatomer protein complex, subunit gamma	8.69×10−8	1.3017
58189	WFDC1	WAP four-disulfide core domain 1	4.19×10−7	1.3844
10093	ARPC4	Tubulin tyrosine ligase-like family, member 3; actin related protein 2/3 complex, subunit 4, 20kDa	7.43×10−7	1.1458
60	ACTB	Actin, beta	1.08×10−6	1.1657
9158	FIBP	fibroblast growth factor (acidic) intracellular binding protein	1.15×10−6	1.1794
56993	TOMM22	Translocase of outer mitochondrial membrane 22 homolog (yeast)	2.35×10−6	1.0773
2537	IFI6	Interferon, alpha-inducible protein 6	2.57×10−6	1.1754
9556	C14orf2	Chromosome 14 open reading frame 2	2.72×10−6	1.1734
55837	EAPP	E2F-associated phosphoprotein	3.11×10−6	1.0270
5889	RAD51C	RAD51 homolog C (S. cervisiae)	5.16×10−6	1.0461

Figure 1

Heat map visualization of the patterns of expressional changes for the top 50 most significantly differentially expressed genes across various datasets.

The upregulated gene with the lowest P-value was ZFR, which is mainly expressed in neural tissue, but also weakly expressed in other tissue types (25,26), suggesting a neuronal function. A recent study identified ZFR as a putative genes associated with hereditary spastic paraplegias by using whole-exome sequencing (27). The downregulated gene with the lowest P-value was COPG1, whose function has yet to be determined.

Functional annotation

To investigate the biological roles of the DEGs in the hippocampi of patients with AD, the present study performed a categorized GO enrichment analysis. GO provides a common descriptive framework and functional annotation of the gene datasets. GO categories are separated into three groups: Biological processes, cellular components and molecular function. The present study examined GO categories separately using the web-based software GENECODIS. The results showed that genes associated with the respiratory electron transport chain (GO: 0022904; P=1.64×10−11) and gluconeogenesis (GO: 0006094; P=2.84×10−5) were significantly enriched among biological processes, while for molecular functions, protein binding (GO:0005515; P=3.03×10−29) and nucleotide binding (GO: 0000166; P=5.41×10−14) were significantly enriched, and with regard to cellular components, genes associated with the cytoplasm (GO:0005737; P=8.67×10−33) and mitochondrion (GO: 0005739; P=1.00×10−23) were significantly enriched (Table III, Fig. 2A).

Table III

GO terms of differentially expressed genes (top 15).

GO ID	GO term	No. of genes	FDR
Biological processes
GO:0022904	Respiratory electron transport chain	14	1.64×10−11
GO:0006094	Gluconeogenesis	7	2.84×10−5
GO:0006006	Glucose metabolic process	9	4.94×10−5
GO:0006200	ATP catabolic process	8	1.15×10−4
GO:0010388	Cullin deneddylation	4	1.18×10−4
GO:0006915	Apoptotic process	19	1.35×10−4
GO:0006096	Glycolysis	6	2.91×10−4
GO:0007264	Small GTPase mediated signal transduction	13	3.44×10−4
GO:0042776	Mitochondrial ATP synthesis coupled proton transport	4	3.61×10−4
GO:0015992	Proton transport	6	3.84×10−4
GO:0007165	Signal transduction	26	7.19×10−4
GO:0006120	Mitochondrial electron transport, NADH to ubiquinone	5	1.01×10−3
GO:0006810	Transport	17	1.26×10−3
GO:0048146	Positive regulation of fibroblast proliferation	5	1.30×10−3
GO:0016071	mRNA metabolic process	10	1.40×10−3
Molecular function
GO:0005515	Protein binding	117	3.03×10−29
GO:0000166	Nucleotide binding	59	5.41×10−14
GO:0005524	ATP binding	37	5.48×10−7
GO:0005525	GTP binding	16	1.32×10−5
GO:0046961	Proton-transporting ATPase activity, rotational mechanism	5	4.49×10−5
GO:0003924	GTPase activity	12	4.55×10−5
GO:0015631	Tubulin binding	5	8.56×10−5
GO:0046933	Hydrogen ion transporting ATP synthase activity, rotational mechanism	4	3.08×10−4
GO:0005509	Calcium ion binding	18	5.87×10−4
GO:0008137	NADH dehydrogenase (ubiquinone) activity	5	6.24×10−4
GO:0003713	Transcription coactivator activity	10	6.81×10−4
GO:0022857	Transmembrane transporter activity	5	8.83×10−4
GO:0005516	Calmodulin binding	8	1.17×10−3
GO:0047485	Protein N-terminus binding	6	1.81×10−3
GO:0003878	ATP citrate synthase activity	2	1.87×10−3
Cellular components
GO:0005737	Cytoplasm	133	8.67×10−23
GO:0005739	Mitochondrion	60	1.00×10−23
GO:0005829	Cytosol	62	6.03×10−16
GO:0005743	Mitochondrial inner membrane	25	1.38×10−15
GO:0005634	Nucleus	101	4.71×10−14
GO:0005856	Cytoskeleton	26	7.29×10−7
GO:0005625	Soluble fraction	17	1.64×10−6
GO:0016020	Membrane	66	3.04×10−6
GO:0005886	Plasma membrane	60	3.94×10−6
GO:0005759	Mitochondrial matrix	12	4.18×10−6
GO:0005654	Nucleoplasm	25	4.28×10−6
GO:0045121	Membrane raft	9	3.28×10−5
GO:0005730	Nucleolus	31	4.98×10−5
GO:0005753	Mitochondrial proton-transporting ATP synthase complex	4	1.58×10−4
GO:0030054	Cell junction	15	2.50×10−4

FDR, false discovery rate; GO, gene ontology.

Figure 2

Functional annotation of significantly enriched differentially expressed genes. (A) The top 10 enriched gene ontology categories for biological processes; (B) The top 10 enriched Kyoto Encyclopedia of Genes and Genomes pathways.

The present study subsequently performed a KEGG pathway enrichment analysis in order to further evaluate the biological roles of the DEGs. A hypergeometric test with P<0.05 was used as the criterion for pathway detection. According to the KEGG analysis, oxidative phosphorylation was the most significant pathway (P=1.61×10−13). Furthermore, pathways involved in Parkinson's disease (P=3.77×10−10) and Huntington's disease (P=2.26×10−8) were also highly enriched (Table IV, Fig. 2B).

Table IV

KEGG pathways of differentially expressed genes (top 15).

KEGG ID	KEGG term	No. of genes	FDR	Genes
hsa00190	Oxidative phosphorylation	17	1.61×10−13	COX4I1, NDUFAB1, UQCRC1, NDUFV2, COX6A1, COX6B1, ATP5L, ATP5J2, ATP6AP1, ATP6V1E1, NDUFA4, ATP5C1, NDUFA9, ATP6V1B2, NDUFS3, UQCRH, ATP5B
hsa05012	Parkinson's disease	14	3.77×10−10	COX4I1, NDUFAB1, VDAC2, UQCRC1, NDUFV2, COX6A1, COX6B1, NDUFA4, ATP5C1, NDUFA9, UCHL1, NDUFS3, UQCRH, ATP5B
hsa05016	Huntington's disease	14	2.26×10−8	COX4I1, NDUFAB1, VDAC2, UQCRC1, NDUFV2, COX6A1, COX6B1, EP300, NDUFA4, ATP5C1, NDUFA9, NDUFS3, UQCRH, ATP5B
hsa05010	Alzheimer's disease	12	5.59×10−7	COX4I1, NDUFAB1, UQCRC1, NDUFV2, COX6A1, COX6B1, NDUFA4, ATP5C1, NDUFA9, NDUFS3, UQCRH, ATP5B
hsa04310	Wnt signaling pathway	10	1.70×10−5	TCF7L1, PORCN, EP300, PPP2R1A, JUN, CTNNBIP1, NFAT5, TBL1X, RBX1, DAAM1
hsa05131	Shigellosis	7	1.97×10−5	ARPC4, NFKBIA, ELMO1, ACTB, ARPC1A, WASF2, CD44
hsa04520	Adherens junction	7	4.80×10−5	TCF7L1, ACTB, EP300, EGFR, WASF2, FYN, SORBS1
hsa04260	Cardiac muscle contraction	7	6.66×10−5	COX4I1, UQCRC1, COX6A1, COX6B1, ATP1A1, SLC9A6, UQCRH
hsa00020	Citrate cycle (TCA cycle)	5	7.87×10−5	ACLY, OGDHL, MDH2, ACO2, SUCLG1
hsa05120	Epithelial cell signaling in Helicobacter pylori infection	6	2.87×10−4	NFKBIA, ATP6AP1, ATP6V1E1, EGFR, JUN, ATP6V1B2
hsa03050	Proteasome	5	3.05×10−4	PSMD13, PSMB3, PSMC1, PSMC3, PSMD8
hsa05130	Pathogenic Escherichia coli infection	5	1.11×10−3	ARPC4, ACTB, ARPC1A, TUBA1B, FYN
hsa00250	Alanine, aspartate and glutamate metabolism	4	1.64×10−3	NIT2, GAD1, ADSL, GOT1
hsa05100	Bacterial invasion of epithelial cells	5	3.27×10−3	ARPC4, ELMO1, ACTB, ARPC1A, WASF2
hsa04810	Regulation of actin cytoskeleton	8	4.03×10−3	ARPC4, SSH3, MYL12B, ACTB, ARPC1A, EGFR, WASF2, ITGB8

KEGG, Kyoto Encyclopedia of genes and genomes; FDR, false discovery rate; hsa, Homo sapiens.

PPI network construction

The present study established the PPI networks of the top 10 upregulated and downregulated DEGs using Cytoscape software. The interaction network included 863 nodes and 1,304 edges. In the PPI network, degrees of interaction were defined to determine the number of neighbors a node directly connected to, and nodes with a high degree of interaction were defined as hub proteins. The significant hub proteins included β-actin (ACTB; degree, 268), hepatoma-derived growth factor (HDGF; degree, 218) and WD repeat-containing protein 82 (WDR82; degree, 87) (Fig. 3).

Figure 3

Constructed protein-protein interaction networks of the top 10 up- and down-regulated DEGs. Nodes represent proteins and edges represent interactions between two proteins. Red- and green-colored nodes represent products of up- and down-regulated DEGs, respectively. Blue nodes denote products of genes predicted to interact with the DEGs. DEG, differentially-expressed gene.

Discussion

The present study aimed to identify altered hippocampus gene expression and their further association with other biological processes that regulate causative factors for AD to provide diagnostic factors or therapeutic targets of AD. An integrated analysis of DEGs from four publicly available GEO datasets of hippocampi from patients with AD was performed. In total, 295 genes were consistently differentially expressed across the studies with 109 upregulated genes and 186 downregulated genes. The upregulated gene with the lowest P-value was ZFR, which is mainly expressed in neural tissue (25,26) and may therefore have a role in neuronal function. A recent study identified ZFR as a putative gene associated with hereditary spastic paraplegias by using whole-exome sequencing (27), and from this, the present study deduced that ZFR may be implicated in the underlying processes of AD, which is required to be confirmed by further experiments. The downregulated gene with the lowest P-value was COPG1, whose function remains to be elucidated.

In line with previous studies, certain genes identified in the present study have been closely associated with the development of AD, including S100A6 and TGFBR3. A study on the roles of S100 family proteins in nervous system function and disease found that mRNA expression levels of six family members (S100A1, S100B, S100A6, S100A10, S100A4, S100A13) displayed a 100-fold range in mouse brains, five of which (S1100A1, S100A6, S100A10, S100A13, and S100B) showed age-dependent increases in adult mice that ranged from 5- to 20-fold (28). S100A6-protein immunoreactivity was found to be specifically located within astrocytes associated to amyloid plaques in an APP/London transgenic mouse model of AD, as well as in the brains of patients with AD. S100A6 was upregulated in the amygdala as well as in hippocampal regions (29). Another study detected that biglycan proteoglycans were upregulated in familial AD, while TGFBR3 was markedly downregulated in sporadic AD fibroblasts. Furthermore, the differential expression of TGFBR3 in familial AD and sporadic AD cells was associated with the severity of AD (30).

In the present study,the results of the PPI network analysis of the top 10 upregulated and downregulated DEGs indicated that the significant hub proteins included ACTB, HDGF and WDR82. ACTB, which encodes β-actin, is a candidate reference gene for normalization of target gene expression in polymerase chain reaction (PCR) analysis due to its high conservation. A previous study determined the mRNA levels of ACTB and other genes in the frontal cortex of patients with AD and control subjects using PCR analysis with SYBR Green technology to identify suitable endogenous reference genes in human post-mortem brain tissues for the expression analysis of potential candidate genes associated with AD (31); according to this study, ACTB was the least suitable candidate with reliable expression among a set of suitable endogenous reference genes due to low expression stability in the frontal cortex of AD (32). Of note, the actin cytoskeleton has been reported to have an important role in AD pathology by mediating synaptic degeneration (32).

In order to elucidate the biological roles of the DEGs in AD, a categorized GO enrichment analysis was performed in the present study. The results showed that the respiratory electron transport chain was the most significantly enriched GO category for biological processes. To further evaluate the biological role for the DEGs, the present study performed a KEGG pathway enrichment analysis. According to the KEGG analysis, the most significantly enriched pathway was oxidative phosphorylation. A previous study provided evidence of neuronal metabolic impairments at the transcriptomic and protein level in the brains of patients with AD (33), which was ascribed to the downregulation of mitochondria-associated genes, in particular, oxidative phosphorylation genes in consistency with the fact that AD is a degenerative disease. Furthermore, the present study found that the pathways of several neurodegenerative diseases, including Parkinson's disease, Huntington's disease and Alzheimer's disease, were also highly enriched according to the KEGG pathway enrichment analysis, which was due to dysregulation of genes associated with mitochondrial energy metabolism, including COX4I1, NDUFAB1, UQCRC1, NDUFV2, COX6A1, COX6B1. This finding validated the integrated analysis methods used in the present study.

It is noteworthy that the present study had several limitations. The heterogeneity of the datasets used may have distorted the analysis, as clinical samples may have been heterogeneous with regard to clinical activity or gender. Furthermore, the effects of varying degrees of severity of AD on the differences in hippocampal gene expression were not taken into account. However, the present integrated analysis of different datasets of hippocampal gene expression in patients with AD may have facilitated the detection of genes that would have been missed in the analysis of a single patient or study cohort. Despite these limitations, the present study provided novel information regarding the molecular mechanisms of AD; however, further analyses are required to confirm the present findings.

In conclusion, the present study performed an integrated analysis, which provided significant insight into the global molecular changes associated with AD pathology. Furthermore, the present study identified DEGs as well as other biological functions, which may contribute to the successful identification of diagnostic factors or therapeutic targets for AD and the development of effective targeted therapies. Further functional studies may provide additional insight into the role of the DEGs in the pathophysiology of AD.