For this study, we have selected five microarray datasets referred here as BL (Blalock et al., 2003), B7 (Burger et al., 2007), R7 (Rowe et al., 2007), B8 (Burger et al., 2008) and K9 (Kadish et al., 2009). They consist of a total of 287 arrays and used two different Affymetrix chip types, RG_U34a and RAE230a. The data represented young rats that were learning unimpaired and aged rats majority of which were learning impaired with some learning unimpaired animals. The BL and K9 studies were similar in design where only the unimpaired young and impaired aged animals were considered for comparison. The B7, R7 and B8 studies were similar in design where both young and aged groups contained controls animals that were learning impaired, e.g., cage controls, stress controls and controls for visual impairment. These datasets were quality checked and normalized using RMA methods (Bolstad et al., 2003; Gautier et al., 2004), had outliers removed and batch effect adjusted using age and spatial learning impairment as covariates (Uddin and Singh, 2013). For each dataset, aged and young samples were separated and assessed further for the presence of array outliers (Supplementary Table S1 and Supplementary Figures S1–S6). Since the WGCNA network construction method is correlation based, before proceeding with network analysis it was made sure that the correlations between genes in each dataset were reasonable as suggested in the literature (Miller et al., 2010). This was done by calculating Pearson’s correlations between the expression levels of each pair of genes in the aged or young preprocessed datasets and by plotting the correlation values in histogram plots (Supplementary Figure S7). All data preparation steps including WGCNA, GO, and statistical analyses were performed in R using appropriate software packages.

Using the preprocessed transformed data (genes in columns and samples in rows), gene networks were constructed for aged and young using the WGCNA R package (Zhang and Horvath, 2005; Langfelder and Horvath, 2008) following a slightly modified protocol based on the approaches previously described (Oldham et al., 2006, 2008; Miller et al., 2008, 2010). The overall network analysis process for a single dataset is described below.

A weighted correlation between two genes represents connection strength between the genes in a network. For each dataset, a network adjacency or connection strength matrix (network data) was created by taking the signed correlations of the gene expression values between each pair of genes raised to a power of beta. Beta is the weight, a soft threshold, and was determined in advance in such a way so that the resulting network follows approximate scale free topology. The values in the diagonal (self-correlation) were converted to zero.

To save computational time, genes were filtered out from a network adjacency matrix based on their connectivity (i.e., only genes with above average median connectivity were kept for network analysis). The overall connectivity for each gene (denoted by k) is the sum of connection strengths (weighted correlation) between that gene and all other genes in the network. It is scaled to lie between 0 and 1 and represents how strongly a gene is connected to all other genes in the network.

Following filtering an adjacency matrix contained genes with reasonably high network connectivity. This adjacency matrix was used to determine a network topological overlap, construct a hierarchical clustering dendrogram of 1—topological overlap, determine network modules using a hybrid tree-cutting algorithm, and to visualize network modules. In a co-expression network, an edge between two genes (nodes) represents a co-expression relationship. For each dataset or module a network interaction file was created from its adjacency matrix, and used in Cytoscape for visualization and analysis.

Network analysis often results in a large number of modules. It is sometimes useful to reduce the number of modules by merging those whose expression profiles are very similar. This was accomplished by merging modules whose member genes were highly co-expressed. To calculate the co-expression similarity of entire modules, their module eigengenes were calculated. The module eigengene is defined as the first principal component of a given module. It can be considered as a representative of the gene expression profiles in a module. The module eigengenes were clustered on their consensus correlation, which was the minimum correlation across the two sets.

A list of genes belonging to each network module was exported to tab delimited text files along with all necessary information. For each module there were two files, the first file contained a list of genes with their gene symbols, mean expression, module names, and intra-modular connectivity. This file was used for GO analysis using The Database for Annotation, Visualization and Integrated Discovery (DAVID)² (Huang et al., 2007, 2009a,b). The second file contained co-expression interaction information between each pair of genes in a module along with the topological overlap and correlation information. This interaction file was used for network visualization and analysis.

