Suicide candidate genes associated with bipolar disorder and schizophrenia: An exploratory gene expression profiling analysis of post-mortem prefrontal cortex

Suicide candidate genes in bipolar disorder and schizophrenia

Between the suicide vs. non-suicide group within schizophrenia, 70 genes were differentially expressed (Table 3). Most of these genes were down-regulated. From the above lists of differentially expressed genes, within diagnostic groups, two genes overlapped (Fig. 1A). Specifically, the phospholipid scramblase 4 (PLSCR4) and empty spiracles homolog 2, Drosophila (EMX2) genes were down-regulated in both suicide groups compared to the non-suicide groups. As negative controls, the normalization control probe set of 100 genes were tested by the same ANOVA model, and no genes met our statistical criteria (Fold Change ≥ |1.3| and FDR < 0.1) between the suicide group vs. non-suicide group within bipolar or schizophrenia diagnostic categories.

Gene symbol	Description	FC1	FDR2	P/A call3
GJA1	gap junction protein, alpha 1, 43 kDa (connexin 43)	-2.0	0.02	P
MT1X	metallothionein 1X	-1.9	0.04	P
MT1M	metallothionein 1M	-1.9	0.06	P
AGXT2L1	alanine-glyoxylate aminotransferase 2-like 1	-1.9	0.03	P
SDC4	syndecan 4 (amphiglycan, ryudocan)	-1.7	0.01	P
SOX9	SRY (sex determining region Y)-box 9	-1.7	0.02	P
NTRK2	neurotrophic tyrosine kinase, receptor, type 2	-1.7	0.02	P
HSPB1	heat shock 27 kDa protein 1	-1.7	0.04	P
EIF5A	eukaryotic translation initiation factor 5A	-1.7	0.06	P
EFEMP1	EGF-containing fibulin-like extracellular matrix protein 1	-1.6	0.04	P
ID4	inhibitor of DNA binding 4, dominant negative helix-loop-helix protein	-1.6	0.02	P
APOE	apolipoprotein E	-1.6	0.03	P
PLSCR4	phospholipid scramblase 4	-1.6	0.03	P
AGT	angiotensinogen (serpin peptidase inhibitor, clade A, member 8)	-1.6	0.02	P
TUBB2B	tubulin, beta 2B	-1.5	0.04	P
MT1E	metallothionein 1E (functional)	-1.5	0.03	P
FAM107A	family with sequence similarity 107, member A	-1.5	0.03	P
SLC4A4	solute carrier family 4, sodium bicarbonate cotransporter, member 4	-1.5	0.05	P
GLUL	glutamate-ammonia ligase (glutamine synthetase)	-1.5	0.04	P
ALDH1L1	aldehyde dehydrogenase 1 family, member L1	-1.5	0.02	P
DDIT4	DNA-damage-inducible transcript 4	-1.5	0.04	P
MT1G	metallothionein 1G	-1.5	0.04	P
SLC1A3	solute carrier family 1 (glial high affinity glutamate transporter), member 3	-1.5	0.08	P
PPAP2B	phosphatidic acid phosphatase type 2B	-1.5	0.05	P
EMX2	empty spiracles homolog 2 (Drosophila)	-1.5	0.05	P
TP53BP2	tumor protein p53 binding protein, 2	-1.5	0.03	P
MT1H	metallothionein 1H	-1.5	0.02	P
APOLD1	apolipoprotein L domain containing 1	-1.4	0.08	P
FGFR3	fibroblast growth factor receptor 3 (achondroplasia, thanatophoric dwarfism)	-1.4	0.03	P
FXYD1	FXYD domain containing ion transport regulator 1 (phospholemman)	-1.4	0.00	P
NOTCH2	Notch homolog 2 (Drosophila)	-1.4	0.03	P
METTL7A	methyltransferase like 7A	-1.4	0.05	P
SLC14A1	solute carrier family 14 (urea transporter), member 1 (Kidd blood group)	-1.4	0.06	P
BBOX1	butyrobetaine (gamma), 2-oxoglutarate dioxygenase 1	-1.4	0.03	P
MT1F	metallothionein 1F (functional)	-1.4	0.05	P
EPHX1	epoxide hydrolase 1, microsomal (xenobiotic)	-1.4	0.02	P
SLC7A11	solute carrier family 7, (cationic amino acid transporter, y+ system) member 11	-1.4	0.05	P
TPD52L1	tumor protein D52-like 1	-1.4	0.02	P
RHOBTB3	Rho-related BTB domain containing 3	-1.4	0.05	P
VIL2	villin 2 (ezrin)	-1.4	0.03	P
NDRG2	NDRG family member 2	-1.4	0.06	P
CLDN10	claudin 10	-1.4	0.06	P
NTSR2	neurotensin receptor 2	-1.4	0.03	A
ITPKB	inositol 1,4,5-trisphosphate 3-kinase B	-1.4	0.10	P
SDC2	syndecan 2 (heparan sulfate proteoglycan 1, cell surface-associated, fibroglycan)	-1.4	0.03	P
IL17RB	interleukin 17 receptor B	-1.3	0.03	P
SPON1	spondin 1, extracellular matrix protein	-1.3	0.02	P
LOC645745	metallothionein 1H-like protein	-1.3	0.03	P
ITGB4	integrin, beta 4	-1.3	0.03	P
MLC1	megalencephalic leukoencephalopathy with subcortical cysts 1	-1.3	0.05	P
ATP1A2	ATPase, Na+/K+ transporting, alpha 2 (+) polypeptide	-1.3	0.05	P
RAB31	RAB31, member RAS oncogene family	-1.3	0.05	P
ATP1B2	ATPase, Na+/K+ transporting, beta 2 polypeptide	-1.3	0.05	P
AHCYL1	S-adenosylhomocysteine hydrolase-like 1	-1.3	0.03	P
CX3CR1	chemokine (C-X3-C motif) receptor 1	2.1	0.02	P
C3	complement component 3	1.6	0.02	P
CRH	corticotropin releasing hormone	1.6	0.05	P
ABCG2	ATP-binding cassette, sub-family G (WHITE), member 2	1.5	0.03	P
DUSP6	dual specificity phosphatase 6	1.5	0.04	P
NPY	neuropeptide Y	1.4	0.03	P
SLC25A23	solute carrier family 25 (mitochondrial carrier; phosphate carrier), member 23	1.4	0.03	A
LAPTM5	lysosomal associated multispanning membrane protein 5	1.4	0.03	P
TNFSF10	tumor necrosis factor (ligand) superfamily, member 10	1.4	0.01	P
TYROBP	TYRO protein tyrosine kinase binding protein	1.4	0.03	P
CD74	CD74 molecule, major histocompatibility complex, class II invariant chain	1.4	0.03	P
P2RY13	purinergic receptor P2Y, G-protein coupled, 13	1.3	0.02	A
CSF1R	colony stimulating factor 1 receptor	1.3	0.03	P
HLA-DRA	major histocompatibility complex, class II, DR alpha	1.3	0.04	A
HLA-DPA1	major histocompatibility complex, class II, DP alpha 1	1.3	0.05	P
A2M	alpha-2-macroglobulin	1.3	0.03	P

