Two Gene Co-expression Modules Differentiate Psychotics and Controls

Schizophrenia and bipolar disorder are highly heritable psychiatric disorders. Associated genetic and gene expression changes have been identified, but many have not been replicated and have unknown functions. We identified groups of genes whose expressions varied together, i.e. co-expression modules, then tested them for association with schizophrenia. Using Weighted Gene Co-expression Network Analysis, we show that two modules were differentially expressed in patients versus controls. One, up-regulated in cerebral cortex, was enriched with neuron differentiation and neuron development genes, as well as disease GWAS genetic signals; the second, altered in cerebral cortex and cerebellum, was enriched with genes involved in neuron protection functions. The findings were preserved in five expression data sets, including sets from three brain regions, from a different microarray platform, and from bipolar disorder patients. From those observations, we propose neuron differentiation and development pathways may be involved in etiologies of both schizophrenia and bipolar disorder, and neuron protection function participates in pathological process of the diseases.


Introduction
Schizophrenia (SCZ) and bipolar disorder (BD) are psychiatric disorders with world-wide lifetime prevalences of around 1% 1-4 . These disorders have been shown to be highly heritable by monozygotic and dizygotic twin studies and by adoption studies: the heritability estimates for SCZ range from 70 to 85% and, for BD, from 60 to 85% 5 . Over the past two decades, genetic studies of thousands of samples have identified hundreds of candidate genes [6][7][8] . Most of these findings have not been supported by genome-wide studies. Moreover, functional genomic roles have not yet been determined for most of the associated genes.
Gene expression can bridge the gap between genetic variation and disease susceptibility as an intermediate phenotype that is regulated by a combination of genetic and epigenetic factors. Gene expression transcription profiling is widely used and has been thoroughly validated by the MicroArray Quality Control (MAQC) project 9,10 . Expression microarray studies test thousands of gene transcripts for differential expression simultaneously. When all these genes are tested for association with disease, it creates a multiple testing problem: to minimize false positives, the multiple test correction sets the significance threshold so high that true positives might be missed as well. Additionally, testing individual genes for disease association also ignores the interaction between genes. So far, hundreds of gene expression microarray studies of psychiatric diseases have been reported, including studies of SCZ, BD and autism. However, the significant gene changes detected by one study are seldom replicated in another, let alone across three or more 11 .
Gene expression network analyses are an alternative approach for analyzing expression data that reduce the sample space tested (that is, the number of hypotheses to be tested). Gene expression networks are constructed from expression data from thousands of genes, and describe the interactions among groups of transcripts. They can be used to observe systematic alterations in expression associated with complex diseases such as psychiatric disorders. Horvath and colleagues have developed an algorithm for creating networks, called Weighted Gene Co-expression Network Analysis (WGCNA), which is widely used 12,13 . This method identifies groups of genes within a network whose expressions are highly correlated. These groups, called modules, can then be compared between cases and controls, among different tissues, species, or other phenotypes or clinical traits [14][15][16][17] .
Oldham et al. were the first to apply this method to expression in brain. They compared the gene expression across different brain regions and demonstrated that the modules reflect the underlying cell-type composition of the regions 18 .
The method was first applied to a psychiatric disease by Torkamani et al., who detected SCZ-associated gene co-expression modules in one microarray data analysis, and found that aging affected gene expression in normal controls, but not in SCZ patients 19 . They made the interesting distinction between constructing networks from case and control data separately and constructing a network from a combination of case and control data. Constructing modules from case and control networks separately allows comparison of basic module structure between the two groups. If modules are not substantially preserved in cases, a conclusion can be drawn that fundamental relationships among genes have been disrupted. If modules are substantially preserved in cases, the case and control data can be combined for network construction, and the detected modules assessed for differential expression between cases and controls. Torkamani found that modules were substantially preserved in cases, and identified five SCZ-associated modules; Gene Ontogeny (GO) enrichment analysis showed those five modules were associated with oxidative phosphorylation, angiogenesis, neuron differentiation, chromatin and nucleosome assembly, and inositol phosphate metabolism, separately.
Voineagu et al. applied WGCNA to another psychiatric disease, autism. In addition to a GO enrichment analysis, they found one of the differentially expressed modules was enriched in disease genome-wide association study (GWAS) signals. They interpreted this as evidence that the module's member genes were causally associated with autism. They also found another module with altered expression but not enriched in GWAS signals and took it as an indication of a non-genetic etiology 20 .
Note that all three of these network studies were based on analysis of a single microarray data set; their findings have not yet been replicated in other data sets at the network level.
To demonstrate the robustness and reproducibility of WGCNA findings, we used multiple microarray data sets to study co-expression networks in brains of psychiatric patients. We constructed gene expression networks, and then identified gene co-expression modules within the networks. The modules were tested for association with SCZ. We assessed whether the disease-associated modules were detected in independent data sets, including sets from several different brain collections and different brain regions. Reproducible, or preserved, modules were then evaluated for case-control differences in the independent data sets.
Since it is widely accepted that SCZ and BD have genetic factors in common 21,22 , we tested whether the gene modules perturbed in SCZ were similarly perturbed in BD. We also tested whether the disease-associated modules were enriched with genetic variants that had been previously associated with disease by GWASs 23,24 .