GO functional Annotation Clustering analysis was performed through DAVID web-services using the gene list for each young network module. In this research, DAVID web-services were accessed programmatically by using an R package called RDAVIDWebService (Fresno and Fernández, 2013). Since gene symbols can be confusing and often fail to produce a perfect match, the corresponding affymetrix IDs were used to query the DAVID database. GO functional annotation information was obtained for all modules in the young and the aged categories.

In DAVID, for each functional cluster an enrichment score is calculated. This enrichment score is the geometric mean (in −log scale) of the p-values of all member annotation terms and is used to rank their biological significance (Huang et al., 2009b). Thus, the top ranked annotation clusters will most likely have consistently lower p-values for their annotation members. The significance of a gene-enrichment p-values for each annotation term is calculated based on a modified Fisher exact test method known as the EASE score (Hosack et al., 2003). The default threshold of the EASE score was set at 0.1.

Network modules for young and aged were compared across studies and platforms for their repeatability using the statistics implemented in WGCNA software package. This was done in two ways: (a) module preservation; and (b) module overlap.

Differential network analysis allows one to compare two different networks side by side, for example, between a control and a disease network. Networks for several interesting modules identified in this research were visualized side by side between the young and aged groups using Cytoscape and compared for their differential co-expression.

Top hub genes were identified by using module eigengene-based connectivity or k_ME values in both the young (learning unimpaired) and aged (predominantly learning impaired) networks. Module eigengene-based connectivity k_ME, also known as module membership, is calculated for each gene. It is defined by correlating each gene’s expression profile with the module eigengene of a given module (Zhang and Horvath, 2005; Langfelder and Horvath, 2008; Miller et al., 2010). Hub genes were validated by assessing their repeatability in networks constructed from independent datasets and by investigating their functions in relevant pathways. In addition, expression patterns of selected hub genes were verified using meta-analysis and forest plots.

Repeatability of the candidate hub genes were assessed as follows. For each module, hub genes identified in the exploratory (R7) networks were checked if they are also identified as hub genes in the validation networks (e.g., B8, K9, or B7) with high k_ME values as well as with t-test p-values ≤ 0.05 (between two networks at a time e.g., one exploratory and one validation). In cases where a module from an exploratory network matched to multiple modules in a validation network, genes from multiple significant modules in the validation network were combined together and then compared to the hub genes in the exploratory network module.

Literature searches were performed using PubMed to explore characteristics and functions of selected ASLI candidate hub genes and their relationship to learning and memory formation.

Based on the quality of data and number of samples (Supplementary Table S1), R7 aged (R7-A) and young (R7-Y) datasets were chosen as the exploratory datasets; B8 young (B8-Y), K9 young (K9-Y), B7 aged (B7-A), and B8 aged (B8-A) datasets were chosen as the validation datasets (Table 1). After preprocessing, other groups did not have sufficient number of samples for WGCNA, they were excluded from this analysis. The networks were constructed for each of the aged and young datasets separately (i.e., B7-A, B8-Y, B8-A, K9-Y, R7-Y, and R7-A). However, GO based functional analysis and visualization was done only for the networks from R7 young and aged exploratory datasets, and the results were validated independently in networks constructed from the validation datasets. All datasets combined, this research included 65 young rat samples that were identified as learning unimpaired, 66 aged samples that were identified as learning impaired, and 29 aged samples that were identified as learning unimpaired.

Based on the scale free topology model fit analysis the soft-threshold power for both R7-Y and R7-A was determined to be 6. This power also results in an approximate straight line relationship in the scale-free topology plots (Supplementary Figures S8, S9). Performing similar analyses, the soft powers for B8-Y and B8-A were determined to be 10 and 8, respectively. For the B7-A dataset, the soft power was 9 and for K9-Y it was 10.