¹FC represents ratio of geometric means between suicides vs. non-suicide.

² P-values were adjusted by FDR.

³Absent (A) or Present (P) detection call by MAS5 statistical algorithm

Real time-PCR (RT-PCR) tested the validity of these two shared, differentially expressed genes in the available subset of the microarrayed samples. In the schizophrenia suicide group, the EMX2 and PLSCR4 expression levels were significantly down-regulated by comparison of mean expression levels (EMX2 t(9) = 2.42, p = 0.02; PLSCR4 t(18) = 3.77, p = 0.0005) in the suicide group compared to the non-suicide group (Fig 1B). The estimated fold changes in the suicide group were -1.51 for EMX2 and -2.16 for PLSCR4 relative to the non-suicide group. These differences were consistent with our microarray data. However, in bipolar disorder, these two genes could not be validated in a subset of tissues available from the original microarray study patients (Fig. 1B). These results highlight the overall findings that few common differentially genes for suicide vs. non-suicide exist between diagnostic groups.

Biological process in the differentially expressed genes

Functional annotation of the differentially expressed genes by Gene Ontology indicated that 9 biological processes were significantly overrepresented (at level 4; P < 0.05) among the suicide candidate genes in the schizophrenia cohort (Fig. 2). The transport genes were the most commonly over-represented biological process in our suicide candidate gene list for schizophrenia, including the glial, high affinity, glutamate transporter (SLC1A3). In contrast, no significantly over-represented biological process group emerged with suicide candidate genes in the bipolar disorder group.

