Research ArticleNeuroscience

Genome-Wide Analysis of a Wnt1-Regulated Transcriptional Network Implicates Neurodegenerative Pathways

See allHide authors and affiliations

Science Signaling  04 Oct 2011:
Vol. 4, Issue 193, pp. ra65
DOI: 10.1126/scisignal.2002282

Abstract

Wnt proteins are critical to mammalian brain development and function. The canonical Wnt signaling pathway involves the stabilization and nuclear translocation of β-catenin; however, Wnt also signals through alternative, noncanonical pathways. To gain a systems-level, genome-wide view of Wnt signaling, we analyzed Wnt1-stimulated changes in gene expression by transcriptional microarray analysis in cultured human neural progenitor (hNP) cells at multiple time points over a 72-hour time course. We observed a widespread oscillatory-like pattern of changes in gene expression, involving components of both the canonical and the noncanonical Wnt signaling pathways. A higher-order, systems-level analysis that combined independent component analysis, waveform analysis, and mutual information–based network construction revealed effects on pathways related to cell death and neurodegenerative disease. Wnt effectors were tightly clustered with presenilin1 (PSEN1) and granulin (GRN), which cause dominantly inherited forms of Alzheimer’s disease and frontotemporal dementia (FTD), respectively. We further explored a potential link between Wnt1 and GRN and found that Wnt1 decreased GRN expression by hNPs. Conversely, GRN knockdown increased WNT1 expression, demonstrating that Wnt and GRN reciprocally regulate each other. Finally, we provided in vivo validation of the in vitro findings by analyzing gene expression data from individuals with FTD. These unbiased and genome-wide analyses provide evidence for a connection between Wnt signaling and the transcriptional regulation of neurodegenerative disease genes.

Introduction

Wnts constitute a large family of secreted, but spatially restricted, lipophilic signaling molecules. They exert wide-ranging biological activities by initiating equally varied signaling cascades (1, 2). The best-studied (“canonical”) Wnt signaling pathway involves the disheveled (DVL) and Axin-mediated inactivation of glycogen synthase kinase–3β (GSK-3β), producing an attendant increase in β-catenin, which activates lymphoid enhancer factor (LEF) and T cell factor (TCF) (TCF/LEF)–dependent transcription. Negative feedback, most prominently through increased expression of the intracellular canonical signaling antagonists AXIN2, naked-1 (NKD1), and naked-2 (NKD2), tightly controls the overall activity of this pathway (13). An emergent property of this complex regulatory machinery is that the activity of the Wnt signaling pathway oscillates in the face of sustained Wnt binding (47).

First identified as potent oncogenes, Wnts encode master regulators of fetal brain development and continue to govern neuronal growth and survival in the adult brain (8, 9). For example, Wnt1 is essential for neural crest induction and proper development of the midbrain dopaminergic system (911) and inhibits apoptosis in various cell types (1215). Furthermore, Wnt1, through the canonical signaling pathway, maintains the multipotency of human neural stem cells (16) and enhances neural progenitor proliferation and differentiation (17).

Wnt directly binds to more than two dozen partners, activating at least six different second messenger cascades in addition to the canonical pathway (1). As the number of signaling pathways known to be initiated by Wnt has increased, so too has the number of identified steps in each pathway and the number of components involved in mediating each step [for instance, β-catenin alone has more than three dozen known binding partners (18)]. Furthermore, the potential for crosstalk between Wnt and other signaling cascades has grown so extensive that these cascades resemble signaling webs more than directed pathways. The richness of Wnt signaling has also frustrated attempts to discern the most biologically important roles for this pathway from a myriad of plausible scenarios. Thus, it is surprising that few studies have systematically studied this pathway from a genomic standpoint.

Wnt signaling has also been implicated in several forms of neurodegenerative disease, most prominently Alzheimer’s disease (AD) and frontotemporal dementia (FTD), through various direct and indirect mechanisms (1923). For example, Wnt modifies the activity of presenilin and microtubule-associated protein tau (tau; MAPT), proteins encoded by two genes that when mutated cause either AD or FTD, respectively. Presenilin deficiency leads to β-catenin stabilization (24), and AD-causing presenilin mutations also disrupt β-catenin function (25), providing a direct link between canonical Wnt signaling and genetic factors that cause AD. Pathological tau phosphorylation is also a hallmark of both AD and a subset of FTD cases, and GSK-3β is a tau kinase (26, 27). GSK-3β–mediated tau phosphorylation accelerates tau-induced neurodegeneration in animal models of tauopathy (21), whereas GSK-3β inhibitors have been touted as potential therapeutic agents in tau-related dementias, such as FTD and AD (28, 29). Wnt signaling therefore may provide a bridge between neurodevelopment and neurodegeneration (30, 31).

To develop an unbiased understanding of the global effects of Wnt signaling in the nervous system from a systems biology perspective, we measured higher-order gene expression patterns in primary human neural progenitor (hNP) cells during the early phases of Wnt signaling. This analysis revealed a complex pattern of transcriptional changes occurring over a 72-hour period, which uncovered a strong connection between Wnt1-dependent signaling networks and expression of genes implicated in the etiology of both AD and FTD, as well as other neurodegenerative conditions. Haploinsufficiency of GRN (granulin), one of the genes down-regulated by the Wnt1 signal, causes a dominantly inherited FTD; however, virtually nothing is known about GRN regulation or how the progranulin (PGRN) protein it encodes affects neuronal survival. As an experimental validation of our network predictions at the level of a single gene, we investigated the potential mechanism of this interaction. We discovered a reinforcing feedback relationship between the control of Wnt and PGRN abundance, where decreased PGRN increased WNT1 expression, whereas Wnt1 repressed GRN expression and decreased PGRN abundance. This further implicates aberrant Wnt signaling in the etiology of FTD and supports investigation of Wnt pathway modulation in the treatment of neurodegenerative disease.

Results

A genome-wide time course of Wnt-induced changes in transcript abundance

The initial aim of this study was to obtain an unbiased view of the time-dependent changes in gene expression caused by Wnt signaling. Because different Wnts can potentially activate multiple distinct signaling cascades, we focused on the effects of Wnt1, which is a consistent activator of the canonical pathway (32, 33). Our experimental design eliminates many of the limitations encountered in previous studies of Wnt-mediated transcription, including using few time points (34), relying on transformed or immortalized cell lines that may poorly reflect normal human tissue (35), or surveying only a small fraction of the transcriptome (36). Specifically, we measured genome-wide mRNA abundance with Illumina Human RefSeq-8 BeadArrays at multiple time points (2, 4, 6, 8, 24, and 72 hours) to optimally capture different early epochs after the initial Wnt signal (36).

After Wnt1 application, the individual time course of these changes in gene expression—which could reflect changes in transcription, mRNA abundance, or both—showed various patterns, including transients, monotonic changes, and oscillations (Fig. 1A). Although the expression of only a minor fraction (0.3 to 20%) of genes was changed at any given time point (Fig. 1B; t test, P ≤ 0.05), the expression of nearly two-thirds of sampled genes changed at some point during the 72-hour time course examined. Overall, the largest changes in gene expression, whether up or down, were apparent at 2, 6, and 24 hours; fewer changes were observed at 4 and 72 hours.

Fig. 1