The genes that remained after preprocessing (Supplementary Table S2) were used to calculate the signed Pearson correlation coefficients for all pairwise comparisons of gene-expression values across all young and aged samples. The correlation matrix for each group was then transformed into a matrix of connection strengths (i.e., an “adjacency” matrix) using a soft power beta as determined above. This resulted in a network adjacency matrix for each dataset, for example, for R7 it generated an 8053 × 8053 matrix.

First, connectivity value for each gene was calculated from the adjacency matrix. Next for each dataset, the average median connectivity k_med was used as a cut-off value to filter out genes with very low connectivity. For R7-Y k_med was 0.46 and for R7-A k_med was 0.54. We selected the average k_med = 0.5 as the minimum connectivity cut-off, which removed 2379 genes, leaving 5674 high connectivity genes for the R7 network analysis (Supplementary Table S2). For B8 and K9 the median k_med was 0.4 and 0.35, which resulted in 5202 and 4796 high connectivity genes, respectively. The number of B7 genes was already low and close to the numbers of other filtered datasets. So, in order to prevent information loss no filtering was done on these B7 genes.

A major goal of gene correlation network analysis is to identify groups of highly interconnected genes (Zhang and Horvath, 2005; Oldham et al., 2006) termed as modules. The expression profiles of genes in a module are highly correlated across the samples. In a co-expression network, modules are identified by searching for genes with similar patterns of connection strengths to other genes, or genes with high topological overlap. The topological overlap values are calculated using the adjacency and connectivity values, which determine which genes will be in which module and form a network. The values range between 1 and 0 representing maximum and minimum interconnectedness. The module identification method in WGCNA is based on using a node dissimilarity measure in conjunction with a clustering method. Since the topological overlap matrix is non-negative and symmetric, it is turned into a dissimilarity measure by subtracting from 1. Genes are hierarchically clustered using the average linkage method, taking 1-topological overlap as the distance measure and modules are determined by choosing a height cutoff for the resulting dendrogram. In the dendrogram, discrete branches of the tree correspond to modules of co-expressed genes. Following these steps, gene network modules for the young and aged samples were identified separately for each dataset using the filtered weighted correlation matrices as prepared above.

Figures 1, 2 show the hierarchical dendrograms of topological overlaps for the 5674 genes in R7-Y and R7-A, respectively. There are several height cut-off algorithms implemented in the WGCNA R package. In this research, the cut-tree hybrid method was chosen to pick a height cut-off and to identify modules, which are shown in the panel below the dendrograms. The default lowest cut-off resulted in six modules in the young network and 15 modules in the aged network. Each module is labeled with a unique color (except gray) for easy visualization and understanding. The color gray is preserved for genes that do not belong to any module.

The aged network resulted in many modules, most with small numbers of genes, for example, 13 of the modules had fewer than 300 genes each and nine of them had less than 200 genes each (result not shown). For better comparison, the number of modules in the aged network was brought closer to that of the young network. This was accomplished by merging the modules (Supplementary Figure S10). In order to keep the module numbers similar to that of the young network, a cut height of 0.4 was chosen that generated seven modules in the aged network (including the gray module; Figure 3).

Since module names/labels in a network were randomly generated, the seven aged modules were matched to the seven young modules to check for similarity and module overlap of gene members (Supplementary Tables S3, S4). Once a significant match was found, modules in the aged network were renamed after the matched young network module names. Table 2 shows the final modules in the young and aged networks along with the number of genes belonging to each module. In addition, Table 2 shows which aged modules are matched to which young modules. The black module from the aged network had genes matching significantly to both the blue and brown modules in the young network. The aged brown, red and cyan modules matched to the green, red and yellow young modules, respectively, while the blue and pink aged modules matched a single turquoise young module. This module matching process is helpful when comparing similar modules between networks, for example, aged vs. young.

For clarity, only the top 500–600 most connected genes and their co-expression interactions were used to create each module network. Co-expression information from all modules were combined and imported into the Cytoscape for visualization. Figure 4 shows all six modules in the R7 young networks where the modules are represented by the color of their respective names (e.g., the blue module is represented by the color blue).