Additionally, we analyzed the list of differentially expressed genes for each diagnosis by the Ingenuity Pathways Analysis (IPA) software to identify biological pathways and networks. We identified distinct signaling networks from suicide candidate genes that included the EMX2 gene in both disorders (Fig. 3). In bipolar disorder, the pathway perspective suggested a signaling network related to both cellular movement and cell to cell signaling, with interactions encompassing 10 differentially expressed, suicide candidate genes (Fig. 3A). By contrast, in schizophrenia patients, the differentially expressed genes were related in a cell death signaling network (Fig. 3B).

Microarray data and Patient Samples

The brain tissues were meticulously collected in a standardized manner via pathologists in the offices of the Medical Examiner in several states with the families' permission under the aegis of the Stanley Foundation Brain Collection (Array Collection plus Consortium Collection) [26]. The selection of specimens, clinical information, diagnoses of patients, and processing of tissues were conducted by Stanley Foundation Consortium as described previously [26]. Gene expression profiling utilized post-mortem prefrontal cortex (Brodmann's Area 46/10) mRNA and Affymetrix Human Genome U133 Set A (HGU133A) using standardized techniques as described [26, 27]. The prefrontal cortex was selected as the region of interest due to its role in executive functioning, impulsivity ("lack of premeditation"), and decision making. Disadvantageous decision making and impulsivity have been found to increase the risk of suicide [28, 29].

The Stanley Foundation's microarray database is an anonymous, de-identified dataset without any protected health information. Patients' demographic variables used in this study are listed in Table 1.

The robust multi-array averages (RMA)-normalized microarray data from four independent studies were downloaded from the SMRIDB. Microarray data from the same platform, Affymetrix Human Genome U133 Set A (HGU133A), were used to avoid platform-to-platform variation. The platform contains 22,215 probe sets. Quality control analyses for each chip were described previously [15]

For the bipolar disorder cohort, the total dataset consisted of 49 suicide completers' gene chips and 58 non-suicide gene chips, while for schizophrenia cohort, the total dataset consisted of 22 suicide completer gene chips and 89 non-suicide chips. Among 45 bipolar samples, there were 22 suicide cases, and 23 non-suicide cases. Among 45 schizophrenia patients, there were 10 suicide cases, and 35 non-suicide cases. Two to three microarray chip datasets were generated from the each patient's sample. These repeated microarray data from each patient were treated as technical replicates.

Statistical analysis of microarray data

Microarray data was analyzed by a statistical method described previously with slight modifications [26]. Briefly the following steps were followed within each diagnostic group (see Table 1). First, the differentially expressed genes between suicide completers vs. non-suicide groups were filtered by average fold change (FC ≥ |1.3|) using the BRB-array tool[30] without covariates. Second, the influence of continuous demographic variables (such as age, post-mortem interval (PMI) and brain pH) with the nominal variable suicide was tested using ANOVA. Then, categorical variables such as sex, smoking, alcohol and drug abuse were tested using chi square tests of association (Statview software SAS, Cary, NC). In addition, correlation analyses of the demographic factors with expression levels of the differentially expressed probe sets from step 1 were performed. Continuous variables were analyzed by Spearman's rank correlation and categorical variables were tested by ANOVA. P-values were adjusted by False Discovery Rate (FDR) in both tests [31]. Third, significant confounding factors were tested as possible covariates for ANCOVA model inclusion with the following criteria: The variable was required to show both 1) significant association with suicide as well as 2) significant correlation with expression levels of the differentially expressed genes. However, no variables met the criteria in both disorder groups. Therefore, no covariates were used in the omnibus ANOVA, using the factor as suicide vs. non-suicide. As an exploratory analysis, a more liberal FDR P-value (< 0.1) of significance was selected for expression level differences as previously described [32, 33], using the BRB array software tool's FDR default setting. For negative controls, we performed statistical analysis with the HGU133A normalization control probe sets using the same ANCOVA models, ensuring the adjustment did not produce "noise" or aberrant false positives. The Microarray Suite, version 5.0 (MAS5) software was used to filter genes with low expression levels as either present or absent, applying the detection call statistical algorithm. This algorithm suggests whether a gene is present or absent.

