Uncovering Molecular Biomarkers That Correlate Cognitive Decline with the Changes of Hippocampus' Gene Expression Profiles in Alzheimer's Disease

Background Alzheimer's disease (AD) is characterized by a neurodegenerative progression that alters cognition. On a phenotypical level, cognition is evaluated by means of the MiniMental State Examination (MMSE) and the post-morten examination of Neurofibrillary Tangle count (NFT) helps to confirm an AD diagnostic. The MMSE evaluates different aspects of cognition including orientation, short-term memory (retention and recall), attention and language. As there is a normal cognitive decline with aging, and death is the final state on which NFT can be counted, the identification of brain gene expression biomarkers from these phenotypical measures has been elusive. Methodology/Principal Findings We have reanalysed a microarray dataset contributed in 2004 by Blalock et al. of 31 samples corresponding to hippocampus gene expression from 22 AD subjects of varying degree of severity and 9 controls. Instead of only relying on correlations of gene expression with the associated MMSE and NFT measures, and by using modern bioinformatics methods based on information theory and combinatorial optimization, we uncovered a 1,372-probe gene expression signature that presents a high-consensus with established markers of progression in AD. The signature reveals alterations in calcium, insulin, phosphatidylinositol and wnt-signalling. Among the most correlated gene probes with AD severity we found those linked to synaptic function, neurofilament bundle assembly and neuronal plasticity. Conclusions/Significance A transcription factors analysis of 1,372-probe signature reveals significant associations with the EGR/KROX family of proteins, MAZ, and E2F1. The gene homologous of EGR1, zif268, Egr-1 or Zenk, together with other members of the EGR family, are consolidating a key role in the neuronal plasticity in the brain. These results indicate a degree of commonality between putative genes involved in AD and prion-induced neurodegenerative processes that warrants further investigation.


Introduction
Gomez Ravetti and Moscato have recently shown that the abundance of five proteins, within a panel that also measured other 115 cytokines and growth factors, can be used to predict the development of clinical Alzheimer's Disease (AD) [1]. The biomarker molecular signature is composed of IL-1a, TNF-a, IL-3, EGF and G-CSF and has the same level of specificity and sensitivity as the original 18-protein signature proposed by Ray et al. [2] in late 2007, who introduced this important dataset in the literature. In the original work, Ray et al. had employed the abundance of 120 signalling proteins in plasma to obtain their 18protein signature set. They used a training set of 83 samples to identify patients that progressed to AD in two to six years. The proposed 5-protein signature has an average of 96% accuracy in predicting clinical AD but it is still linked to the joint measurement of 120 protein abundances.
In this paper, we are revisiting the quest of finding biomarkers of AD. However, this time we aim at finding biomarkers in hippocampus tissue samples which would complement the results of the previous studies on plasma biomarkers. This study will now give a different perspective on the progression of the disease, keeping a systems biology and functional genomics approach. Towards this end, we have chosen to rely on an informative experimental design and dataset contributed by Blalock et al. [3]. We believe that their dataset may help us to locate, either directly or indirectly, other biomarkers of interest that could eventually be detectable in plasma.
Blalock et al. analysed samples from 35 patients with four different levels of AD severity: control, incipient, moderate and severe; for this paper we used only 31 samples for which information is available online. The label assigned to each sample (its ''level of severity'') was decided after considering two important scores, those provided by the MiniMental State Examination (MMSE) and the Neurofibrillary Tangle count (NFT). The MMSE score is based on a questionnaire that aims at measuring the level of cognitive impairment of a patient. The questions are aimed at evaluating different aspects of cognition, such as orientation, shortterm memory (retention and recall), attention and language. A normal score can range from 24 to 30, mild cognitive impairment on the interval 20 to 23, moderate AD between 10 to 19, and the rest (from 0 to 9) are all considered severe AD cases.
As previously mentioned, Blalock et al. [3] also used the NFT score to assign a severity label to each sample. The NFT score is a well established method for the neuropathological diagnosis of AD [4]. The score is usually based on the average counts of neurofibrilary tangles considering different regions of the brain. A NFT score is a recognised indicator of AD, nevertheless, it is not completely effective as there is evidence that NFTs were also identified in healthy aging brains [5,6,7,8].
The analysis by Blalock et al. [3] focused on the identification of ADrelated genes (ADG) and incipient ADG (IADG) using a methodology based on the correlation of the genes with NFT and MMSE scores. In turn, they identified putative biological processes and signalling pathways that are significantly present in those gene lists. Our analysis takes a different direction. While still based on the same dataset, we are attempting to map the progression of the disease, finding biomarkers linked to disease severity, by identifying the genes associated with the divergence of the gene expression profile of a sample with the gene expression average profile of the Figure 1. This plot illustrates that the third step of our methodology, the use of the Jensen-Shannon divergence, does not appear to give an interesting separation of the samples in the absence of a previous feature selection step. For this graph, all 22,215 genes were considered in the calculation of the average profile of the samples in the ''Control'' and ''Severe AD'' classes. The square root of the Jensen-Shannon divergences to the ''Control'' and ''Severe AD'' average profile are computed, respectively giving, for each sample, its x and y coordinates in this plot. Observe that most of the ''Control'' samples have values lower than 0.12, with two exceptions. This result is expected, as the probability distribution function of the ''Control'' class was used. However, most of the samples from AD patients (having either ''Incipient AD'', ''Moderate'' or ''Severe'' labels), show a divergence with the Control average gene expression profile. Figure 2 shows the important contribution provided by the feature selection step. doi:10.1371/journal.pone.0010153.g001 Figure 2. This plot illustrates that after application of the feature selection steps, followed by the computation of the gene expression profile's average profile of the samples in the ''Control'' and ''Severe AD'' classes (now on a set of 1,372 probes), the samples are now more clearly separated. Here, all ''Control'' samples have the square root of the Jensen-Shannon divergences to the average gene expression of the ''Control'' samples (x-coordinate) smaller than 0.12 (almost all severe AD have x-coordinates greater than 0.15). In addition to that, most samples labelled ''Severe AD'' are located on the same region. Both results are expected. However, it is interesting that in this (x,y)-plot most samples that are labelled ''Incipient AD'' or ''Moderate AD'' seem to ''bridge'' between the regions that have most of the ''Control'' samples and the region that have most of the ''Severe AD'' group. This result is interesting as no samples from ''Incipient AD'' nor ''Moderate AD'' have been used in the first three steps of our methodology. In essence, the work is a ''test set'' indicating that it is reasonable to expect that some genes in the genetic signature of 1,372 probes have information about a putative ''progression'' trend of the disease, from the ''Control'' to the ''Severe AD'' profile. In what follows, correlations across all the samples with these divergences are used as a method to try to identify those gene profiles that are most correlated with the progression from ''Control'' to ''Severe AD''. doi: 10.1371/journal.pone.0010153.g002 ''Control'' group. Analogously, we are interested in identifying the genes that seem to best correlate with the ''convergence'' to the average profile of the ''AD Severe'' group of samples. The difference between Blalock et al.'s [3] methodological approach to data analysis and ours is very important. We aim to uncover genes that correlate with the divergence of the gene expression profiles, instead of relying only on correlations with the NFT and MMSE values.
Our objective is to uncover genes which are highly correlated to the progression of the disease. With this objective in mind, we will concentrate the first part of our analysis on the two most extremely separated classes, the sets of samples that have been labelled as ''Control'' and those labelled ''AD Severe''. This important initial decision was made based on the fact that the four classes are, in some sense, arbitrarily defined as specific thresholds for the MMSE and NFT scores that were decided ad hoc. Therefore, we decided to first focus on the transitional patterns that can be identified from a ''normally aging'' to an ''AD-severe'' gene expression profile in hippocampus. With this approach, we also avoid selecting genes that diverge from the normal-aged profile by causes other than AD, as we expect that the severity scale in AD has a higher probability of being correct in the ''Severe AD'' cases (since they have high values of NFT and low MMSE scores, clearly a joint combination highly appreciated as a disease hallmark). This approach has an additional advantage. Using this particular dataset and with focus on the effects of incorrect diagnoses, two publications indentify four possible misdiagnoses between control and incipient AD [9,10]. In our case, the samples that have been labelled either ''Incipient AD'' or ''Moderate AD'' play the role of a ''Test set'', as they are not used to select probes for establishing a molecular signature, thus avoiding misdiagnoses problems.