Wnt1-induced changes in the human neural progenitor transcriptome vary across time. (A) Representative genes showing common patterns of changes in gene expression over 24 hours, normalized to initial values. (B) Genes with significantly (P ≤ 0.05; t test) increased (magenta) or decreased (cyan) message abundance after Wnt1 application. (C and D) qPCR validation of microarray data. (C) Bars are range-normalized changes in expression (2 hours versus control) showing high correlation in directionality of expression, array (ILMN) versus qPCR (ρ = 0.48, n = 24 genes; one-tailed P ≤ 0.01). (D) Time series plots showing high correlation between qPCR- and microarray-quantified transcript. (E) Canonical Wnt reporter: Wnt1 applied to hNPs stably expressing a LEF/TCF-dsGFP reporter (n = 4 cultures) (139) causes a time-dependent increase in canonical activity that peaks between 4 and 6 hours. (F) Overlap in gene expression changes across time. (Upper left) Plot showing the fraction of genes whose message increased at 2 hours, but consistently decreased at subsequent time points. (Top, left to right) Venn diagrams showing the overlap of genes whose message increased (magenta) at both 2 and 8 hours or decreased (cyan) at 6 and 24 hours. (Bottom) Venn diagrams of genes whose message increased at 2 hours and then down-regulated at 6 hours (left), down at 6 hours then up at 8 hours (middle), or up at 8 hours and down again by 24 hours (right).

At 2 hours, increases in gene expression predominated, whereas at 6 and 24 hours, most differentially expressed genes showed decreased mRNA abundance compared to untreated cells. After the initial wave of increased transcript abundance at 2 hours, we found that quantities of most of these mRNAs had declined at 6 hours, then increased again at 8 hours, only to return to near-baseline activity by 24 to 72 hours. Next, using real-time quantitative polymerase chain reaction (qPCR), we examined whether we could confirm the initial burst of increased gene expression and the overall patterns of changes in gene expression recorded by the microarrays. First, we examined changes in the expression of 25 genes at t = 2 hours. The genes that compose this sample showed a wide range of changes in expression at t = 2 hours rather than changing in lock-step with the bulk of genes at this time (Fig. 1B). Eighteen of 24 genes (75%) identified as changing in gene expression by microarray showed a consistent pattern of changes when analyzed by qPCR (Fig. 1C). After range normalization, the microarray and reverse transcription–PCR (RT-PCR) data exhibited a significant Pearson correlation coefficient (ρ = 0.48, n = 24; one-tailed P ≤ 0.01). Next, we examined the changes in abundance for a subset of six mRNAs across all time points (Fig. 1D). From time point to time point, the direction of expression changes was the same 83% of the time, with a high collective correlation of 0.68 ± 0.14 (median ± SD) between qPCR and microarray results. Together, these data indicated a high degree of reproducibility of the microarray-based measurements of changes in mRNA abundance across time.

Because Wnt1 is the archetypical ligand for the canonical pathway, we tested the temporal relationship between activation of this pathway and the observed pattern of transcriptional changes. In hNPs stably expressing a TCF/LEF reporter (16, 37), signaling through the canonical pathway peaked 4 to 6 hours after Wnt1 application, returning to baseline by 24 hours (Fig. 1E). This monophasic rise and fall in activity of the canonical pathway coincides predictably with the initial changes in gene expression (Fig. 1B), but contrasts with those changes at later time points. These data suggest that canonical signaling may initiate a more complex cascade of changes in either transcription or mRNA stabilization, so that later events are mediated by both canonical and noncanonical signaling. Consistent with this hypothesis, we observed substantial changes in the expression of components of the noncanonical Wnt signaling pathways (fig. S1).

Successive waves of increases and decreases in transcription could arise from a set of coherently oscillating genes or represent a more complex summation of asynchronously changing sets of genes. The latter interpretation is supported by a closer inspection of the specific genes changing at each time point. For example, of the 2848 genes showing increased expression at 2 hours, less than half showed a decrease in mRNA abundance at later time points, with most of those showing decreases at 6 hours (Fig. 1F, upper left trace, and table S1). There is also a limited overlap of the genes participating in the successive bursts of increased expression. For example, only 153 genes showed increased expression at both 2 and 8 hours, representing 5 and 14% of the total number of genes changed at those time points, respectively. Together, these data suggest the presence of multiple superimposed transcriptional programs that sum to form apparent waves of altered gene expression.

Within this complex picture, several time-dependent functional patterns emerge, including the expected Wnt1-induced transcription of many genes implicated in the canonical Wnt signaling pathway, including WNT5B, FZD (frizzled) 2, 4, and 8, AXIN1, GSK3B, CTNNB1 (β-catenin), TCF7L2, NKD1, and AXIN2 (table S2). The products of many of these genes stimulate proliferation, repress differentiation, or are direct feedback regulators of the Wnt pathway itself (for example, Axin2). These apparent waves of Wnt1-mediated transcription are consistent with the oscillatory-like patterns of Wnt target gene expression that underlie operation of the embryonic vertebrate segmentation clock (38, 39).

Pathway analysis

We performed a Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (40) as a first step in unbiased characterization of the function of differentially expressed Wnt1-induced genes. This analysis confirmed early enrichment (t = 2 hours) for genes involved in Wnt signaling (Fig. 2A, Table 1A, and table S2). Two-thirds of the genes in the canonical Wnt pathway represented in the KEGG differentially regulated after the Wnt1 signal, compared to the control condition, validating the activation of Wnt signaling in these experiments.

Fig. 2

Wnt1 modulates genes implicated in Wnt signaling and AD. (A) KEGG pathway: graphical summary of the diversity of Wnt signaling–related genes (canonical and noncanonical), significantly enriched (FDR ≤ 5%) at 2 hours after Wnt1 application (n = 59 of 151 KEGG Wnt genes). Significantly increased (magenta) or decreased (cyan) mRNA abundance (t test; *P ≤ 0.05). (B) KEGG AD pathways highlighting genes (magenta) whose message was significantly increased by Wnt1 at t = 2 hours (t test; *P ≤ 0.05). (See fig. S9 for scalable version of this figure.)

Table 1

GO analysis of changes in Wnt1-mediated gene expression reveals time-specific enrichment. (A) Selected KEGG pathways (129) or (B) selected highest-level (level 5) GO biological processes showing significant categorical enrichment at each time point (FDR ≤5%). GO analysis was performed in DAVID (40) on all significantly changed genes (t test, P ≤ 0.05), irrespective of the direction of the change. ATP, adenosine triphosphate; ncRNA, noncoding RNA; MAPK, mitogen-activated protein kinase; GTPase, guanosine triphosphatase; VEGF, vascular endothelial growth factor.

View this table:

Our analysis also revealed enrichment for genes and pathways involved in any of the following three broad categories (Table 1, A and B, and table S3): (i) cell proliferation or death, (ii) energy metabolism, and (iii) regulation of the transcriptional or translational machinery (representative KEGG pathways at early and late time points; figs. S2 and S3). Because more than 90% of genes that show changes in gene expression at each time point do so in the same direction, GO analysis was performed on all differentially expressed genes irrespective of the direction of the change. Few differences were noted when these changes were analyzed with respect to the direction of the change in expression (table S3B).

KEGG pathway analysis also revealed a significant enrichment of neurodegenerative pathways (Table 1A) including the AD (Fig. 2B), Huntington’s disease (HD) (fig. S3), and Parkinson’s disease (PD) pathways, consistent with the hypothesized role of Wnt in neurodegenerative disease (21, 23, 4146). Moreover, we saw enrichment for genes associated with AD in the National Institutes of Health Genetic Association Database (table S4). Indeed, nearly one of three of the genes in the KEGG AD pathway were among the earliest group of genes showing increased Wnt1-mediated expression (t = 2 hours), a biologically meaningful and statistically significant result (see Materials and Methods). Because this analysis was based on an unbiased and genome-wide derived overlap, these data provide a first tier of systematic evidence for a causal connection between Wnt signaling and regulating expression of genes implicated in neurodegenerative disease.

Time series analysis