Biological significance analysis of the network modules was performed using the functional annotation clustering analysis in DAVID that utilizes the GO and other biological pathway information databases. DAVID functional annotation clustering analysis was used through the RDAVIDWebService tool in R to identify the most relevant (overrepresented) biological terms associated with each module gene list. The DAVID database offers extended annotation coverage with over 40 annotation categories, including GO terms, protein-protein interactions, protein functional domains, disease associations, bio-pathways, sequence features, homology and many more (Huang et al., 2009b). However, for reasons of simplicity and to better understand the biological significance of the network modules identified above, only the biological processes (BP), molecular functions (MF), and cellular components (CC) GO terms and all KEGG Pathway terms were included in the functional annotation clustering analysis.

Affymetrix probe set identifiers of all the genes belonging to a network module (Table 2) were used as the input gene list. The total number of genes from the RAE230A array for the R7 dataset (after preprocessing and filtering) was 5674, and was used as a background population. Rattus norvegicus was used as species. The function getClusterReportFile(…) in RDAVIDWebService was used with default parameters to retrieve all relevant information. Next getClusterReport(. . .) function was used to extract the functional annotation chart file, which was saved as a text file and later analyzed. An enrichment score cutoff of 1.0 was used to minimize the number of clusters that were returned.

Table 3 shows the summary result of GO analysis for the young modules. The most significant GO functional categories represented by the genes belonging to each module are also shown in Figure 3. The results show that, in general, each module is highly enriched with genes functioning in broad but distinct GO functional categories or biological pathways with highly significant enrichment scores.

Network modules were validated by assessing their preservation and overlap across datasets. This was done by comparing the modules’ gene expression data as follows: R7 young vs. R7 aged; R7 young vs. B8 and K9 young; R7 aged vs. B8 and B7 aged.

In order to assess the changes in co-expression patterns of the young as they age and how the aging would affect learning impairments, we compared several interesting network modules between young (learning unimpaired) and aged (predominantly learning impaired) networks generated from the R7 data. This comparative investigation involved visualizing them side by side, comparing expression patterns between networks, and searching for key genes. In addition, it involved identifying the key genes’ functions and pathways that can help explain the learning differences as well as the aging effect that had been observed between the young and aged animals. Differential expression levels for the top 5674 genes in the R7 data were calculated by using the limma package in Bioconductor. The log fold changes of expression differences between young and aged for all genes were saved as a tab delimited text file, and later loaded as node attributes in Cytoscape for each module.

The module that is most relevant to this article is the yellow (“learning and memory”) module. Figures 8A,B presents the differential co-expression networks of this module between young (A) and aged (B) rats, which demonstrates a clear difference in expression patterns between the young and the aged genes. The majority of the genes in the aged yellow network show lower expression compared to the young. In addition, the comparative analysis demonstrates differential co-expression for many genes between the two networks (i.e., some genes display more co-expression interaction than others and this varies between the young and the aged networks). The results allow one to identify a number of key ASLI genes for further investigation (see below).

In a co-expression network, genes that are highly connected with many other genes are called hub genes. These genes show significant correlation with the module eigengenes and have high within-module connectivity. After closely studying the networks in young (learning unimpaired) and aged (predominantly learning impaired), we have identified a set of key hub genes in each module. Some of the hub genes in the yellow module in R7 are shown in Table 4. Some of them are already known as learning genes and were identified in our previous meta-analysis (Uddin and Singh, 2013). Table 5 shows the number of significant AY (aged vs. young) meta-analysis genes that are also members of different modules in the R7-Y network. Particularly, it shows that there are 165 AY significant meta-analysis genes in the yellow module. Effect size estimates from the meta-analysis for the above ASLI candidate hub genes are summarized in Supplementary Table S5. In addition, we have created individual forest plots for some of these hub genes, which are presented in Supplementary Figures S11–S25.