Samples and quality control
Cerebellum (CB) and parietal cortex (PCX) brain tissues were obtained from Stanley Medical Research Institute (SMRI) 25 ,. They came from the SMRI's Neuropathology Consortium and Array collections, and included 50 SCZ samples, 50 BD samples, and 50 unaffected control samples. Expression data for these samples came from the NCBI Gene Expression Omnibus (GEO) database (GSE35978). One of the two prefrontal cortex (PFC) brain expression data sets came from Dobrin's group using SMRI samples (PFC-SMRI) 25 , while the second came from the Victorian Brain Bank Network (PFC-VBBN) 26 and was obtained from at GEO (GSE21138, sample information and preparation in Supplementary File).
ComBat, a batch effects adjustment program that we have previously shown to be the best available, was used to remove batch effects from the data sets (see Supplementary File for detailed preprocessing steps) 27,28 .

Gene network construction and module detection
We used weighted gene co-expression network analysis (WGCNA) 13 to identify modules of co-expressed genes within gene expression networks. To construct the network, the absolute values of Pearson correlation coefficients were calculated for all possible gene pairs. Values were entered into a matrix, and the data were transformed so that the matrix followed an approximate scale-free topology (see Supplementary File for detailed information). A dynamic tree cut algorithm was used to detect network modules 29 . WGCNA and the dynamic tree cut algorithm were implemented in R 12,29 .
We ran singular value decomposition (SVD) on each module's expression matrix and used the resulting module eigengene, which is equivalent to the first principal component 30 , to represent the overall expression profiles of the modules.

Module preservation statistics
We utilized the module preservation statistic Z summary , to assess the module preservation from different expression data sets 33 (see Supplementary File for formal definition). Unlike the cross-tabulation test, Z summary not only takes into account the overlap in module membership, but also the density and connectivity patterns of modules. In addition, for our study, network-based preservation statistics only require that module membership be identified in the original data sets, reducing the variation coming from various parameters setting to identify new modules in validation data sets. We converted the probe-level measurements into gene-level measurements to make data from different platforms comparable. The probe within a gene that had the highest coefficient of variation was used to represent that gene. Overall, 8497 genes were retained in our preservation calculation.

Differential expression test for modules
We used multiple linear regressions on the modules' eigengenes to remove the effects of sex, age, pH and PMI from the SMRI samples and VBBN samples. The residual eigengene values were then used to test the disease association using Pearson's correlation test. We used false discovery rate (FDR) for multiple testing correction 49 .
The TGen-BD GWAS data came from the genome-wide association study of the Translational Genomics Research Institute's (TGen) BD data (http://www.tgen.org/), which comprised 1,190 BD cases and 401 controls. Sample genotyping was conducted using the Affymetrix GeneChip Mapping 5.0K Array. We performed imputation using MaCH v1.0 to increase the density of interrogated SNPs 31 , with HapMap data as reference. Overall, 2,593,107 SNPs in GAIN-SCZ, 3,281,319 SNPs in GAIN-BD, and 2,542,706 SNPs in TGen-BD were included after imputation.
Imputed GWAS data was used to test whether the two disease-associated modules were enriched in SCZ/BD association signals 32 . We used a previously-reported procedure to run the enrichment test 32 . First, the max −log(P-value) of a SNP located between 20kb upstream and downstream of a gene was assigned to represent the gene, then the gene set's enrichment scores (ES) were calculated based the gene's rank. SNP level permutation was applied to generate the distribution of the ES and then the distribution was normalized. For multiple gene sets, FDR was calculated by joining all the distributions of ESs, each for one gene set, generated by permutation. As the difference of gene length distribution between genes in one module was significant (p=5.47e-27), we also applied a permutation procedure to verify whether there was a gene length bias in the genetic signals enrichment test 43,50 , and there was no bias (see Supplementary Methods section 8).