Results
The results have been obtained using four steps in tandem: 1) abundance quantization of gene expression values and filtering of probes (this step is supervised by using the samples labelled either ''Control'' or ''Severe AD''); 2) a feature selection algorithm to refine the probe selection based on numerical solution of a combinatorial optimization problem (the (alpha,beta)-k-Feature Set methodology); 3) a correlation analysis (that requires the computation of Jensen-Shannon divergences). Finally, a fourth step involves the pathway and Gene Ontology analysis of the results.
The first two steps only used the samples labelled either ''Control'' or ''Severe AD''. The third step requires several procedures and uses all of the samples. We first compute an   average gene expression profile for the classes ''Control'' and ''Severe AD''. This step is followed by the computation of the square root of the Jensen-Shannon divergence [11] of the gene expression profile of each sample with the average profiles of the classes ''Control'' and ''Severe AD''. Finally, we perform a correlation analysis of each gene expression profile (now across all samples) with the results of the square root of the Jensen-Shannon divergence (we do it twice, one for the ''Control'' and the other for the ''Severe AD'' case). With this information, and using stateof-the-art pathway analysis and text mining tools, as a result of our final analysis step, we provide a comprehensive list of results of the differentially regulated genes, patterns of up (down)-regulation and the pathways that seem to be implicated in the progression of AD. We refer to the Methods section for a completely reproducible and in-depth explanation of our methodology.