The candidate ASLI hub genes were checked for their repeatability in networks constructed independently from B8, K9, and B7. The results are summarized in Table 6. Details of the hub gene validation data are available in Supplementary Tables S6–S11. The results show that a number of hub genes from the yellow module are repeated in one or more independent datasets in B8, K9 or B7 with a p-value ≤ 0.05. From the R7 yellow module Prkacb, Scn2b, Cntn1, Pten and Ndfip2 were found present as hub genes in the K9 network; Dlgap1 was found in the B7 and B8 networks; and Camk1g was found repeated in the B7 network. Notably, many of these hub genes were in the list of top 20 mean KME values in other networks, but their p-values were not significant, for example, Dlg3, Mapre1, Dpp6, Stxbp1, Impact and Mapk1.

One useful property of a co-expression network is module. In a module the expression patterns of the genes are mutually correlated (Langfelder and Horvath, 2008). The focus on co-expression modules, each consisting of possibly hundreds of genes with common co-expression across samples, allows for a biologically motivated reduction of data while also alleviating the problem of multiple comparisons (Levine et al., 2013). Further, just as correlated genes tend to have similar biological functions, on a larger scale, modules tend to contain genes with similar biological functions (Lee et al., 2004).The results obtained in this research and the follow up network analysis support these hypotheses. For example, the use of WGCNA reduced R7 data into a few biologically meaningful co-expression modules. The follow up GO analysis and literature search results were persuasive enough to indicate that each module gene set likely serve a distinct major biological function, thus, pointing to the widely held notion of “co-expression to co-functionality”. It is important to note that the networks and modules constructed from R7 microarray data were based on the gene expression patterns alone (i.e., there was no prior knowledge of the genes’ function at the time of network construction). Once the networks were divided into modules and their module-wise GO functional analysis was performed, it was indeed observed that each module pointed to a broad but distinct category of biological function, and genes in each module shared similar subcategories of functions all converging to the broad functional category of the module (Table 3, Figure 4). Particularly, the genes in the yellow module showed significant enrichment in GO functions and pathways related to learning and memory formation in the brain. This is expected as the young and aged rats used in the original research were tested for their memory performances where the experimental young animals displayed clear learning unimpairment and the majority of the aged animals demonstrated learning impairment. The inclusion of various controls, including the aged animals that demonstrated learning unimpairment, served as controls e.g., for stress, physical ability as well as for factors that may contribute to learning irrespective of age. It is likely that the learning unimpaired aged rats might also display sign of learning impairment as they progress towards further aging. Although, the other modules are enriched with functions not directly related to learning and memory, they are critical for normal neuronal processes such as communication, growth, development and maintenance. For example, genes in the brown module are significantly enriched in functions contributing to the various cellular processes and communication, the green module genes in developmental processes, and the red module genes in oligodendrocyte development.

Thus, alteration of these modules’ normal module-wise functions at old age through altered gene expression, as observed in the datasets, has the potential to affect normal functioning of learning and memory formation process. Preservation of these modules were not only validated across networks created from independent datasets, but also the gene members of these modules demonstrated significant module membership (module overlap) across the independent networks (Figures 6–8).

Gene co-expression analysis studies in multiples species, tissues and platforms have shown that co-expressed genes tend to be functionally related (Williams and Bowles, 2004; Obayashi et al., 2008; Oldham et al., 2008). In order to investigate, whether observed clusters or modules of co-expressed genes are of functional significance, Lee and Sonnhammer (2003) observed that genes involved in the same biochemical pathways tend to be clustered together in a number of eukaryotic genomes. By a heuristic generalization known as “guilt by association”, it has been computationally established that functionally related genes are organized into co-expression networks, in practice assisting functional annotation of uncharacterized genes (Michalak, 2008). For example, physically interacting proteins in yeast were found to be encoded by co-expressed genes (Ge et al., 2001; Wuchty et al., 2006). These observations likely have inspired the development of co-expression network analysis methods. Gene network modeling using co-expression approaches provide insight into cellular activity as genes that are co-expressed often share common functions (Piro et al., 2011). Such networks have been widely used to study many diseases and phenotypes because of their ease of use and their ability to provide more biologically meaningful results (Gargalovic et al., 2006; Chen et al., 2008; Min et al., 2012; Zhou et al., 2014; Holtman et al., 2015; Maschietto et al., 2015; Rickabaugh et al., 2015; Spiers et al., 2015; Ye and Liu, 2015).