Our initial differential expression analysis revealed complex patterns of changes in gene expression. Analyzing this level of biological complexity across time is complicated by the nature of the time series data itself. For example, time series data are often mathematically nonstationary and show autocorrelation, creating nonlinearities or discontinuities in the data, which limit the use of many statistical techniques that assume fixed or Gaussian probability distributions (47). A further complication of analyzing time series data is that expression patterns in pairs of genes are often not monotonically correlated, rendering inappropriate commonly used nonparametric tests of association, such as the Spearman rank correlation (48, 49).

To reach a more complete biological understanding of Wnt1’s effects, we analyzed these data in three separate ways, using distinct statistical techniques that are relatively insensitive to the underlying data’s distribution: (i) dynamic time warping (DTW) (5052), (2) independent components analysis (ICA) (53, 54), and (3) mutual information–based network analysis (55).

A systems-level analysis of Wnt1 expression trajectories using DTW

The signal processing method of DTW provides a flexible and robust metric for measuring the difference between two waveforms (5052). DTW creates this distance measure by locally compressing or stretching (warping) one trace to best match the other and then summing the distances of individual aligned elements. In this way, DTW provides high power to detect differences in time series by accounting for the trajectory of changes in expression over time, unlike a purely statistical technique, such as analysis of covariance (ANCOVA). To better establish the applicability of DTW to microarray time series data, we analyzed how noise or phase shifting affected DTW-based gene ranking or GO enrichment, first with simulated data and then as applied to our Wnt1 data set. First, we found that DTW distance increased linearly with random noise and in proportion to the length of random traces being compared (fig. S4, A and B). Next, we assessed the effects of phase shifting of identical waveforms on the resulting DTW distance (DTWdist) and compared this to the Spearman correlation, a more common metric in microarray analysis. Because zero is the lower boundary of DTWdist, increased phase shifting resulted in monotonically increasing distances, unlike the correlation coefficient, which crossed the zero point (fig. S4, C and D). This relationship reduces the ambiguity in interpreting DTW-measured distances that are zero (or very small), an advantage of DTW over correlations measures (that is, a DTW distance of zero necessarily indicates that two traces are highly similar, not merely in phase). When comparing less idealized transients of gene expression, DTW distance again conforms better to the common sense notion of association than does either the Spearman or the Pearson correlation (fig. S4E).

We used DTW to identify those genes most strongly influenced by Wnt1 signaling by measuring the warping distance between the time course of changes in gene expression in cells treated with Wnt1 and that in untreated cells (longer DTW distances represent traces that are more Wnt1-specific). When these distances were ranked from most to least Wnt1-specific (high to low), they initially fell off exponentially (with an exponential decay constant, τ = 0.94, and correlation coefficient, R2 = 0.99; fig. S5A), where one decay constant is equivalent to 300 genes. Using this natural breakpoint, we limited further analysis to the top-ranked 300 genes (table S5A). Sensitivity analysis revealed that this ranking predictably degraded with the addition of increasing noise (fig. S5, B and C).

KEGG pathway analysis of the native data revealed an enrichment for cancer-related genes (fig. S6), consistent with Wnt1 being a proto-oncogene. Likewise, a GO analysis of these genes demonstrated enrichment for genes implicated in neurogenesis, angiogenesis, or axon and dendrite growth (table S5B), and genes associated with neuropsychiatric disease (disorders with both psychiatric and neurodegenerative features) including the dementia-related genes APP (amyloid precursor protein), MAPT, and GRN (table S5C).

A functional ICA analysis of the Wnt1 expression network

Our studies of differential gene expression and ontological enrichment rest on the following assumptions: (i) Perturbations (such as drug application) alter existing biological processes by changing patterns of gene expression across time. (ii) Each gene product may be involved in many different biological processes (for instance, β-catenin has different functions at cell surface and nucleus). (iii) Multiple genes may subserve a single biological process. (iv) Only a few genes are essential for a given biological process, although many other genes may contribute to it. On the basis of this biological model (fig. S7A), we can formulate an analytical model (fig. S7B) for the purposes of computation. By viewing the overall pattern of gene expression obtained from microarray analysis as resulting from the sum of many different biological processes, our analysis becomes a special case of blind source separation, that is to say, recovery of a set of mixed signals, with little or no information about source signals or the mixing process. ICA, which extracts non-Gaussian statistical structures from high-dimensional data (54, 5659), is among the most widely used techniques for blind source separation and has been applied to the analysis of microarray data (58).

ICA can decompose the observed expression data into the most likely original independent sources of variation [a visual demonstration of ability of parallel decomposition ICA (pICA) to recover separate signals from highly convoluted data is illustrated in fig. S7C]. When applied to gene expression data, these independent components (the original sources of variation defined by ICA) will ideally correspond to a set of biologically coherent processes (that is, convergent physiological effects). For example, a “neuronal maturation” component might be dominated by groups of co-regulated genes involved in coordinating processes like neurotransmitter packaging, neurite outgrowth, and synapse formation.

We first performed pICA on the expression data from Wnt1-treated cells, pooled across all seven time points and replicates, producing seven statistically independent components (ICA1-7). In each component, each gene is valued by its statistical contribution to that component, termed gene loading. In each component, genes with gene-loading values below a threshold of 3.0 were discarded from the component, creating modules (ICM1-7) that captured only the most relevant genes in that component. Examples of normalized time series data for the top 10 genes, by gene-loading value, in each of seven independent components are displayed in Fig. 3A. In creating these modules (that is, loading-thresholded components), we used a cutoff far more stringent than those previously reported so as to be conservative in our analysis (58).

Fig. 3

Unsupervised parallel independent component analysis (pICA) blindly separates gene expression patterns by biological function. (A) Normalized, mean expression time course for the top 10 genes in each ICA component, ranked by gene loading. (B) GO enrichment via DAVID was performed on each independent component module after thresholding at a gene-loading level of 3.0. Colors delineate individual modules. Listed are the top nonredundant level 5 biological processes (129), disease ontologies (underlined), or KEGG pathways (boldface), with associated P values (n.s., not significant; n = number of genes per module). (C) Overlapping MiMI interactome networks built using the top 20 odd-ranked genes versus the top 20 even-ranked genes from ICM2 [Top Odd versus Top Even (left) or Bottom Odd versus Top Even (right)] networks built from the bottom 20 genes versus Top Even genes. The Top Odd network recovered significantly more Top Even genes (turquoise circles; n = 24) than did the Bottom network (magenta circles; n = 8). (D) Genes in each ICM were probed against the Broad Institute’s Molecular Signatures Database. The most highly enriched data set is presented for each module, as well as a representative sampling of other significantly enriched data sets. (See fig. S10 for scalable version of this figure. See table S8 for a full listing of GSEA-identified data sets.)

The individual genes that contributed most strongly to each component were subjected to GO and KEGG analysis. Using ICA, we observed modular ontologic enrichment, consistent with biologically meaningful gene clustering (Fig. 3B and table S6). For example, KEGG annotations delineated a single Wnt signaling–enriched component and a single AD-enriched component (ICM2). However, we also observed some redundancy, with multiple components enriched for genes related to altered translation or ribosomal function. Sensitivity analysis to evaluate the robustness of this technique revealed that the ability of ICA to recover relevant biological processes degraded gradually in the face of increasingly noisy expression data (table S7). For example, ICA still recovers a Wnt module even when the noise present in the underlying data is increased by a full SD.