Probe selection and Jensen-Shannon divergence computations based on class information
We start our analysis with a baseline comparison, which we have chosen to include for illustrative purposes. Figure 1 provides an example of the importance of performing an initial probe/gene selection step. The example serves as an argument for the necessity of the first two steps of our method. We have normalized each individual gene expression profile, and we have computed the average gene expression profile for classes ''Control'' and ''Severe AD'' (following the same procedure we will use in the third step of our method, but in this case using all probes in the array).
We have used the square root of the Jensen-Shannon divergence of a pair of samples (a pair of gene expression profiles) as our measure of ''dissimilarity'' between them. The square root of the Jensen-Shannon divergence quantifies the difference between two probability distribution functions (PDFs) and it is a metric (we refer the reader to the Methods section for a mathematical definition and a discussion of its properties). Figure 1 plots the divergence of each sample with the average expression profile of the classes 'Control' and 'Severe AD'; sqrtJSD(P, P P C ) denotes the square root of the Jensen-Shannon divergence between sample P and the average profile on the 'Control' class ( P P C ). Analogously, sqrtJSD(P , P P S ) denotes the square root of the Jensen-Shannon divergence between sample P and the average profile on the 'Severe AD' class ( P P S ). The advantage of using the probe/gene selection steps, which reduces the number of genes to the most informative ones, will be evident when we later compare Figure 1 with Figure 2. However, Figure 1 already shows some interesting patterns. For instance, we can observe that a high percentage of the samples from AD patients (having either 'Incipient AD', 'Moderate' or 'Severe' labels) show sqrtJSD(P, P P C ) values greater than 0.115, which indicates measurable divergence with the Control average gene expression profile. Figure 2 presents the same procedure, but only after the feature selection step has significantly reduced the number of probes fom 22,215 to 1,372. We refer to the Methods section for details. In Figure 2, an arguably more coherent arrangement can be observed. As expected, the group of control samples (in green) have lower values of sqrtJSD(P, P P C ) and higher values of sqrtJSD (P, P P S ). Obviously, the opposite behaviour is observed for the samples belonging to the severe cases. What cannot be expected, however, is a layout of the samples that could provide evidence of a continuous ''progression'' of the disease. The Figure shows that the samples of 'Incipient AD' are close to the control group and the 'Moderate AD' samples are closer to them and also link to severe AD. A priori, since those samples had not been used for probe selection, they could have been in any position in the (sqrtJSD(P , P P C ), sqrtJSD(P, P P S ) plane. Finally, Figure 3 presents the results of the MMSE score as a function of the sqrtJSD(P , P P C ), showing an inverse correlation between them. A similar situation happens between MMSE and sqrtJSD(P, P P S ), but in this case low MMSE scores correspond to low values of sqrtJSD(P,Ps), giving a positive correlation. It is this interplay between positive and negative correlations that has enabled us to find interesting biomarkers. In the next subsection, we explain how these correlations were used to identify probes that ''diverge from'' their values in the ''Control'' group and ''converge to'' the values in the ''Severe AD'' group.

Gene correlation analysis
The third step employs a correlation analysis to select the group of probes that are the most strongly correlated. Intuitively, the idea is fairly straightforward as illustrated in the following ''Gedankenexperiment'' (a thought experiment). Assume, for argument's sake, that the MMSE of each patient P is not actually phenotypical information assigned to each sample. Instead, assume that the MMSE values are the microarray probe expression of some gene. In this ''thought experiment'', let MMSE(P) be the expression of this hypothetical gene probe on sample P, and fDataset be the set of values it has for each sample. The correlation of the sampleordered set of values {MMSE(P)} with the set of sample-ordered values {sqrtJSD(P, P P C )} is negative, indicating that, in general, this hypothetical MMSE probe reduces its values as the whole gene expression profile of sample P diverges from the average ''Control'' profile ( Figure 3). Analogously, there exists a positive correlation of the set of values {MMSE(P)} with the values of the set {sqrtJSD(P, P P S )}. This indicates that the values of MMSE tend to be reduced as the profile of sample P ''converges to'' the average profile of samples in the ''Severe AD'' group. We have computed these correlations for all probes in the signature, which are given in the supplementary material (File S2 sheet 'correlation Analysis') and are the basis for our analysis.
We also refer the reader to Figure 4, which presents the computed correlations. Tables 1 and 2 present the one hundred most correlated probes (in absolute values). In the supplementary material (File S2 sheet 'correlation Analysis'), the correlation of each of the 1,372 probes that were selected by our method is given (and annotated, including Affymetrix and Stanford's Source outputs) to facilitate further analyses.
As the objective is to detect the probes correlated with the progression of AD, we will select those probes with high absolute correlations values with both groups, an indication of a   divergence of the average control profile together with a convergence to the severe AD profile; these correlations computed over all sample types. We need to check both groups according to their correlations to the average profile. The first group of probes we are interested in are those that have a positive correlation with the sqrtJSD(P, P P C ) and a negative correlation with sqrtJSD(P, P P S ). The probes in this group are those probes with under-expression in the non-disease sample but are overexpressed in the severe AD cases. The second group has the opposite behaviour, the probes' expression values have a negative correlation with sqrtJSD(P, P P C ) and a positive correlation with sqrtJSD(P, P P S ). This pattern can be visualised in Figure 4, where the elliptical shape of the dispersion of the probes in this scatter plot indicates that our methodology has preserved all the significant probes for both classes and that there are no probes (after the filter) presenting a high correlation simultaneously with the control and severe AD profiles.
On these values a new selection criterion is applied, as we wanted to identify the group of probes that have strong correlations to both groups in absolute value. This symmetry of our argument stems from the interest in understanding the biology of the progression of the disease. For identifying disease biomarkers we may just concentrate in finding the probes that present an upregulation trend when progressing from ''Control'' to ''Disease''. However, here we would also like to identify those probes that become increasingly downregulated, which, in turn, would help us to identify significantly dysregulated biological pathways (as members of the pathway will be either up or downregulated). Towards this end, we rank the probes in the order given by their Euclidean distance from the origin of coordinates in Figure 4. We selected an arbitrary cut-off value of fifty probes (the selected probes are marked in red). These fifty probes are also identified by their Gene Symbols in Figures 5 and 6.
Calculating the distance of each probe to the origin, on the sqrtJSD system of coordinates, we further selected the 50 most distant probes and analysed their behaviour.   material presents the distance to the origin of the 1,372 probes analysed. In Table 3, it can be seen which genes have some putative annotation that links them to AD (17 genes out of 48). Figure 7 shows the heat map of the 50-probe signature, where the probes and patient samples are ordered by considering the similarity of their gene-expression values only. It can be observed that the Memetic Algorithm (MA), a high performance combinatorial optimization ordering method [12] for microarray datasets introduced in 2007, ordered most of the patients with or without an incipient level of AD on the left and the more severe cases on the right. When ordering the probes' gene expression, the MA perfectly sorted the groups previously described. We refer to [12,13] for details of the MA. The supplementary material (File S2 '1372 norm. +heat map+GO') presents the heat map of the 1,372 gene-probes, with samples and probes sorted by the MA.

Transcription factors analysis of most correlated probes
The signature of 50 probes we present in Figure 7 has 48 different genes (some probes are related to the same gene). The two repeated genes in this 50-probe list are ATP5C1 (ATP synthase, H+ transporting, mitochondrial F1 complex, gamma polypeptide 1) and PPIA (peptidylprolyl isomerase A (cyclophilin A)) [14,15,16,17], a calcineurin regulatory protein. A recent study that used RT-PCR to examine tissue from 90 AD and 81 control human brains reports that cyclophilin is reduced in AD (both for females and males as compared with their gender-matched groups) [18]. We note here that the cutoff of 50 probes circumscribes the initial description a little, but most of the later discussion uses information from the whole signature to identify dysregulated pathways. Figure 8 presents the heat map of the 1,372-probe signature. The probes were sorted with the MA but the samples remain in the same position as obtained previously with the 50probe signature.
We analysed this list of genes using GATHER [19], an online tool for annotating signatures. Forty-one genes out of fifty have a motif for EVI1 (ecotropic viral integration site 1) and thirty-nine of them have a binding motif with V$TCF1P_Q6 (TCF1: transcription factor 1, hepatic; LF-B1, hepatic nuclear factor (HNF1), albumin proximal factor). The same analysis can be done if we divide the set of genes in two groups. The first group has positive correlation with the control profile and are overexpressed in AD; the second group has a positive correlation with the severe profile, and tend towards being underexpressed in AD (see Table 3). Table 4 presents the overrepresented motifs. We note, however, that we believe that the best results to identify putative overrepresented regulatory motifs can be obtained using the whole signature of 1,372 probes, and we will present the results of this investigation after presenting the case of the most correlated probes.
Another interesting pattern emerged when analysing the KEGG Pathways of the 50-probe signature using GATHER and PATHWAY Studio [20]. Using GATHER, three KEGG Pathways appear significantly represented, Amyotrophic lateral sclerosis (ALS), Oxidative phosphorylation and ATP synthesis. Using PATHWAY Studio, we automatically built the ''commonregulators'' diagram by selecting a filter that only considers protein interactions and binding. The resulting diagram is presented in Figure 9. As can be seen from the figure, we have chosen a circular  membrane layout and our previously uncovered 5-protein signature [1] (IL1-a, TNF-a, IL-3, EGF and G-CSF) in plasma (plus IL-6) appears to have a strong relationship with CSF1 (colony stimulating factor 1 (macrophage)), the most positive correlated gene with the control profile (see Table 1). It is also worth mentioning, that CSF1 was found differentially expressed in blood of AD and Control subjects and belongs to the 18-protein signature uncovered by Ray et al. [2] in 2007.
Five of the 50 most correlated probes correspond to genes already mapped to KEGGs Alzheimer's disease Pathway KEGG:05010 and together with LDHA they link to impaired metabolism and the ''novel glucocorticoid hypothesis'' We have observed that five genes, which are the most correlated probes with our putative signature for disease severity, can be mapped to the AD pathway of the , and PPP3CA (protein phosphatase 3 (formerly 2B), catalytic subunit, alpha isoform), the last one also known as Calmodulin-dependent calcineurin A subunit alpha isoform. In all cases, the probes showed a reduction of expression with AD severity, which may indicate a sign of impaired mitochondrial functions and energy uptake [35,36].
In addition to these five, we observed the reduced expression of the glycolytic enzyme LDHA, which may also indicate another challenge for energy metabolism in these neurons. Although glucose is generally considered to be the only substrate for brain energy metabolism, moncarboxylates have also been hypotheised as alternative substrates [37]. Laughton et al. report segregation in the hippocampus, with LDHA present in astrocytes and not in neurons. Instead, it is pyruvate dehydrogenase that is present in neurons but not in astrocytes and as a consequence of this study they support the argument that a metabolic compartmentalization exists in the human cortex and hippocampus where lactate produced by astrocytes could be oxidized by neurons [37]. We have also observed a reduction in expression of a probe that corresponds to PDHA1 (Pyruvate dehydrogenase (lipoamide) alpha 1, 200980_s_at) with increasing AD severity. The reduction of PDH expression, and the concurrent increase in pyruvate carboxylase gene expression, was discussed by Landfield et al.
[38], who argue that: ''These changes suggest that reduced pyruvate flux through PDH and decreased oxidative metabolism of glucose may develop early in AD. Interestingly, the inactivation of PDH is also a major pathway through which glucocorticoid activity acts to conserve glucose, and apparently, to induce insulin resistance [65,66]. Thus, our data are consistent with the possibility that GC effects on this and other important target pathways in brain are enhanced in both aging and AD. If so, such alterations in glucocorticoid efficacy may have implications for AD pathogenesis as well as for the increased risk of AD associated with normal aging.'' Our results seem to indicate that LDHA might also be discussed within the extended metabolic pathways that serve as the basic framework of this novel, more complex hypothesis [   . In comparison to SV2B this new mRNA encodes for the same protein but it has an elongated 39-untranslated region (39UTR) that contains several AU-rich (AUR) cis-acting elements which are probably involved in posttranscriptional regulating of SV2Bb translation. In conclusion, alteration of SV2B(b) expression appears to be involved in processes of neuronal degeneration'' (see also [67]). We note that SV2B is only expressed in vesicles that undergo calcium-regulated exocytosis [68] and is a regulator of synaptotagmin 1 [69], which is a synaptic calcium sensor with a role in neurotransmitter release previously studied in AD [70,71,72,73,74,75]. We present a number of genes related to synaptic function and neuronal plasticity which are increasingly down/up regulated later on the manuscript and on the supplementary material (File S3 Sheet 'Synapse').
Analysis of the 1,372-probe signature reveals alterations in calcium and insulin signalling Using GATHER, we have identified 32 genes in the Calcium signalling pathway http://www.genome.jp/dbget-bin/show_pathway? hsa04020 (p-value,0.009). They are ADCY2, ADORA2B, AGTR1, ATP2A3, ATP2B1, ATP2B2, ATP2B4, AVPR1A, CALM1, CALM3, CREBBP, GNA14, GNAS, GRM5, HTR2A, ITPR1, ITPR2, LHCGR, NFATC1, PHKA2, PLCB1, PLCE1, PPP3CA, PPP3R1, PRKCB1, PTAFR, SLC25A6, SLC8A2, SYK, TBXA2R, TNNC2, and TTN. We cannot do enough justice in this manuscript to the several different hypotheses that point at imbalances/deregulation in calcium signalling and AD pathology. Instead, we contribute to these interesting discussions with our findings of genes related to this pathway within this group of 32 genes. The gene symbols in boldface can be mapped to the KEGG Pathway hsa04080, Neuroactive ligand-receptor interaction; those in italics can be mapped to KEGG Pathway hsa04310, Wnt Signalling. Being aware of the existing interest on Wnt Signalling and AD, we went back to the list of genes present in our (alpha,beta)-kfeature set signature and we identified others that can also be linked to Wnt signalling, like CSNK1G3, CSNK2A2, FRAT1[76,77,78,79, Figure 7. Heat map of the 50-probe signature and the transcription factors with best p-values, for the whole set of 50 probes and for the two groups considered. The samples and probes were sorted using the memetic algorithm given in [12], using the Euclidean distance. The transcription factors were obtained using Chang and Nevins' GATHER system to interpret genomic signatures [634]. The coloured cell and the number 1 indicate that the transcription factor has a binding motif with the gene for that row. The levels of severity as defined by Blalock et al. [635] are indicated in the first line: (0) Control, (1) Incipient AD, (2) Moderate AD and (3)  In addition, most of the remaining 32 genes in the Calcium signalling pathway can be mapped to KEGG Pathway hsa04070, Phosphatidylinositol signalling system (CALM1, CALM3, ITPR1,  ITPR2, PLCB1, PLCE1, PRKCB1), and Gap Junction (ADCY2,  GNA14, GNAS, GRM5, HTR2A, ITPR1, ITPR2, PLCB1,  PRKCB1).
This fact suggested that we should check how many genes were mapped to these pathways. We found that Phosphatidylinositol signalling system was indeed the third pathway with most ''hits'' in our signature, and also with other 12 genes (CDIPT, CSNK1G3 PIK3C3, PIK3R1, PIK3R4, PI4KB, PIP5K1A, PIP5K1C, PIP4K2C, PTEN, SKIP and TTK) which brings the total number to 19. We have also found (CCND3, CSNK1A1, CSNK2A2, CTBP1, CTBP2, FRAT1, FZD5, PPARD, PPP2CA, PPP2R2B, RBX1, SMAD3, TBL1X, TCF7L1, TCF7L2, VANGL1) bringing the total to 22 genes. We refer the reader to the supplementary material (File S3 Sheet 'Phosphatidylinositol signalling') for inspection of the individual pattern of expression of all these genes.
Together  173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199] in AD pathogenesis. Figures 10,11,12,13, and 14 illustrate down(up)-regulation of genes in these signalling pathways (Calcium signalling, Neuroactive ligand receptor pathway, WNT, Phosphatidylinositol and Insulin signalling, respectively). Figure 15 shows the expression of probes corresponding to genes for which there are known associations to synaptic function and neuronal plasticity. We refer the reader to the supplementary material (File S3) for more searchable information.
Transcription factors analysis of 1,372-probe signature reveals significant associations with the EGR/KROX family of proteins, MAZ, and E2F1 The analysis of the 1,372-probe signature indicates that they can be linked to putative transcription factors that have been previously implicated in AD and other neurodegenerative diseases. Using GATHER, we have observed that there is a strong association with motif V$KROX_Q6 (p-value,0.0004) with 719 out of 1294 genes in our signature; V$MAZ_Q6 (p-value,0.001, with 1003 genes); and V$E2F1_Q6_01 and V$E2F1_Q3_01 (with p-values which are smaller than 0.002 and 0.009 respectively). Of the 1294 genes associated with the 1,372 probes (by GATHER), more than half of them (656) have a motif for V$E2F1_Q6_01 and 603 have a motif for V$E2F1_Q3_01.
The involvement of the EGR/KROX (immediate early genes) family of proteins in the pathogenesis of Alzheimer's disease was first suggested in [221]. Studies of the behavioural consequences of stress have shown a link between the activation of the glucocorticoid receptor mediated response and EGR1, one of the members of this family [222]. It has been recently proposed that different members of the EGR/KROX family have dif-ferent roles in learning and memory and cognitive functions [223,224,225,226,227,228]. Mutant mice experiments showed that EGR1/KROX24 is required for the consolidation of longterm memory, while it is EGR3 the one linked to short-term memory [229], with EGR2 having perhaps other type of phenotypic characteristics not yet mapped [230]. In rat hippocampus, EGR1 decreases with aging [231]. In a recent study, it has been shown that initial playbacks of novel songs transiently increase EGR1 but that the observed response selectively Figure 9. 'Common-regulators' 50-probes' signature. The figure was obtained using Pathway Studio [569]. The program received as input the 50-probes displayed in Fig. 7 and automatically searched all the known putative common regulators relationships. The highlighted proteins are the 5protein signature (IL1-a, TNF-a, IL-3, EGF and GCSF) of [1]. We have also highlighted IL-6 (discussed in [1] in the context of results of classifiers that also use it) and CSF1, Colony-stimulating factor 1, (macrophage). doi:10.1371/journal.pone.0010153.g009    habituates after repetition of the stimulus, with a different expression profile after one day [232] (see [233] and also [234] in which the homolog of NEFM, one of our biomarkers of reduced expression with increasing 'AD severity' called NF-M, is showed to be involved in the development and/or maturation of the oscine song control system).
We found the following connection between EGR/KROX, E2F1 and MAZ transcription factors that makes their concurrent finding notable. A recent study of microRNA signature of prioninduced neurodegeneration [64] has shown that EGR1, E2F1 and MAZ might be also implicated in the putative deregulation of immune response related genes by miRNAs via modulation of transcriptional regulators in scrapie-infected mice. We leave these findings for the next section of the manuscript where we will discuss them and present a list of common differentially expressed genes in these two neurodegenerative processes.
The 1,372-probe signature contains a significant number of genes differentially expressed that are linked to synaptic function and neuronal plasticity The existence of several genes among the most correlated ones (NRXN1, SV2B, NEFM, etc.,) motivated us to try to identify which genes were present in the 1,372-probe signature that are also related to synaptic function and neuronal plasticity. We have identified 42 probes that can be divided into two groups, those that seem to be increasingly downregulated with AD severity (CABP1 [235,236,237,238,239,240,241,242,243], CADPS2 [244,245,246,247,248,249], COLQ [250], DMD [251,252,253,254,255,256] [322,323,324,325,326,327,328,329,330,331,332,333,334,335,336,337,338,339,340], SV2B [68,69,341,342,343,344,345,346,347,348,349,350,351,352,353,354,355,356,357,358,359]) and those that present an upregulation pattern (CASK [360,361,362], CDK5R1 [363,364,365,366,367,368,369,370,371,372,373,374,375,376,377,378,379], CHRNA1, CHRNA9, CHRNB3, CTBP2, DLG1/SAP97 [380,381,382,383,384,385,386,387,388], DLGAP2, GABRA5 [389,390,391,392,393,394], GABRQ [395], GLRA3 [396,397,398], GRIK3/GLUR7 [399], HOMER3 [400], ICA1 [401], ITGB1 [402,403], MCTP1 [404,405], PPP1CC [406], SNPH [407,408,409,410,411,412,413,414], SSPN [415], SYNC1, and USH1C [416,417,418]). The reader can consult the supplementary material (File S2) for the individual expression patterns of these genes. If, in agreement with Klemmer et al. [362], consider synapses as the most complex cellular organelle, with approximately 1500 proteins interacting in an activity dependent manner, we can argue that we must be inclusive with our list of references to help other researchers map the literature of their functions. Our aim is that experts can use this information to find ways of building novel testable hypotheses of AD neuronal plasticity impairment in the hippocampus. Our approach here has been to map what is currently known, and link it with the current biomedical literature, to facilitate experts that understand processes in detail.
We have already discussed some of the increasingly downregulated genes, another important candidate for further study is NRG1 (Neuregulin 1), a gene that has already been linked to several neuronal diseases. It is a candidate for susceptibility to schizophrenia and bipolar disorder (see [419,420,421,422,423,424,425,426,427,428,429,430,431,432,433,434] and references therein). There have been reported links of NRG1 with AD. BACE1 (beta-Site APP-cleaving enzyme) is necessary for the cleavage of the amyloid-beta precursor protein, and BACE1 participates in the proteolytic processing of NRG1 [435,436], and there exists some concerns about BACE1 inhibition as a potential therapeutic intervention due to its interaction with NRG1 and potential effects on remyelination [437]. In particular, NRG1 has been reported as a possible biomarker in cerebral spinal fluid, since its levels have been reported to be significantly increased in AD. Pankonin et al. suggest that: ''While (NRG1) is not detected in human serum, a novel neuregulin antagonist activity was identified in human serum that could have prevented its detection. These results suggest that human neuregulin is selectively targeted from cortical neurons to white matter extracellular matrix where it exists in steady-state equilibrium with cerebral spinal fluid where it has the potential to serve as a biological marker in human neuronal disorders'' [438]. NRG1 seems to collaborate with the ERBB4 receptor, and Li et al. propose that together they control glutamatergic synapse maturation and plasticity [439]. A single nucleotide polymorphism in NRG1 has also been associated as a risk factor to positive symptoms of psychosis in a proportion of late-onset AD [440]. With this evidence it is clear that NGR1 [439,441,442,443,444,445,446,447] as well as the whole panel presented here are excellent candidates for further studies due to their well supported role in synaptic function in health and disease states.
Other probes which present an upregulation trend that we would like to highlight are BCL2 [493,494], FYCO1 [495,496], PAX6 [111,497,498,499] (Figure 17), and QKI [500] (Figure 18). The increase of expression of these probes, together with SOX2, is intriguing as they are related to differentiation from stem cells and are considered critical in neurogenesis [501,502,503,504,505,506,507,508,509,510]. Our results support the combined use of them in tracking AD progression in this tissue. In addition, we have previously mentioned the relevance of EGR1 in coordinating a large number of genes that seem to be differentially expressed in this study. EGR1 also appears with a marked upregulation in severe AD patients (we refer to the supplementary material File S2 Sheet '1372 norm. +heat map+GO' for its gene expression profile). We found that this link is very important, as the homologues of EGR1, zif268, Egr-1 or ZENK, together with other members of the EGR family, are consolidating a key role in the neuronal plasticity in the brain [226,230,511,512,513,514,515,516,517,518,519,520,521,522,523,524,525,526,527,528,529,530,531,532,533,534,535,536,537,538,539,540,541,542,543,544,545,546,547,548,549] and links with AD and cognitive decline progression are starting to be reported [514,515,550,551,552,553,554].
At the same time, prospective studies should encompass some other genes which appear downregulated with increasing AD severity. Top of the list is perhaps LDB2/CLIM1 (LIM domain binding 2), recently pointed as a marker (with LMO4 [555,556]) of the control program of the development of neuronal subtype diversity of the cerebral cortex [557]. TRIM36 is another interesting candidate for further studies [558]. A gene that shares the same trend of dowregulation is CAMK1G (calcium/ calmodulin-dependent protein kinase IG) [559,560,561,562,563,564]. When analysing prefrontal cortical tissue from mice with inducible deletions of BDNF (Brain-derived neurotrofic factor), Glorioso et al. employed microarray gene expression profiling to show that there were alterations to early-immediate genes (including EGR1) and CAMK1G [563]. They conclude their manuscript stating that: ''while altered BDNF expression may not represent the primary disturbance in AD, changed expression of, or altered responsiveness to BDNF (and subsequently decreased SST levels) may represent a critical feature of Alzheimer's disease progression.'' VSNL1 (Visinin-like protein 1) [565], a CA++ sensor protein is also down-regulated (see Figure 19), a finding which is paralleled in the work of Youn et al. [566], who found similar changes in hippocampus. . By examination of the promoter regions of putative microRNA targets, they found that some transcription factor motifs were significantly enriched, E2F-1 (p-value = 6.01610 214 ), KROX (p-value = 9.34610 214 ), MAZ (p-value = 2.23610 211 ) and PAX6 (p-value = 1.76610 29 ). Our identification of EGR1/ KROX-24 and PAX-6 as upregulated with AD progression, and the identification of motif V$KROX_Q6, V$MAZ_Q6, V$E2F1_Q6_01, V$E2F1_Q3_01 as enriched in our signature were two contributing factors that motivated us to explore any further similarities that we could find.
In [64], an analysis of the predicted target genes of their microRNA signature, linked with differentially expressed genes in scrapie-infected mice [65] as well as two other publications [567,568], led Saba et al.
[64] to identify a network of de-regulated immune response-related genes. Additionally, they identified the putative transcription regulator genes that are targets of miRNAs similarly de-regulated. In essence, a possible hierarchy of deregulations of microRNAs, which, deregulated transcription factors that then, modify 1282 target genes. A Gene Ontology analysis also indicated that the ''data sets were found to be in the significant enrichment for genes involved in cell death, regulation of the cell cycle, nervous system development and function and cell signalling pathways.'' As a consequence, we have investigated if some of the 1,282 putative target genes of the miRNA signature of prion induced neurodegeneration also appear in our lists. Of those 1,282 genes we immediately noticed that there were 9 genes listed in our list of the 50 most correlated genes (Table 3). These genes are BCL11A, CSF1, DLG5, FOXO1, HBEGF, NRXN1, SERTAD2, SNRK and ZBTB20. Two of these genes, CSF1 (colony stimulating factor 1 (macrophage)) and HBEGF (heparin-binding EGF-like growth factor) appear to be conspicuous mediators of cytokine and growth factor signalling as Figure 9 illustrates (we obtained this network using Pathway Studio [569] as described in the previous section), and CSF1 and HBEGF seems to be increasing with AD severity. In opposition, the probe corresponding to NRXN1 (Neurexin 1, 209915_s_at) has decreasing expression ( Figure 20). Although no connection has been found between NRXN1 and AD yet, this gene has been implicated in autism [570,571,572  alcoholism dependence [582,583,584], and mental retardation [585]. SERTAD2 (SERTA domain containing 2), mentioned in the previous section, is also known as Transcriptional regulator interacting with the PHD-bromodomain 2, TRIP-Br2, a member of the TRIP-Br family of transcriptional regulators, required for the transduction of mitogenic signals and the execution of seruminducible E2F-mediated cell cycle progression [473]. In our data, the probe for SERTAD2 is increasing with AD severity. It has also been reported that overexpression of SERTAD2 is sufficient to transform murine fibroblasts and promotes tumorigenesis in athymic nude mice due to the deregulation of the E2F/DPtranscriptional pathway thanks to the upregulation of the key E2Fresponsive genes [474]. FOXO1 (Forkhead box O1) also appears upregulated with increasing AD severity, and has been reported as a negative regulator of EGR1 expression via the activation of the PI3K/Akt/Forkhead pathway [586]. The expression of FOXO1 is also induced by E2F1 [587]. The product of this gene has also been reported as a survival factor in deprivation-induced neuronal cell death [588,589] (see also the review in [590]). Although FOXO1 has not been previously implicated in AD, an exception may exist. van Der Heide et al. describe in [591] how the Forkhead transcription factors are involved in insulin signalling. The ''PI3K route'' is a name given to common signal transduction cascade that links neuronal survival, synaptic plasticity (and, as a consequence, learning and memory) [592]. This ''PI3K-Akt-FOXO1 mechanism'' and its role in neurons warrant the current intensive investigation [593,594,595,596,597,598,599,600]. From this group of 9 genes, seven of them (NRX1, SERTAD2, SNRK, HBGEF, FOXO1, CSF1, BCL11A) and QKI have been predicted to be targeted by mmu-mir128 by two or more microRNA prediction tools. We found this to be a connection that is worth exploring. Lukiw and Pogue have reported that following metal-induced reactive oxygen species production (by iron and aluminium-sulfate at nanomolar concentrations) upregulates miR-128 in human neural cells in primary culture [601]. They also report that, together with miR-9, mir-125a, mir-128 is upregulated in AD brain. In the previously cited reference Lukiw reported that: ''miR-9, miR-124a, miR-125b, miR-128, miR-132 and miR-219 are abundantly represented in fetal hippocampus, are differentially regulated in aged brain, and an alteration in specific micro-RNA complexity occurs in Alzheimer hippocampus.'' The expression of probes corresponding to PP2A and PP2B catalytic subunits (i.e. PPP2CA, Protein phosphatase 2 (formerly 2A), catalytic subunit, alpha isoform, and PPP3CA, Protein phosphatase 3 (formerly 2B), catalytic subunit, alpha isoform, Calcineurin A1) shows increasing downregualtion with the progression of AD., see Figure 21. This finding supports a role for downregulation of PPP2CA, PPP3CA in AD pathology .
Finally, in addition to the presence of hyperphosphorylated tau, the accumulation of Amyloid-beta (Abeta) peptide in brain tissue is a hallmark of AD [602]. The identification of the genes involved in the proteolytic processing of APP (beta-amyloid precursor protein), which in turn produces Abeta, is a subject of intense research. Researchers are currently looking at the alterations of APP cellular localization and endocytic trafficking as one mechanism that can modify the processing of APP to Abeta. LRPs are known to regulate APP's endocytic trafficking [603,604,605,606], and seem to be a hub of a number of mounting evidences on processes that link to cholesterol metabolism and atherosclerosis [607]. In our selected panel of 50 proteins we have one member of this family, LRP10 (low density lipoprotein receptor-related protein 10), as one of the most correlated gene expression profiles. In our list of 1372 gene probe signature we also have another member of this family, LRP1B (low density lipoprotein-related protein 1B (deleted in tumors)) [608], While LRP10 appears to be positively upregulated with cognitive decline an inverse relationship is observed for LRP1B.
LRPs are also known to linked to APP via a mechanism that involves the alternative splicing of APBB3/Fe65L2 [609,610,611]. Tanahashi and Tabira have proposed that the splicing of APBB3/ Fe65L2 alters the ability to bind with APP and low-densitylipoprotein-receptor-related protein. They propose that the secretion of beta-amyloid peptide Abeta40 and Abeta42 is increased following the overexpression of APBB3, but there are no visible changes of half-life and maturation of APP, or the secretion of secreted APP [612]. In our dataset, we observe APBB3 expression being upregulated with the increasing cognitive decline, following the same pattern of LRP10.
Polymorphisms on these genes have previously been linked to AD. Tanahashi, Asada and Tabira have reported an association between a polymorphism in APBB3/Fe65L2 and early-onset AD [612] (the link between APBB3 and AD is being increasingly explored, we refer to [613,614,615,616] for further references). Using 500K SNP microarray technology, Poduslo, Huang and Spiro have identified haplotypes in LRP1B as significant for successful aging without cognitive decline in a study involving individuals that were 85 years old or older, had MMSE scores greater than 26, no history of dementia in their families, and no major illnesses (i.e. no cardiovascular problems, diabetes, obesity, or major cancer diseases) and most of them had normal cholesterol levels. Their genome-wide association screening compared these individuals with those that have late-onset AD [617]. Poduslo et al. have suggested that if the decreased production of Abeta42 in successful aging is due to the haplotypes they describe, then LRP1B may be a new target for treatment of AD [608,617], Taken together these results indicate that integrative bioinformatics analytic Figure 19. The expression of a probe for VSNL1 (Visinin-like protein-1) shows increasing downregualtion with AD severity. VSNL1, a neuronal calcium sensor that has received recent attention in AD [636,637,638,639] has also been linked to model systems of schizophrenia, where it has been found upregulated in hippocampus [640]. A previous result by Schnurra et al. raised the possibility that the redution of VSNL1 expressing neurons indicate a selective vulnerabilty of these cells, since they observed that VSNL1 expression enhanced hyperphosphorylation of tau protein (in contrast with nontransfected or calbindin-D28K-transfected cells) [641]. In 2001, Braunewell et al. had already reported the reduction of VSNL1immunoactive neurons in the temporal cortex of AD patients as compared with controls [642]. doi:10.1371/journal.pone.0010153.g019 approaches will be needed to elicit the interactome of LRPs and their role in AD.

Conclusions
This re-analysis of the microarray dataset hippocampal gene expression contributed by Blalock et al. has shown that there exist a relatively large number of probes (1,372) that present a clear pattern of either up or down regulation with increasing AD severity. The signature reveals alterations in calcium, insulin, phosphatidylinositol and Wnt-signalling. Among the group of most correlated gene probes with AD severity we found some linked to synaptic function, neurofilament bundle assembly, neuronal plasticity and inflammation.
A transcription factors analysis of 1,372-probe gene expression signature reveals significant associations with the EGR/KROX family of proteins, MAZ, and E2F1. The gene homologous of EGR1, zif268, Egr-1 or ZENK, together with other members of the EGR family, are consolidating as key players in short and longterm memory and neuronal plasticity in the brain. We have also uncovered a large consensus of this gene expression signature with current genes putatively involved in AD progression. Our results also indicate a degree of commonality between putative genes involved in AD and prion-induced neurodegenerative processes that warrants further investigation.

Dataset
In this contribution, we have used a MIAME compliant, Affymetrix gene expression dataset that is public available and was contributed by Blalock et al [3] in 2004. We thank the authors of that publication for making this useful dataset available to the research community at large allowing further exploration and reanalysis.
The dataset is available from GEO Dataset Browser, accession number GDS1297 (http://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE1297). The Affymetrix human GeneChip, HG-U133A, containing 22,283 targets was used. The dataset is deidentified and the methods for disease classification, based on MMSE and NFT scores, are described in full detail by Blalock et al. in Ref. [3].  . It is possible to observe that one of the probes for NRXN1 (Neurexin 1, 209915_s_at) has decreasing expression with increasing AD severity. We have found no previous evidence of a connection of NRXN1 and AD, but this gene has been previously implicated in autism [570,571,572,573,574,575,576], schizophrenia [577,578,579,580,581], nicotine and alcoholism dependence [582,583,584], and mental retardation [585]. doi:10.1371/journal.pone.0010153.g020 and mental condition to determine if their were eligible for the study. When a mutual agreement existed, the individuals were visted in their homes to review and sign the informed-consent document (which was approved by the University of Kentucky Institutional Review Board). Participants also signed a donor card, and the visit also aimed to establish their baseline mental-status testing. Elegibility for the purpose of the study included having a Mini-Mental State Exam score above 24 [619], passing a series of cognitive tests, and a previous history of absence of neurological disease [620], as well as neither substance abuse nor major psychiatric illnesses. All eligible volunteers were 60 years of age or older and satisfactorily performed normal activities of daily living. The Wechsler Adult Intelligence Scale (Vocabulary) was also applied to exclude significant other medical diseases that could affect cognition and elegible participants must had no previous history of head injury with loss of consciousness.
The research participants that were deemed eligible also signed a form (in addition to the consent document) indicating their agreement to donate their brain to the Sanders-Brown Center on Aging. A full description of the methods used can be found in Brain Donation in Normal Aging Procedures, Motivations, and Donor Characteristics from the Biologically Resilient Adults in Neurological Studies (BRAiNS) Project [621].
Blalock et al. [3] categorized the samples in four groups, with a labelling that indicates different ''levels of severity''. These labels were decided based on the MiniMental State Examination (MMSE) and the Neurofibrillary Tangle count (NFT) of each sample [622]. Samples are then separated in the types 'Control', 'Incipient AD', 'Moderate AD' and 'Severe AD'.

Methodology
Our analysis method consisted of four steps: abundance quantization and filtering of probes; a feature selection algorithm to refine the probe selection; a Jensen-Shannon divergence computation; and finally, a correlation analysis. Each of these steps is described below.
As mentioned in the Results section, we only used the samples labelled as ''Control'' or ''Severe AD'' for feature selection, thus we have a two-class probe/gene selection task. We did not use the samples labelled as ''Incipient AD'' or ''Moderate AD'' for the probe selection steps. Those samples were only used in the final step, at the time of computing the correlation of the gene profile, across all samples, with the Jensen-Shannon divergences computed for the ''Control'' and ''Severe'' classes as explained later in this section.
For the first step, the quantization of the expression values, as well as for the initial data pruning, we used Fayyad and Irani's algorithm [626]. The heuristic algorithm minimises the feature-class entropy and discards genes according to the Minimum Description Length principle. The application of Fayyad and Irani's algorithm not only filters several thousand genes, it also provides thresholds for each probe remaining in the dataset. These quantized values of gene expression leave us with an instance of a combinatorial optimization problem, the (a, b)-k-Feature Set problem [13,627,628].
The (a, b)-k-Feature Set problem is a combinatorial optimisation problem introduced by Cotta, Sloper and Moscato [628] in 2004 to address the problem of feature selection in high-dimensional datasets. We solve an instance of this problem numerically using an integer programming formulation. This approach has been previously employed to obtain molecular biomarker signatures in Alzheimer's Disease [1,629], models of Parkinson disease [630], prostate cancer [631], electrode selection in EEGs [632], and elsewhere. To obtain mathematically proven optimal solutions of the integer programming formulation, the CPLEX commercial optimization solver was used. As in previous contributions of our group, we found gene expression signatures corresponding to values of a maximum and b maximal [1,13,627,628,633]. We refer the reader to these previous contributions for a detailed explanation of the methodology.
At this point, we have a selection of 1,372 probes, a set which we denote as V. For each sample m and probe i [ V, let f im be its expression value. We now define a probability distribution function (PDF) for each sample. For sample, m its PDF We can now compute an average PDF profile for samples in the ''Control'' and ''Severe AD'' groups, denoted by P C and P S respectively. Let C and S be the set of samples with the labels ''Control'' and ''Severe AD'' respectively. The average profile P c~f p p i (c) ,Vi [ Vg, is then: where N C represents the number of samples in class C. P P S is analogously defined.
The Jensen-Shannon divergence between two sample PDFs, i.e. samples l and k (P (l) and P (k) ) is defined as where Figure 21. It is well known that the square root of the JSD (sqrtJSD) is a metric, which means that for a given set of PDFs the following four properties are satisfied: i. sqrtJSD P (l) ,P (k) À Á §0, ii. sqrtJSD P (l) ,P (k) À Á sqrtJSD P (k) ,P (l) À Á , iii. sqrtJSD P (l) ,P (k) À Á 0uP (l)~P(k) , iv. sqrtJSD P (l) ,P (k) À Á zsqrtJSD P (k) ,P (m) À Á §sqrtJSD P (l) ,P (m) À Á : Once the sqrtJSD between each patient and the two average profiles ( P P C and P P S ) has been computed, the genes most correlated with these metrics can be uncovered. We used the Spearman rank correlation, which is a well-known non-parametric method, and can thus be used even when the data does not satisfy assumptions about normality, homoscedasticity and linearity.

Supplementary Material
Supplementary 'File S1' provides a glossary of each gene referenced in this paper including synoms and refrences to iHOP (http://www.ihop-net.org/).
The results referenced in this manuscript are provided in supplementary 'File S2' and 'File S3' in Microsoft Excel format.

Supporting Information
File S1 IHop Glossary of Genes.