Gene modules in parietal cortex
We first analyzed the parietal cortex (PCX) brain gene expression in 50  Networks were constructed in two different ways: first, we constructed one network from control data and one network from case data, then identified modules in the control network and assessed their preservation in the case network. Second, we constructed a network from the combined case and control data, and identified modules within it.

Structure of co-expression modules in schizophrenia cases and controls
Case and control sample gene networks were constructed separately to detect whether there was any gene co-regulation disruption or creation in cases relative to controls. We assessed module preservation using a permutation-based preservation statistic, Z summary , developed by Langfelder et al., which assesses whether the connectivity level and pattern of a module in one data set is preserved in another 33 . The authors suggest the following significance thresholds: Z summary < 2 implies no evidence for module preservation, 2 < Z summary < 10 implies weak to moderate evidence, and Z summary > 10 implies strong evidence for module preservation. In our study, all modules detected in the control data had Z summary greater than 10 in the case data ( Fig. 1), suggesting well-preserved membership and connectivity of the control modules in the SCZ cases. This is consistent with two previous gene network-based psychosis studies 20,34 , where the case modules had no obvious perturbations relative to control modules.
To test the reliability of the module construction results, we compared the modules identified from control samples to published gene expression networks. Cross-tabulation showed that our control modules have similar module membership to modules previously reported by Oldham et al. 18 (Supplementary Fig. 1). Slight differences, including some modules in our controls collapsing into one module in the Oldham data, might have been due to the differences in samples and platforms, data-filtering steps and/or parameters used for network construction and module detection.

Differential expression of modules in schizophrenia patients
Since there was no significant difference in module structure between cases and controls, modules identified in a gene network constructed from both case and control data were analyzed for differential expression in SCZ patients. Twenty-four modules were detected ( Supplementary Fig. 2). We used an eigengene to summarize each module's expression profile 30 (see Methods for formal eigengene definition). After multiple linear regressions on the eigengenes to remove the effects of sex, pH and post-mortem interval (PMI), three covariates that can affect measures of gene expression 11 , the residual module eigengenes for each individual were tested for disease association. The eigengenes of two PCX modules, referred to as M1A and M3A, were significantly associated with SCZ after multiple testing correction ( Fig. 2A). M1A comprised 490 genes, while M3A comprised 106. NOTCH2 and MT1X, respectively, were the hub genes of M1A and M3A, meaning that, of all the genes in their modules, these two genes had the strongest correlation with the module's eigengene, as well as being highly connected to the other genes in their modules (Supplementary Tables 5  and 6).

Replication of module structure
We tested co-expression module preservation across different microarray data sets. We again used the module preservation statistic, Z summary , to compare modules. Data sets included prefrontal cortex (PFC) tissues and cerebellum (CB) tissues, also from the Stanley Medical Research Institute (SMRI) but performed by a different research group 25 , and another PFC data set, where samples came from the Victorian Brain Bank Network (VBBN) 26 .
As measured by the Z summary statistic, six of the 18 modules identified in the SMRI PCX data were strongly preserved in the SMRI PFC data set, including M1A. Six modules were moderately preserved, including M3A, and four modules were not preserved ( Supplementary Fig. 3). Similarly, M1A was strongly preserved and M3A moderately preserved in the SMRI CB and VBBN PFC data sets (Table 1; Supplementary Fig. 4 and 5).