Microarray data captures functional relationship among genes that can provide biologically relevant information. In traditional microarray data analysis, however, these relationships remain essentially unexplored. Thus, a modular approach to gene function through WGCNA provides a sensible way to extract such functional information from large microarray datasets in a biologically meaningful way. Particularly, this analysis has shown that specific learning associated functional gene modules can be identified through co-expression network modeling where genes in the module show significant enrichment in learning and synaptic plasticity related GO functions.

Hub genes play a central role in the structure of co-expression networks as they are often relevant to the function of regulatory networks. The ability to efficiently transit cellular signals within and between co-expressed clusters is facilitated by “hubs”, which are connected to a large number of nodes (Gaiteri et al., 2014). Analysis of the yeast protein-protein interaction network revealed that highly connected nodes are more likely to be essential for survival (Jeong et al., 2000; Carter et al., 2004; Han et al., 2004). The co-expression networks of the yellow “learning and memory” module (Figure 8) display a tight interrelationship of a large number of nodes with some hub genes. What is most interesting is that the co-expression of these hubs and nodes, as demonstrated in this network analysis, is not a random aggregation of some genes. Literature review suggests that the correlated expression pattern of the hub genes in the yellow networks (Figure 8A) may in fact be highly coordinated, and inside the young rats’ hippocampus they may be serving a common purpose. The combined effect of the functions of the hub genes that are co-expressing together in individual modules may in fact contribute to the co-functionality of the whole module. The purpose could be to maintain the functional integrity of the normal process of learning and memory formation mechanisms, which are disrupted in the aging brain. Indeed, the side-by-side comparison (Figure 8) of the yellow module networks between young and aged rats demonstrates a clear difference in expression patterns. The majority of the genes in the aged yellow network show lower expression compared to the young. In addition, the comparative analysis demonstrates differential co-expression for many genes between the two networks. These genes display more co-expression interaction than others and the number of interaction varies between the young and the aged networks. For example, the gene Dlg3 has 63 co-expression connections with other genes in the young network while only one in the aged. Similarly, Dpp6 has only two co-expression connections with other genes in the young network while 68 in the aged (Table 6). Since, such genes with many co-expression connections with other genes in a network often play a role as hubs; we have short listed 19 of these genes as candidate ASLI hub genes from both the young and aged yellow module networks. These genes include Camk1g, Cdk5r1, Cntn1, Dlg3, Dlgap1, Dpp6, Eif5, Gabrg1, Impact, Kcnab2, Mapk1, Mapre1, Ndfip2, Ppp2r2c, Prkacb, Pten, Rasgrp1, Scn2b and Stxbp1 (Table 4). The results show that some of these hub genes are already known as key learning and memory genes and have well established roles in memory functions. While for others, information is emerging indicating their direct or indirect role in learning and memory. Below we summarize what is already known from the literature about the molecular mechanisms of learning and memory formation and how the candidate ASLI hub genes from this research fit into that big picture.

Differential co-expression refers to changes in gene-gene correlations between two sets of phenotypically distinct samples (de la Fuente, 2010). Changes in gene-gene correlation may occur in the absence of differential expression, meaning that a gene may undergo changes in regulatory pattern that would be undetected by traditional differential expression analysis (Gaiteri et al., 2014). The fact that the altered regulatory patterns observed within tissues across phenotypic states in manners that are reflected in altered co-expression networks has been shown in aging mice (Southworth et al., 2009), across corticolimbic regions in major depression (Gaiteri et al., 2010) and between miRNA’s in Alzheimer’s disease (Bhattacharyya and Bandyopadhyay, 2013).