A power analysis estimated the sample sizes for detection of a 1.3 fold change in a gene with a significance criterion of P-value = 0.001 and a power of 0.90 using a previously described method [13]. This analysis estimated a minimum sample size of 27 cases per group for comparing suicide completers vs. non-suicide groups within bipolar disorder and 21 samples per group for comparing suicide completers vs. non-suicide completers within schizophrenia.

Real-time quantitative PCR

Total RNA from the dorsolateral prefrontal cortex (Brodmann area 46) of the Array Collection was used for this experiment. Complementary DNA was synthesized from DNA-free RNA with a random hexamer primer and Superscript III First-Strand Synthesis System according to the manufacturer's protocol (Invitrogen). Using a 384-well format with the Prism7900HT real-time detector (ABI), 2 μl aliquots of (10×) QuantiTect Primer Assay (validated primers to the specific gene of interest; Qiagen), 10 μl (2×) QuantiTect SYBR PCR Master mix (Qiagen), and 8 μl cDNA were mixed together for 20 μl total reaction volume and pipetted into single wells of the 384 PCR plate. Amplification conditions were: (1) 1 cycle for 2 min at 50°C, (2) 1 cycle for 15 min at 95°C, and (3) 45 cycles for 15 s at 94°C, 30 s at 60°C and 30 s at 72°C and fluorescence was measured during the 72°C step for each cycle as recommended by the manufacturer. The β-2 microglobulin (B2M) was chosen as an endogenous control for the normalization of target genes as it was consistently expressed in microarray samples. A total of 32 samples per each disorder were used for this experiment and run in duplicate. In the bipolar disorder cohort, there were 14 suicides and 18 non-suicide cases. In schizophrenia cohort, there were 5 suicides and 27 non-suicide cases. These samples were matched by age, race, gender, PMI, brain pH, side of the brain and quality of RNA. Reactions were quantified by the comparative Ct method using SDS2.2 software (ABI). This RT-PCR data was also statistically analyzed by the amplification plot method using the Data Analysis for Real Time PCR (DART-PCR) approach [34]. This method identified outliers in amplification efficiency by ANOVA and calculated mean expression levels. Statistical differences in expression levels between groups, namely suicide completers vs. non-suicides within a diagnostic group, were tested by one-tail, t-test with unequal variances (Microsoft Excel) as described [35, 36]. To estimate average fold changes between groups, the mean expression values from the DART-PCR approach were used. The alternative 2(-Delta Delta Ct) method [37] for estimating fold change (using all the data without exclusions) verified this fold change estimate. As both methods gave similar estimates, only the DART-PCR approach estimates from mean expression levels were reported.

Functional annotation

The differentially expressed genes were functionally annotated using the DAVID integrated database query tool[38] and by the over-representational analysis method [39]. Functional annotations were based on biological process of Gene Ontology (GO) Consortium[40] at level 4. P-values less than 0.05 were considered significant.

Pathway Analysis

Biologically relevant networks were drawn from the lists of genes that were differentially expressed in bipolar disorder and schizophrenia. This data was generated through the use of Ingenuity Pathways Analysis (IPA) [41], a web-delivered application that enables the visualization and analysis of biologically relevant networks to discover, visualize, and explore relevant networks. Expression data sets containing gene identifiers (Affymetrix probe set ID) and their corresponding expression values as fold changes were uploaded as a tab-delimited text file. Each gene identifier was mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base. These genes, called Focus Genes, were then used as the starting point for generating biological networks. To start building networks, the application program queries the Ingenuity Pathways Knowledge Base for interactions between Focus Genes and all other gene objects stored in the knowledge base, and generates a set of networks. The program then computes a score for each network according to the fit of the network to the set of focus genes. The score indicates the likelihood of the Focus Genes in a given network being found together due to random chance. A score of greater than 2 indicates that there is a less than 1 in 100 chance that the Focus Genes were assembled randomly into a network due to random chance. The scores of the networks generated from the lists of differentially expressed genes were 24 for Bipolar disorder (Fig 3A) and 40 for Schizophrenia (Fig 3B).