Replication of disease association
We tested the modules' hub genes for disease association, since hub genes by definition are highly correlated with their modules' eigengenes 35 . After controlling for the effects of covariates, NOTCH2, the hub gene of M1A, was significantly associated with SCZ in VBBN and SMRI PFC data (p=0.049 and 0.022, respectively), but not in SMRI CB data (p=0.583). MT1X, the hub gene of M3A, was consistently upregulated in SCZ in VBBN PFC, SMRI PFC and SMRI CB data (p=3.5E-03, 0.026 and 1.3E-04, respectively) (Fig. 2B,  2C). These results suggest that M1A gene network perturbation between SCZ and controls occurs in cerebral cortex but not in cerebellum, while M3A perturbation happens across all brain regions.

Modules shared between schizophrenia and bipolar disorder
We were interested in whether the M1A and M3A modules associated with SCZ were also associated with BD. One BD data set with two regions studied was tested, PCX and CB from SMRI. The control samples here overlapped with the control samples used to construct the SCZ case and control network and modules. We first confirmed that the M1A and M3A in PCX were well-preserved in the BD case-control data sets (Table 1; Supplementary Fig. 6 and 7), then we tested whether NOTCH2 and MT1X gene expression were associated with BD. NOTCH2 showed an association trend in PCX (p=0.054) but not in CB (p=0.117), while MT1X showed significant association in both regions (p=0.015 in PCX and 1.5e-3 in CB) (Fig. 2B, 2C), indicating M1A in PCX and M3A in two brain regions were altered in BD, and were shared by two diseases.

Characterization of two disease-associated modules
NOTCH2, a member of the Notch gene family, had the highest intramodular connectivity in the SMRI PCX M1A module. The Notch signaling pathway plays a key role in cell to cell communication, and was reported to play multiple roles in central nervous system development [36][37][38] . Functional enrichment analyses were performed with The Database for Annotation, Visualization and Integrated Discovery (DAVID) (Fig. 3, Supplementary Table  7) 39 . Consistent with the hub gene's function, M1A was highly significantly enriched with gene ontology terms related to neuron differentiation and neuron development genes (p=7.50e-08 and 3.63e-06, FDR q=1.27e-04 and 6.14e-03, respectively).
We also found that 46 genes in M1A overlapped with a synaptic gene group associated with SCZ identified by Lips et al. 40 . These overlapping genes are interesting because ten of them have previously been reported to be associated with either SCZ or BD, which may suggest a relationship between synaptic dysfunction and these two psychiatric diseases. These genes include ABLIM1, APOE, AQP4, SLC1A2, SLC1A3, SLC4A4, GABRG1, NTRK2, ADD3 and MAOA.
Functional enrichment analysis showed M3A to be highly enriched in gene ontology categories related to metallothioneins (MT) and metal binding site functions (p=2.41e-05, FDR q=0.031) (Fig. 4, Supplementary Table 8). MT genes respond to several stimuli, including metals, oxidative agents, inflammation and stress; they influence cognition, protect against neurotoxicity and play a role in the astrocytic response to CNS injuries 41 . The expressions of our MT gene family members, including the hub gene MT1X, and MT1E, MT1F and MT2A have previously been reported to be differentially expressed in SCZ 42 .

Enrichment of genetic association signals
To determine whether the expression alterations to modules M1A and M3A had genetic bases, we tested the modules for enrichment with SCZ and BD genetic association signals. M1A showed significant enrichment of signals from both SCZ and BD GWASs (GAIN, SCZ, p<0.001, FDR q=0.009; GAIN-BD, p<0.001, FDR q=0.0015; TGen, BD, p<0.001, FDR q=0.002) (Supplementary Fig. 8, 9 and 10) (see Supplementary materials for further descriptions of GWAS data sets). Those significances were retained after controlling for gene length bias 43 . M3A was not enriched in signals from either disease. We also used other GWAS data, from non-psychiatric disease type 2 diabetes (T2D), as a negative control for our enrichment analysis. Neither M1A nor M3A were enriched with signals from a T2D GWAS.

