Two gene co-expression modules differentiate psychotics and controls.

Schizophrenia (SCZ) and bipolar disorder (BD) 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, that is co-expression modules, then tested them for association with SCZ. Using weighted gene co-expression network analysis, we show that two modules were differentially expressed in patients versus controls. One, upregulated in cerebral cortex, was enriched with neuron differentiation and neuron development genes, as well as disease genome-wide association study 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 BD patients. From those observations, we propose neuron differentiation and development pathways may be involved in etiologies of both SCZ and BD, 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-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-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].

Subjects and Methods

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, to assess the module preservation from different expression data sets[33] (see Supplementary File for formal definition). Unlike the cross-tabulation test, Z 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].

Module-based GWAS signal enrichment test

GAIN-SCZ was the genome-wide association study of the GAIN SCZ data set, which comprised 4591 cases and controls (1217 European-American cases, 1442 European-American controls, 953 African-American cases, 979 African-American controls). The data was downloaded from dbGaP (http://www.ncbi.nlm.nih.gov/projects/gap/cgibin/study.cgi?study_id=phs000021.v3.p2). GAIN-BD was a genome-wide association study of the GAIN BD data set, which comprised 3261 cases and controls (1079 European-American cases, 1081 European-American controls, 415 African-American cases, 686 African-American controls). The data was downloaded from dbGaP (http://www.ncbi.nlm.nih.gov/projects/gap/cgibin/study.cgi?study_id=phs000017.v3.p1). Whole genome genotyping of GAIN data was done with the Affymetrix Genome-Wide Human SNP Array 6.0.

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).

Results

Gene modules in parietal cortex

We first analyzed the parietal cortex (PCX) brain gene expression in 50 SCZ patients and 50 normal controls from the Stanley Medical Research Institute (SMRI) using the Affymetrix Human Gene 1.0 ST Array[25]. After a series of sample and array-level quality control measures (see Supplementary methods for details), we retained 45 SCZ patients and 46 normal controls with measures of 19,884 transcripts.

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, 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 < 2 implies no evidence for module preservation, 2 < Z < 10 implies weak to moderate evidence, and Z > 10 implies strong evidence for module preservation. In our study, all modules detected in the control data had Z 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.

Figure 1

Gene network modules from PCX from schizophrenia patients and normal controls are well preserved

Z 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 < 2 implies no evidence for module preservation, 2 < Z < 10 implies weak to moderate evidence, and Z > 10 implies strong evidence for module preservation.

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).

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.

Replication of findings in different data sets

Replication of module structure

We tested co-expression module preservation across different microarray data sets. We again used the module preservation statistic, Z, 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 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).

Table 1

Preservation statistic of SMRI-PCX-SZ M1A (top) and M3A (bottom) on VBBN-PFC-SZ, SMRI-PFC-SZ, SMRI-CB-SZ, SMRI-PFC-BD, SMRI-CB-BD data sets.

Test Modules	Brain banks	Platform	Diseases	Brain regions	Zsummary [*]
M1A	VBBN	Affy U133	Schizophrenia	PFC	17
	SMRI	Affy U133	Schizophrenia	PFC	21
	SMRI	Affy HG 1.0	Schizophrenia	CB	29
	SMRI	Affy HG 1.0	Bipolar disorder	PCX	29
	SMRI	Affy HG 1.0	Bipolar disorder	CB	27
M3A	VBBN	Affy U133	Schizophrenia	PFC	3
	SMRI	Affy U133	Schizophrenia	PFC	6
	SMRI	Affy HG 1.0	Schizophrenia	CB	8
	SMRI	Affy HG 1.0	Bipolar disorder	PCX	7
	SMRI	Affy HG 1.0	Bipolar disorder	CB	9

VBBN, victorian brain bank network; SMRI, Stanley medical research institute; PFC, prefrontal cortex; CB, cerebellum; PCX, parietal cortex.

Z < 2 implies no evidence for module preservation, 2 < Z < 10 implies weak to moderate evidence, and Z > 10 implies strong evidence for module preservation.

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-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).

Figure 3

Global gene module characterization-M1A

A. Heatmap of genes in M1A. Samples are listed in columns and genes in rows. Samples bar under the hierarchal clustering tree were marked as red indicating cases and as blue indicating controls. Normalized expression values ranged between -2 and 2, as shown in the color legend under the heatmap. B. Highly connected genes in M1A network. C. Top five enriched gene ontology categories of highly corrected genes in the M1A module.

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].

Figure 4

Global gene module characterization-M3A

A. Heatmap of genes in M3A. Samples listed in column and genes in rows. Samples bar under the hierarchal clustering tree were marked as red indicating cases and as blue indicating controls. Normalized expression values ranged between -2 and 2, as shown in the color legend under the heatmap. B. Highly connected genes in M3A core network. C. Top five gene ontology categories enriched in the M3A module.

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.