As a further biological validation of ICA-based clustering, we sought to integrate our newly delineated expression modules into established gene-protein interaction networks, using the public interaction database MiMI (http://mimi.ncibi.org/MimiWeb/main-page.jsp) as a generator (60, 61). Starting with the gene ranking in the Wnt-specific ICA module ICM2, we created three interaction networks seeded with the (i) top 20 odd-ranked genes (genes ranked 1st, 3rd, 5th…19th; top odd), (ii) top 20 even-ranked genes, or (iii) the bottom 20 odd-ranked genes. We compared the degree of overlap between the odd versus even (top) networks and the top versus bottom networks, with the expectation that top genes (even or odd) should extract networks more similar to each other than to the network formed from lower-ranked genes (bottom). Compared to the bottom network, the top odd gene network recovered three times as many genes from the even network, including the seed gene granulin (GRN; Fig. 3C). Using a different experimental modality, these data confirm functional associations between highly ranked genes in our ICMs.

As a complement to the interaction network building, we sought to place these ICA-derived gene modules into the context of known changes in expression. We probed each gene set against the 3429 gene expression profiles available through the Broad Institute Gene Set Enrichment Analysis (GSEA) Molecular Signatures Database (Fig. 3D and table S8). The resulting gene set enrichment analysis was highly complementary to the GO analysis reported above. Consistent with Wnt1’s role in oncogenesis, the genes in each module were frequently overrepresented in cancer-related data sets. Similarly, there is good agreement between GSEA signatures and the more specific GO categories identified above (Fig. 3B). For example, module 1 was enriched for GO-annotated angiogenesis-related genes, and the genes from this module were overrepresented among genes expressed by regenerating blood vessels. In contrast, genes present in the Wnt-specific module (ICM2) showed the greatest overlap with genes showing increased expression in the brains of Alzheimer’s patients (Fig. 3D; Fisher’s exact test, P ≤ 5 × 10−11), strengthening the possible connection between Wnt1 signaling and dementia (Fig. 2).

The two methods, DTW and ICA, address two fundamentally different biological questions and were used to provide convergent lines of evidence for identifying genes mediating the Wnt1 response. DTW identified genes specifically influenced by Wnt1 treatment by comparing the trajectories of gene expression for Wnt1-treated hNPs versus untreated hNPs. As presented, DTW did not explicitly provide information about which genes might act together to mediate similar biological processes. Our ICA decomposed multivariate expression data from Wnt1-treated hNPs into parallel statistical components, a linear representation of maximally non-Gaussian source signals that are as statistically independent as possible. In effect, this parsed genes into groups that potentially serve similar biological processes. One would expect that the genes common to both lists would be enriched for those with the strongest biological relevance, even if the magnitude of the overlap is modest. That these lists do overlap indicates that the cohort of common genes plays a fundamental role in mediating Wnt1 signaling. Specifically, we observed a significant overlap (17%; hypergeometric test, P ≤ 10−14) in the gene lists obtained by DTW and ICA (table S9). Of the 14 genes identified by DTW that have a genetic association with neuropsychiatric illness (table S5C), 4 genes were shared by the Wnt-specific ICA module (ICM2): APP, GRN, AQP4, and CTSD, the first two being associated with early-onset dementia.

Evidence from our GO, MiMI interactome, and GSEA-based analyses pointed to the biological importance of the ICM2 module, because it was enriched for Wnt- or dementia-related genes, prompting us to more fully investigate the interplay among its constituent genes. First, we performed a weighted gene coexpression analysis (WGCNA) of the genes that compose ICM2 (6265). WGCNA performs unsupervised clustering of coexpressed genes into modules on the basis of topological overlap (TOM), a measure of the degree to which a given pair of genes share common neighbors (64). When applied to ICM2, WGCNA clustered the constituent’s expression profiles into four distinct submodules, each being described by a characteristic module eigengene (Fig. 4, A and B). These eigengenes are composite gene expression vectors derived from the first principal component of the measured mRNA abundances that are used to concisely describe the overall variability of the data (66, 67). Moreover, GO analysis of each submodule revealed unique patterns of significant ontologic enrichment, despite the small number of genes present in each submodule (Fig. 4C).

Fig. 4

Combined topological overlap–based clustering and dynamic Bayesian network construction links Wnt1 signaling with changes in dementia-related genes. (A and B) WGCNA clustering of ICM2 genes: Wnt1-stimulated expression time courses for the genes that compose the ICM2 module were averaged and then subjected to TOM-WGCNA–based clustering (A). This produced four submodules (Mustard, Brown, Blue, and Turquoise). (B) Submodule eigengenes, where singular value decomposition was used to extract a characteristic first principal component eigengene for each submodule. The y axis is eigengene expression. (C) GO analysis reveals functional uniqueness of individual submodules. (D and E) Dynamic Bayesian network (DBN) depicting causal relationships, within each module. (D) Overview of the DBN-based causality network. Edge color codes the original submodule. Node color indicates those genes identified by DTW analysis (magenta; n = 23). Outlined diamonds denote those genes whose expression was increased in the brains of Alzheimer’s patients (n = 20) (70). DLK1 forms the primary hub in this network. [Note: SORT1, like DLK1, is a binding partner for PGRN (73).] (E) A more detailed view of the DLK1 hub and its associated genes, revealing a significant overlap (hypergeometric P ≤ 0.001) with DTW-identified genes and a strong enrichment for genes with increased message in AD brains. Transcriptional drivers (blue dots), targets (magenta arrowheads). (See fig. S11 for scalable version of this figure.)

Advanced Bayesian methods allow one to assess causality among genes with covarying expression patterns. To further explore these relationships in ICM2, the Wnt-dementia module, we constructed a dynamic Bayesian network (DBN) for each WGCNA-clustered submodule, using the G1DBN package (68, 69). We then coded the resulting network to reflect overlap with either Wnt1-modulated genes previously identified by DTW or those ICM2 genes with increased expression in the hippocampus of Alzheimer’s brains (70). As an internal validation of this approach, we observed that the genes identified as “hubs” in the network were primarily drivers of, rather than targets of, “spoke” genes (Fig. 4, D and E). The two most prominent hubs were the delta-like 1 homolog (DLK1), which encodes a Notch antagonist, and the stem cell marker promonin-1 (PROM1). As illustrated, DLK1 forms a nexus connecting the “Blue” and “Turquoise” submodules. Isolating the genes forming connections to DLK1 (Fig. 4D) produced a subset of 103 genes, with statistically significant enrichment for both DTW-identified genes and genes with increased message in the hippocampus of Alzheimer’s patients (70). Closer inspection revealed that this submodule contains several canonical Wnt effectors (including DVL3, APC2, and TSC2), as well as two genes implicated in causing dementia [APP and GRN, which binds both DLK1 and SORT1 (sortilin-1; Turquoise submodule)], strengthening the link between Wnt1 signaling and dementia (7173).

Coexpression network analysis

Our previous network analysis revealed a highly probable link between Wnt1 signaling and dementia but was limited to a small number of genes. To put these observed changes into a more complete systems biology context, we performed a mutual information–based network analysis (55, 74) as a means of systematically exploring the higher-order structure of Wnt1-induced changes in gene expression across a much larger fraction of the genome. Mutual information–based network analysis provides a robust measure of the statistical association of gene expression patterns between gene pairs, but makes no assumptions about the underlying probability distributions (55, 7478). The first step in our analysis was to dimensionally reduce the data (and thereby decrease the data set) by pooling the genes identified by DTW, ICA, or strong deviation from statistical normality (that is, a Gaussian distribution), a common feature of genes that change expression over time, a procedure that yielded 4398 genes. This reduced data set comprises the union of the most specifically Wnt1-regulated genes and those showing the greatest changes in expression over time, thereby providing the context for Wnt1-specific changes. Next, we constructed the Wnt-induced mutual information–based network (see Materials and Methods and fig. S8), and then pruned it on the basis of MRnet, the information-theoretic method of maximum relevance–minimum redundancy(55). This revealed several prominent Wnt-related hubs, two of which contained known dementia-related genes (Fig. 5, A to C). For example, five Wnt genes, Wnt4, Wnt5B, Wnt6, Wnt8A, and Wnt8B, out of seven Wnts represented in this set, colocalized to a nine-gene cluster centered on the hub gene COL25A1. COL25A1 is a brain-specific membrane-bound collagen that inhibits fibrillization of β-amyloid and associates with senile plaques in AD brains (7982). Another hub that revealed a relationship between Wnt signaling and regulation of dementia-related genes (Fig. 5C) consisted of a nine-gene cluster centered on APOA4, which is associated with AD (8385) and is a paralog of the dementia gene APOA1 (8588), and included two genes involved in Wnt/β-catenin signal transduction, DVL2 and CTNNBL1. This cluster also contains the AD-related genes, presenilin1 (PSEN1) and death receptor 6 [DR6; tumor necrosis factor receptor (TNFR) superfamily 21: TNFRSF21], the FTD gene GRN, and TROY (TNFRSF19), which has been implicated in failure of the adult nervous system to regenerate (89). This grouping is all the more pronounced in that PGRN reduces inflammation by binding to TNFRs and that PGRN degradation is blocked by APOA1 (90).

Fig. 5

Wnt1 induces a gene expression architecture that correlates many well-known dementia genes with Wnt-related signal transduction. (A to C) Microarray-based gene expression data were generated from Wnt1-treated or untreated hNPs, followed at seven time points over 72 hours and repeated in triplicate. Displayed is the subset of the MINA (mutual information–based network analysis)–based network that contains connections that exceeded a threshold of 1.1 bits. (A) Multiple Wnts cluster around COL25A1, an AD-related gene. (B) CTNNB1 (β-catenin) and CXCR4 hubs. (C) Neighborhood containing a dense cluster of dementia [magenta; presenilin (PSEN1), GRN, APOA4, DR6 (death receptor 6)] and Wnt transduction–related genes. Nodes are color-coded to reflect genes implicated in neural development (orange), dementia (magenta), Wnt signaling (cyan), or diseases distinct from dementia (yellow). [Note: CCNT1 (cyclin-T1) is a binding partner of PGRN.] (D) Loss of APOA4 dysregulates dementia hub genes: qPCR of hippocampal gene expression among wild-type and APOA4-null mice reveals significant changes in the expression of connected genes (blue), relative to glyceraldehyde-3-phosphate dehydrogenase (GAPDH), but not the other Wnt-dementia hub gene COL25A1 (orange). Values are fold changes in gene expression calculated with the ΔΔCt method (125) [*P ≤ 0.05; n = 4 brains, PCR-ANCOVA (125)]. (See fig. S12 for scalable version of this figure.)

Experimental validation of network predictions

To test whether the in vitro network architecture predicted in vivo relationships—specifically the clustering of dementia- and Wnt-related genes—we examined the effects of loss of the hub gene APOA4. Using qPCR, we compared the expression of genes predicted to be highly correlated with APOA4 in the hippocampi of APOA4 knockout to wild-type mice. Loss of APOA4 substantially altered the expression of many neighborhood genes (that is, genes one to two degrees separated from APOA4), including PSEN1, TNFRSF21 (death receptor 6), TNFRSF19 (TROY), DVL2, and DVL3, validating network predictions in vivo (Fig. 5D).

Thus, these unbiased expression analyses validated, at the whole-genome level, previous hypothesis-driven work, suggesting a connection between Wnt signaling and neurodegenerative disease (43, 45, 46, 9194). Furthermore, both analysis of DTW-differential expression and mutual information–based network analysis revealed a higher-order relationship between GRN and Wnt pathway genes in hNPs. Little is known about the biological function or regulation of GRN, except that haploinsufficiency causes a dominantly inherited form of FTD (9598). Therefore, we further investigated the connection between Wnt and GRN.

Experimental characterization and validation of the Wnt-GRN (PGRN) relationship

We knocked down PGRN in hNPs using lentivirally transduced RNA interference (RNAi) hairpins (Fig. 6, A to C), observing a 50% decrease in PGRN abundance, using gene-specific RNAi, and no off-target effects on PGRN, using control RNAi (Fig. 6B). Decreasing PGRN abundance significantly altered the abundance of the mRNAs encoding several different Wnts, including WNT1 (Fig. 6C), which was increased nearly fourfold by PGRN knockdown.

Fig. 6

PGRN knockdown increases Wnt1 expression. (A) Lentiviral vector using CAG promoter is not silenced in hNPs and knocks down PGRN: hNPs were infected with lentiviral vector expressing DsRed2 and RNAi against GRN. The color plate illustrates hNPs immunolabeled for PGRN (green) and DsRed2 (red), under control conditions. Scale bar, 25 μm. (B) Western blot analysis of PGRN: At 2 weeks, cell lysates were collected, and PGRN abundance was analyzed by Western blot, demonstrating PGRN knockdown (n = 4). (C) qPCR for Wnt1 mRNA showing that GRN knockdown significantly alters Wnt expression: Bars indicate fold change in Wnt1 expression when normalized to either the general internal control β-actin (reflecting total cells) or the vector-derived DsRed (reflecting RNAiGRN expression; means ± SEM). n = 3; t test, *P ≤ 0.05. (D and E) Canonical Wnt signaling dose-dependently inhibits PGRN transcription and expression: (D) Wnt1 significantly represses baseline activity of GRN promoter in HEK293 cells transfected with a GRN luciferase reporter [n = 4; analysis of variance (ANOVA), *P ≤ 0.05]. (E) Western blot analysis reveals linear dose-dependent reduction of PGRN abundance in hNPs 4 hours after bath application of Wnt1 or lithium (Li) [n = 4; doses indicated; Pearson R2 = 0.99 (Wnt) and 0.94 (lithium)].

To determine whether Wnt’s effects on PGRN depended on changes in expression, we assessed Wnt1 regulation of GRN promoter activity. We found that Wnt1 significantly repressed the activity of a human GRN luciferase reporter overexpressed in a human cell line [human embryonic kidney (HEK) 293]. The time course of changes in reporter activity was parallel to that of the Wnt1-dependent changes in GRN expression observed in primary hNPs. Wnt1 repression of GRN-luciferase activity was maximal at 6 hours (Fig. 6D), consistent with studies demonstrating that β-catenin abundance peaks 4 to 6 hours after Wnt application in HEK293 cells (36).

Next, we examined the relationship between PGRN and Wnt1 signaling in primary hNPs. We found that both Wnt1 and lithium, a GSK-3β antagonist that mimics canonical Wnt signaling, decreased PGRN abundance in hNPs in a dose-dependent manner (Fig. 6E). Along with the microarray and qPCR results above, these data indicate that PGRN and Wnt reciprocally regulate each other’s expression in hNPs.

Finally, we sought to validate the connection between Wnt signaling and GRN expression by means of expression data collected from GRN-haploinsufficient demented patients. This data set (99) is composed of whole-genome expression data from the frontal cortex, hippocampus, and cerebellum of patients with either sporadic FTD or GRN-haploinsufficient FTD, or controls of similar sex and age distribution. We observed that various Wnt-related ligands and receptors, and all four LEF/TCF transcription factors implicated in the canonical pathway, were differentially expressed in tissue from either the frontal cortex (FTX), the hippocampus (HIP), or both, of GRN-haploinsufficient FTD patients, compared to their expression in tissue from the FTX and HIP of sporadic FTD patients (Fig. 7A). Moreover, WNT1 expression was increased in tissue from the FTX and HIP of GRN-haploinsufficient patients compared to tissue from the cerebellum, a region not affected by FTD (FTX: log-fold ratio, 1.6; P = 0.006; HIP: log-fold ratio, 2.42; P = 0.053). In contrast, WNT1 expression was not significantly increased in the FTX or HIP of sporadic FTD patients.

Fig. 7

Canonical Wnt transcription factors are among genes most specifically correlated with frontocortical expression in GRN-haploinsufficient FTD. (A) Wnt signaling receptors and ligands differentially expressed in the frontal cortex of GRN-haploinsufficient FTD patients versus patients suffering from sporadic FTD (Bayesian t test, P ≤ 0.05). (B) Disease eigengene: This is a graphical representation of the frontal cortex GRN-specific vector used to correlate disease state and anatomical location with expression of individual genes in the full FTD data set. The data set comprises whole transcriptome gene expression in the frontal cortex (magenta), hippocampus (blue), or cerebellum (orange) of patients with either GRN-haploinsufficient FTD (dark) or sporadic FTD (medium) or controls (light) (99). (C) KEGG Wnt/β-catenin signaling–related genes among the total set of 392 genes significantly correlated either positively (orange; increased message) or negatively (blue; decreased message) with the synthetic FTX-GRN eigengene depicted in (B). Genes binding β-catenin (bold). Significance was based on Spearman rank correlations between each gene and the disease-specific synthetic eigengene (ρ ≥ 0.40, *P ≤ 0.001, n = 56 arrays).

Finally, we sought to identify the changes in gene expression that correlate most specifically with GRN-mediated FTD and determine whether Wnt-related genes were overrepresented in this set. To accomplish this, we first created a synthetic eigengene to represent those genes exclusively exhibiting differential expression in the frontal cortex of GRN-haploinsufficient FTD patients, but not in sporadic FTD patients, nor in the hippocampus or cerebellum or any patients (Fig. 7B). By comparing idealized pattern of gene (eigengene) expression against measured mRNA abundance of the other genes on the array, we found that 392 genes were highly correlated with GRN haploinsufficiency and frontocortical localization. DAVID-based GO analysis (40) of this gene set revealed significant enrichment for genes involved in Wnt receptor signaling (EASE P = 0.001). Ingenuity-based analysis yielded a similar enrichment for Wnt/β-catenin signaling genes (P = 0.004, Fisher’s exact test). Furthermore, among this gene set, we found several β-catenin binding transcription factors, including LEF1 and TCF7L2, two of the four best-described targets of canonical Wnt signaling (Fig. 7C). These data validate and extend the expression data, confirming a connection between Wnt signaling and neurodegenerative pathways, including those associated with PGRN deficiency in humans.

Discussion

This study was designed to expand understanding of the role of Wnt signaling in human neural development using hNPs as a model system. Because Wnts can activate various intracellular signaling pathways, we used an unbiased systems biology approach to identify the Wnt1 targets most relevant to neural stem cells. We found that Wnt1 altered transcript abundance of a large number of genes in an oscillatory fashion. Although Wnt1 affected the expression of only a minority of genes at any given time point, it altered the expression of more than half of the hNP transcriptome during the 24-hour period after treatment. We showed that Wnt1 broadly affected most noncanonical Wnt signaling cascades, despite being the archetypal ligand for initialing canonical (β-catenin) signaling. Moreover, Wnt1 activation of these canonical and noncanonical Wnt pathways appeared to widely influence expression of genes mediating cellular metabolism (for example, steroid, purine, and ribosome synthesis) in addition to their more well-known functions, proliferation or cell fate. Using fetal progenitors as a model system, we found that Wnt1-responsive genes showed a strong collective association with sets of genes implicated in neurodegenerative diseases associated with aging, most prominently AD and FTD.

We used three complementary statistical workflows—DTW, ICA-WGCNA-DBN, and MI-MRnet—to generate a consensus network view of the Wnt1 transcriptional program. Individually and collectively, DTW, ICA-WGCNA-DBN, and MI-MRnet suggested a connection between Wnt1 and PGRN expression and, by extension, a potential connection between Wnt signaling and GRN-mediated FTD. Using a humanized inducible model of GRN deficiency, we found that PGRN and Wnt1 reciprocally regulated each other’s expression. Moreover, we found that the expression of numerous Wnt signaling genes was increased in the frontal cortices of GRN-haploinsufficient patients compared to other brain regions, or the frontal cortices of individuals without GRN mutations, indicating that this connection is clinically relevant. The variety of altered Wnt signaling genes suggests that they mediate a complex network of biological processes in GRN-haploinsufficient neurons. Further in silico and experimental work will be required to delineate these processes and which Wnt signaling genes are most relevant to them.

The Wnt1 gene expression network

Transcriptional oscillations have become increasingly recognized as a key feature of synchronized developmental programs (6, 7). We observed that a single application of Wnt1 induced multiple waves of changed expression in cells derived from early fetal neural progenitors over a 72-hour period. This is akin to what has been observed with the transcription factor hairy and enhancer of split1 (Hes1), which drives variable period oscillations in transcription within neural progenitors of the mouse telencephalon (100). Similarly, Wnt3a produces oscillations within developing limb buds (57) that control the rhythmic production of somites (precursors of the vertebrae) (6, 101), thereby linking the segmentation clock and signaling gradients within presomitic mesoderm. We identified a set of oscillating genes implicated in a wide range of biological pathways, including pathways relating to angiogenesis, immune modulation, cell proliferation, energy metabolism, and neuronal maturation (36, 38, 39). Thus, we propose that Wnt might serve a dual role in the developing nervous system: Wnt may act within neuroblasts to initiate the set of parallel developmental programs that are required for the formation of a mature neuron. In addition, Wnt secreted by neuroblasts could synchronize neural maturation with the development of necessary supporting tissue, such as vasculature and glia (102).

Regulation of PGRN expression

PGRN, like Wnt, was first described for its role in oncogenesis, and later shown to have neurotrophin or growth factor–like activity (103106). This led to the hypothesis that GRN haploinsufficiency leads to a relative deficit of PGRN and thereby the death of neurons within specific brain regions, such as the frontal cortex. However, more recent evidence suggests that active suppression of neuronal GRN expression, not overall reduced amounts of PGRN, precipitates frontal lobe degeneration (99, 107). Our findings that Wnt1 and PGRN reciprocally regulate each other’s expression provide a biological model that could potentially explain this active, neuron-specific decrease in PGRN.

We found that suppression of GRN expression, essentially mimicking FTD-associated human GRN haploinsufficiency, increased WNT1 expression, whereas Wnt1 suppressed GRN expression in otherwise normal neuronal cells. Hypothetically, a haploinsufficiency-induced reduction in PGRN could set up a negative feedback loop that actively suppressed GRN expression (beyond the decrease in GRN resulting from its haploinsufficiency). Wnts can act as autocrine or paracrine signals because they exhibit limited diffusion away from the source (108, 109). Therefore, increased neuronal Wnt production would selectively suppress expression of GRN in neighboring neurons, rendering them incapable of compensating for haploinsufficiency (that is, unable to increase expression from their unaffected allele). This model could explain several clinical pathophysiological observations including (i) the concentration of PGRN in the blood of individuals haploinsufficient of GRN levels can be far below the 50% that would be expected based on the decrease in gene dosage (110112), and (ii) cells in different brain regions of affected patients vary in their ability to effectively maintain sufficient levels of PGRN synthesis (99, 107).

Whether increased Wnt signaling improves survival of GRN-haploinsufficient neurons remains an open question. However, it is possible that these cells engage Wnt as an alternative pro-survival signal. Compensating for the loss of PGRN-mediated tropism by means of Wnt could partially normalize function early in the disease process. However, it could ultimately contribute to cortical atrophy by further decreasing GRN expression. A similar model of “compensatory engagement” has been proposed to explain the indolent course and ultimate cognitive system failure in chronic neurodegenerative diseases like Alzheimer’s (113, 114). Because Wnt1 asynchronously activates functionally distinct transcriptional programs, it is possible that the pro-survival effects of Wnt1 could be pharmacologically separated from its repression of GRN expression. Therefore, for individuals with GRN haploinsufficiency, targeting the right subset of Wnt1 target genes very early in life could conceivably prevent the onset of disease, without affecting the protective effects of Wnt signaling.

In summary, our data support a role for Wnt signaling in clinical dementia in general, and FTD more specifically. Together, the combination of functional genomic, bioinformatic, and cell biological data supports a possible mechanism of GRN dysregulation in haploinsufficient patients. These data lay the preclinical foundation for developing Wnt modulators as potential treatments for PGRN-mediated FTD and other neurodegenerative conditions.

Materials and Methods

Primary cell culture

We generated human fetal neural progenitors as previously described (16, 115, 116) and propagated them in human neural stem cell proliferation medium consisting of Neurobasal A, 10% BIT 9500 (Stem Cell Technology), fibroblast growth factor 2 (FGF2) (20 ng/ml; PeproTech), epidermal growth factor (EGF) (20 g/ml; PeproTech), leukemia inhibitory factor (LIF) (2 ng/ml), and heparin (2 μg/ml). To eliminate the potential confounding effects of FGF2 and EGF (117), we replaced proliferation medium with minimal growth factor medium [Neurobasal A, 5% BIT 9500, and heparin (2 μg/ml)] 48 hours before Wnt1 application. After this growth factor washout, the medium was again replaced 100% to eliminate autocrine factors that might have accumulated (16). hNPs were divided into treatment and control groups, and recombinant Wnt1 (400 ng/ml; PeproTech) was applied to the treatment groups. hNPs were harvested at time = 0, and then both Wnt1-treated and untreated control NPs were harvested at 2, 4, 6, 8, 24, and 72 hours. Three biological replicates were collected at each time point after Wnt1 treatment. For PGRN knockdown experiments, hNPs expressing pLCR-RNAiPGRN were differentiated in medium supplemented with 1% fetal calf serum (FCS; Hyclone), 500 μM all trans-retinoic acid (Sigma), 10 μM forskolin, NT3 (20 ng/ml), and brain-derived neurotrophic factor (BDNF). Functional confirmation of PGRN-related effects was performed on hNP lines derived from separate individuals than those whose cells were used for the time series microarray experiments. HEK293 cells were propagated in Dulbecco’s modified Eagle’s medium (DMEM), 10% FCS, and 1× Antibiotic-Antimycotic mix (Invitrogen). Recombinant human Wnt1 (PeproTech) was used at 200 ng/ml.

Viral preparation, infection, and knockdown

Lentivirus for overexpression or inducible knockdown of gene expression was produced as previously described (16, 118, 119) (see Supplementary Methods for details). Hairpins directed against GRN were designed, cloned into either pLCR or pTRIPZ (Open Biosystems) using the PCR-shagging protocol (120), and then tested for their ability to achieve greater than 50% knockdown. PGRN knockdown in TRIPZ-infected cells was induced by addition of doxycycline (2 μg/ml). APOA4 knockout mice were generated by homologous recombination and maintained as previously described (121, 122).

Immunodetection

Immunoblotting and immunocytochemistry of whole-cell lysates or cultures were performed by standard methods, essentially as previously described (123) (see Supplementary Methods).

Reporter assays

GRN promoter analysis was performed in HEK293T cells transected (1:50) with Renilla luciferase plasmid (pRL-EF) and either constitutively active firefly luciferase or pCMV-Tag4a PGRN expression plasmid (Switchgear Genomics) and then processed according to the manufacturer’s instructions (Dual-Luciferase Assay; Promega). Wnt activity was assessed essentially as previously described (16) (see Supplementary Methods for details).

Microarray analysis

Total RNA was harvested from hNPs with the RNeasy kit (Qiagen) according to the manufacturer’s protocol. Spectrophotometry [NanoDrop; A260/280 (absorbance at 260/280 nm), −1.8] and the Agilent Bioanalyzer were used to determine RNA concentrations and RNA quality. Total RNA (200 ng) was amplified, labeled, and hybridized to HumanRef-8v2 and v3 Expression BeadChip (Illumina). Data were preprocessed through BeadStudio (Illumina) to produce raw output files. All further processing was conducted in the R statistical computing environment (http://www.r-project.org/). When indicated, background correction, variance-stabilizing transformation (VST), and robust spline normalization (RSN) were applied using the Bioconductor package LUMI (http://www.bioconductor.org). Probes that varied in sequence between platforms were eliminated, leaving 18,396 probes that were examined in all experiments described herein.

Time point differential expression analysis

Microarray-based measurements of mRNA expression were log2-transformed, and then at each time point, the values from untreated control hNPs were subtracted from the measured values for Wnt-treated hNPs. By applying this procedure on a gene-by-gene basis, we created a data set in which genes whose expression was unaffected by Wnt1 treatment would be expected to have a normalized mean expression of zero. Significance of differential expression at individual time points was calculated with a zero-centered, one-way t test. Differential expression ratios and P values on the in vivo Human FTD data set [Gene Expression Omnibus (GEO) ID: GDS3459] were calculated in R with the Limma (Linear Models for Microarray) differential expression analysis package (124). The uncorrected P values were computed from a Bayesian moderated t test (124). For qRT-PCR, complementary DNA (cDNA) was generated from total RNA using random hexamers and SuperScript III (Invitrogen). Real-time PCR was conducted with SYBR-Rox SuperMix (Bio-Rad) and an HT-7900 (Applied Biosystems) and statistically analyzed in R with the t test or ANCOVA methods adapted from available SAS code (125).

GO analysis

GO analysis was performed with DAVID Functional Annotation Bioinformatics Microarray Analysis software (40), except for the sensitivity analysis where topGO was used (126, 127). In DAVID-based analyses, the reported P values are derived from the Expression Analysis Systematic Explorer (EASE) score probability, a modified Fisher’s exact test that is more conservative than the standard Fisher’s exact test in examining P values of gene lists (128). Where indicated, significance of overrepresentation was adjusted for multiple comparisons to control the false discovery rate (FDR) by means of the Benjamini-Hochberg step-down procedure for calculating the FDR in the case of independent tests or the approximated FDR tools provided in DAVID (40). All level 5 GO Biological Processes (BPs), the highest-level terms in the BP GOgraph structure (129), or KEGG pathways that remained significantly overrepresented (FDR = 5%) were reported. Ingenuity-based analysis reports enrichment probabilities based on the Fisher’s exact test.

Dynamic time warping

Whole time course analysis used a DTW-based approach (50, 51, 130) to calculate the distance between Wnt1-treated and control time series gene expression data on a gene-by-gene basis, using all 18,396 probes. In brief, DTW locally compresses or stretches (warping) one waveform to best match a reference waveform. The more warping that is required, the greater the difference between the two waveforms. Using this scheme, greater distances between trajectories for a given gene indicate a larger effect of Wnt1 treatment. Calculation of warping paths used the full data window, unconstrained endpoints, and the Euclidean distance to measure path differences (130). DTW sensitivity analysis: (Random traces) DTW distance was computed between randomly generated traces containing between 10 and 40 time points, and the same time course with increasing amounts of added Gaussian noise (μ = 0, 0 ≤ σ ≤ 1). DTWsimilarity (sDTW) was computed as follows: sDTW = 1/(1 + DTWDist).

Mutual information–based network analysis

Wnt1-induced networks were inferred in four steps on a refined gene list. (i) Dimensional reduction: The refined gene list was created by pooling the following: ICA-identified genes, the top 300 DTW-identified genes, and those genes with significantly nonnormal distributions, as assessed by either the Jarque-Bera skewness test or the Anscombe-Glynn kurtosis test (P ≤ 0.05). (ii) Network building: We calculated the weighted gene coexpression network on the basis of the entropy normalized mutual information (see below) distance (131) between each pair of genes in our refined list. First, data were discretized; that is, continuous variables were converted into discrete ones by equal frequency binning. Second, corrected marginal entropy (HX, HY) and joint entropy (HXY) were calculated from raw discrete marginal and joint probabilities with a James-Stein–type shrinkage estimator, which is robust for small sample sizes (132, 133). Raw mutual information (MI) measurements between gene pairs were calculated from the marginal and joint entropies according to the equation MI(X,Y) = H(X) + H(Y) − H(X,Y), expressed in natural bits (nats), where 1 nat = log2(e) binary bits. These were transformed into a maximum entropy-normalized dissimilarity metric (entropy-normalized MI) (131). This metric was calculated according to D'=1I(X,Y)max(Hx,Hy). Next, D′ was scaled according to MI(X,Y) = log2(e)*D′(Max − Min). (iii) Network pruning: We eliminated low probability connections from this network with the “maximum relevance, minimum redundancy” (MRnet) criteria (133). (iv) Final refinement: networks were further pruned by eliminating edges that fell below 1 bit of information (fig. S8). The resulting sparse network approximated a power-law distribution of connections, typical of the scale-free topology of gene-expression and protein networks (134, 135).

Parallel independent component analysis

Eigenarrays corresponding to parallel statistically independent components were calculated using an R implementation of the fastICA algorithm (136). pICA was used to extract seven independent components from the expression data for the Wnt1-treated hNPs. Because the ICA algorithm requires searching the maxima of a target function, the final output may reflect local rather than absolute maxima.

Modeling suggests that the ordering of genes in individual components is relatively, but not absolutely, stable for those genes whose gene loading exceeds a threshold of 3. To overcome this limitation (56), we created consensus components by reseeding the fastICA 250 times and then median-averaging the gene loading for each gene in each component. Because the fastICA output is not ordered, but the ranking of gene loadings for a given component was >90% identical from run to run, we used the Damerau-Levenshtein distance to match component lists to create seven consensus independent components.

Other network analyses

Topological overlap–based WGCNA networks and dynamic Bayesian network were performed as described (6265) in R (http://www.r-project.org/).

Disease eigengene analysis

The FTD data (99) were downloaded from the gene expression omnibus Web site (GEO ID: GDS3459) and then log2-transformed before further analysis. The disease-specific eigengene was calculated to reflect perfect up-regulation in the specimens from the frontal cortex of GRN-haploinsufficient patients, but unchanged expression in other regions and other disease states. The putative disease eigengene is itself a vector of length equal to the number of samples in the data set. Each element is limited to the values of {−1,0,1}, where these values correspond to either up-regulation, no change, or down-regulation in a given sample. Then, this vector is pairwise-correlated with the expression of every gene using the Spearman rank correlation. [This measure is more robust than the previously used Pearson correlation (62, 63, 135), being less sensitive to outliers and other nonlinearities in the data.] A Spearman correlation (ρ ≥ 0.40) was chosen, which is equivalent to a raw *P ≤ 0.001, as calculated with R statistical software. The underlying algorithm used by the software calculates the P value as follows: The distribution of Spearman’s rank correlation coefficients (ρ) in the null case (zero correlation) can be approximated by several methods; however, for larger sample sizes (n ≥ 30), it is nearly the t distribution (137, 138). Using the relationship tρ*(n2)/(1ρ2) (137, 138), we calculated that a ρ = 0.4, n = 56 produces a Student’s t ≈ 3.21. In turn, this t value and n = 56 are roughly equivalent to a P = 0.001, using tables readily available in most introductory statistics textbooks.

Supplementary Materials

www.sciencesignaling.org/cgi/content/full/4/193/ra65/DC1

Methods

References

Fig. S1. Wnt1 differentially modulates Wnt pathway genes across time.

Fig. S2. Wnt1 induces significant early enrichment of genes involved in survival signaling, energy metabolism, and biosynthesis (anabolism).

Fig. S3. Wnt1 induces significant later enrichment of genes involved in blood vessel development, axonogenesis, and Huntington’s disease.

Fig. S4. Dynamic time warping (DTW) sensitivity analysis of simulated expression data.

Fig. S5. DTW sensitivity analysis of experimental data.

Fig. S6. DTW identifies Wnt1 effects on cancer-related genes.

Fig. S7. Demonstration of robust blind source separation using parallel independent component analysis (pICA).

Fig. S8. Overview of mutual information–based network inference and sensitivity to thresholding effects.

Fig. S9. Scalable version of Fig. 2: Wnt1 modulates genes implicated in Wnt signaling and Alzheimer’s disease (AD).

Fig. S10. Scalable version of Fig. 3: Unsupervised pICA blindly separates gene expression patterns by biological function.

Fig. S11. Scalable version of Fig. 4: Combined topological overlap–based clustering and dynamic Bayesian network construction links Wnt1 signaling with changes in dementia-related genes.

Fig. S12. Scalable version of Fig. 5: Wnt1 induces a gene expression architecture that correlates many well-known dementia genes with Wnt-related signal transduction.

Table S1. Overlap of gene changes across time.

Table S2. Wnt1-induced immediate changes in known Wnt-related genes.

Table S3. Wnt1-mediated enrichment of GO and KEGG pathways varies across epochs.

Table S4. Wnt-1 induces changes in genes implicated in AD by genetic association studies.

Table S5. GO analysis of whole time series DTW-identified genes.

Table S6. Extended ICA gene list and ontology analysis table.

Table S7. ICA-GO sensitivity analysis.

Table S8. Gene Set Enrichment Analysis of ICA module genes probed against the Broad Institute’s Molecular Signatures Database (MSigDB).

Table S9. Unique, independent statistical analyses identify similar Wnt1-mediated changes in gene expression.

References and Notes

  1. Acknowledgments: We thank G. Konopka, K. Winden, R. Versano, and J. Miller for their technical expertise and invaluable discussions and J. Davis-Turak, F. Gao, and G. Coppola for assistance with BeadStudio data preprocessing. Funding: This work was supported by the John Douglas French Alzheimer’s Foundation (E.M.W.), National Alliance for Research on Schizophrenia and Depression (E.M.W.), National Institute of Mental Health (NIMH) K08MH74362 (E.M.W.), NIMH R01MH060233 (D.H.G.), National Institute on Aging 5R01AG026938 (D.H.G.), and NIH/National Institute of Neurological Disorders and Stroke Neurobehavioral Genetics Training Grant 5T32NS048004-05 (E.R.). Author contributions: D.H.G. and E.M.W. conceived the study. E.M.W. wrote the analytic software, collected the data, and wrote the initial draft of the manuscript. E.M.W. and D.H.G. provided financing, edited, and had final approval over the manuscript. D.L., E.R., and G.E.O. were responsible for qPCR, tissue culture, in vitro knockdown, and GRN expression experiments. E.M. and H.R. were responsible for the APOA4 knockout mouse colony and harvesting of brain tissue. Competing interests: The authors declare that they have no competing interests. Data availability: Microarray expression data files can be obtained from the NIH Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/).
View Abstract

Stay Connected to Science Signaling

Navigate This Article