Discussion
We looked for case-control differences at the level of gene co-expression regulation, rather than the level of individual genes. We found strong evidence of such differences in one expression data set, and then demonstrated robust gene network preservation and disease association across five additional different data sets. These data sets came from different brain banks and different brain regions, and were produced on different expression microarray platforms; they contained data from patients with SCZ and BD, as well as normal controls. To the best of our knowledge, this is the first time gene expression alterations associated with disease have been replicated across so many gene expression data sets. This reproducibility suggests that our results are reliable and that the WGCNA network-based method increased the statistical power of our expression study.
M1A was enriched with neuron development and neuron differentiation genes, among others. Using module preservation statistics, we determined that the M1A module was highly preserved in the five other data sets, and was consistently up-regulated in SCZ and BD cerebral cortex relative to controls, but not in cerebellum. This result is consistent with Torkamani et al's finding that modules enriched in neuron differentiation are differentially expressed between SCZ patients and normal controls 19 . A test for enrichment with GWAS signals suggested that M1A expression changes may have a genetic basis. Based on those observations, we speculate that the etiology of psychiatric diseases could be related to developmental dysfunction in the cerebral cortex, extending the neurodevelopmental hypothesis of SCZ 44 , which has recently received further support from imaging studies 45 and deep sequencing analysis 46 .
M3A was enriched with metallothionein genes. It was moderately conserved in other data sets and consistently showed differential expression in both cerebral cortex and cerebellum between SCZ and BD patients, relative to controls. M3A was not significantly enriched in GWAS signals. Expression changes in these MT genes in SCZ and BD has been previously reported 42,47 , but, so far, no positive genetic associations of these genes have been reported. Considering that the major physiological function of MT genes is binding of heavy metals to protect the central nervous system (CNS) against neurotoxicity or other injuries 41 , the change in the expression of M3A may be a downstream effect of disease or may be related to environmental insult. Alternatively, loci regulating M3A or its member genes may exist, but not have been detected by GWAS.
One advantage of module-based association studies is that they substantially reduce the multiple testing correction burden inherent in genome-wide association studies 51 and microarray studies, particularly in brain expression studies where sample size is always limited by tissue availability. Microarray studies test tens of thousands of genes at one time, which require that p-values be very small to reach genome-wide significance. However, WGCNA produced tens of modules of highly correlated genes to be tested for association. In our case, we found M1A to be significant in most of our data sets, but when we ran an individual gene-based analysis, the hub gene of M1A, NOTCH2, was nominally significant but did not survive the multiple testing correction (data not shown). WGCNA may also reveal the function of an otherwise poorly characterized gene, if that gene is a member of a module highly enriched in a particular function. In addition, by comparing results from different data sets, we demonstrated that findings at the network/pathway level are consistent.
Our results suggested that the observed similarities between SCZ and BD are manifested at the transcriptional level. Shared genetic elements have previously been demonstrated 21,22,48 , and some expression studies reported overlapping aberrant individual genes at the transcriptional level 47 . Our study explored the shared pathology defect at a higher level, revealing a larger scale gene pathway-based perturbation. Further functional work can focus on those genes, especially the hub genes, which may reveal the molecular bases of the symptoms shared between these two diseases.
In summary, our study revealed that 1) SCZ and BD shared altered expression of neuron differentiation and neuron development genes (M1A) and metallothionein genes (M3A), which may indicate common etiology or pathology; 2) M1A's up-regulation in cerebral cortex could be regulated by genetic variants already found to be associated with SCZ and BD risk; 3) M3A expression change was found in both cerebral cortex and cerebellum, and was not detectably regulated by genetic variants.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.

Figure 1. Gene network modules from PCX from schizophrenia patients and normal controls are well preserved
Z summary is the summary preservation statistics, using the control modules as reference modules. Y-axis represents preservation statistics for the corresponding module in the case data sets, and x-axis is the gene numbers in each module. The dashed blue and green lines indicate the thresholds Z=2 and Z=10, respectively. Z summary < 2 implies no evidence for module preservation, 2 < Z summary < 10 implies weak to moderate evidence, and Z summary > 10 implies strong evidence for module preservation.

Figure 2. Module and hub genes' association test results
Eigengene-based test detected 24 modules in PCX (A), listed on the x-axis. The y-axis indicates the -log10 of association p value. The red line represents the p=0.05 threshold. M3A and M1A modules were significant after multiple test correction, with FDR q value<0.05 (green bars with red stars on the top). NOTCH2 in M1A (B) and MT1X in M3A (C), with disease in replicate data sets. Green bars represent the significance of association in SCZ and control networks, while red bars represent the significance of association in BD and control networks. The purple line marks the significance threshold of p=0.05.