What becomes apparent is that the differential expression and differential co-expression may point to distinct cellular mechanisms involved in ASLI, which may be working in different ways in the cell. For example, our differential expression meta-analysis has identified a large number of genes showing significantly altered expression in the aged rats compared to young rats (Supplementary Tables S1, S2 in Uddin and Singh, 2013). These genes include many immediate early (e.g., Arc) or late phase genes (during gene expression) as well as other genes contributing to aging and ASLI. Major functions disrupted by these genes include cell viability, axonogenesis, quantity and synthesis of IP3 and formation of cells.

On the other hand differential co-expression analysis presented here identified a set of modules each with distinct functions. In addition, it has identified a set of candidate ASLI hub genes in one of those modules. From the known function of these hub genes as explained above, it is evident that many of these genes function as kinases and phosphatases in the neuronal information flow process, starting from the synaptic junctions/synapses to the nucleus to activate various transcription factors. Though scattered in different networks, meta-analysis has also identified few hub genes functioning as kinases or in ion channels. Thus the hub genes may be triggering one or more mechanisms that activate other key factors in a number of pathways, which set the stage for the expression of several immediate early or late phase genes, which again most likely activate the expression of majority of the differentially expressed genes. Learning in the young animals most likely induces such mechanisms that synchronously regulate transcription of multiple genes, and may potentially generate co-expression relationships.

Another important observation to note is that all the learning related genes identified in the differential expression and pathway analyses (and genes they generally interact with) are scattered in different networks and pathways (Appendix 6.3.1 to Appendix 6.4.4 in Uddin, 2015). In contrast, the current differential co-expression analysis identified many known learning genes (or genes that appear to be contributing to learning and memory functioning) that are highly concentrated and co-expressed in the yellow “learning and memory” module.

Interestingly, the candidate ASLI hub genes are expressed at a comparatively lower level, with small differences in expression (e.g., effect size) between young and aged samples. For example, the Prkacb hub gene was not known to be a learning gene (Supplementary Tables S3, S4 in Uddin and Singh, 2013). Like Prkacb, most of the hub genes failed to show significant effect size or differential expression values and remained undetected in the meta-analysis (Supplementary Tables S1, S2 in Uddin and Singh, 2013). This fact highlights the importance of alternate analysis like WGCNA to identify genes that are not detected using the traditional methods. Similar observations have been demonstrated in past studies (Rhinn et al., 2012). For example, the alpha synuclein gene variant “aSynL”, containing a long 3′sUTR, was identified as the most differentially coexpressed gene in several Parkinson’s disease datasets. However, aSynL was not highly differentially expressed and thus would have likely been overlooked by traditional microarray analysis (Gaiteri et al., 2014). Thus, through the identification of modules and hubs, differential co-expression analysis can be used to prioritize specific phenotype-related important molecules.

Another very interesting property of co-expression networks is the network connectivity. Our findings (Supplementary Figure S26) support the newly emerging hypothesis (Oldham et al., 2006; Miller et al., 2008) that differential connectivity is different from differential expression. During the network construction process, I selected genes with high connectivity and filtered out all low connectivity genes (Supplementary Table S2). The observation is that the resulting network modules represent a set of highly connected genes as hubs that were virtually absent in the differentially expressed top gene list (in Uddin and Singh, 2013) and vice versa. In fact, it has been reported that gene-gene correlations in disease can occur with or without changes in expression (Hudson et al., 2009). In addition, differentially expressed genes in some complex psychiatric diseases can have low connectivity, which reside on the periphery of co-expression networks for neuropsychiatric disorders such as depression, schizophrenia and bipolar disorder (Gaiteri and Sibille, 2011; Gaiteri et al., 2014).