Statistical epistasis and progressive brain change in schizophrenia: an approach for examining the relationships between multiple genes.

Although schizophrenia is generally considered to occur as a consequence of multiple genes that interact with one another, very few methods have been developed to model epistasis. Phenotype definition has also been a major challenge for research on the genetics of schizophrenia. In this report, we use novel statistical techniques to address the high dimensionality of genomic data, and we apply a refinement in phenotype definition by basing it on the occurrence of brain changes during the early course of the illness, as measured by repeated magnetic resonance scans (i.e., an 'intermediate phenotype.') The method combines a machine-learning algorithm, the ensemble method using stochastic gradient boosting, with traditional general linear model statistics. We began with 14 genes that are relevant to schizophrenia, based on association studies or their role in neurodevelopment, and then used statistical techniques to reduce them to five genes and 17 single nucleotide polymorphisms (SNPs) that had a significant statistical interaction: five for PDE4B, four for RELN, four for ERBB4, three for DISC1 and one for NRG1. Five of the SNPs involved in these interactions replicate previous research in that, these five SNPs have previously been identified as schizophrenia vulnerability markers or implicate cognitive processes relevant to schizophrenia. This ability to replicate previous work suggests that our method has potential for detecting a meaningful epistatic relationship among the genes that influence brain abnormalities in schizophrenia.

Introduction

Schizophrenia is conceptualized as a genetically complex disorder that is probably the result of multiple alleles of small effect in many different genes, at least some of which are likely to interact with one another.[1] Despite the widely-recognized plausibility of this conceptualization, progress in delineating the genetic and neural mechanisms of schizophrenia has been limited by the difficulty of developing adequate models that can account for genetic complexity, particularly epistasis. Epistasis has several connotations: one in which a gene suppresses the effect of another; and a number of genes interacting together to influence the expression of the phenotype. Here, we use epistasis in the latter sense.

Phenotype definition has also been a major challenge for research on the genetics of schizophrenia. Schizophrenia is characterized by a variable mixture of psychotic, disorganized, and negative symptoms; it also has a broad variety of clinical outcomes, ranging from prolonged periods of remission to a disabling treatment-refractory course that is sometimes accompanied by cognitive decline.[2] Most current genetic studies do not take the complexity of the phenotype into account. Genome-wide association studies (GWAS) have been disappointing because of poor replicability, a problem likely to be due in part to weak or imprecise phenotype definition.[3-15] GWAS only require that subjects meet specified diagnostic criteria, such as those put forth in DSM IV or ICD 10. This limited approach to phenotype definition runs the risk of grouping clinically (and probably genetically) diverse patients under a single rubric, thereby reducing the likelihood that replicable and accurate associations will be found. In this report we describe several approaches that we have implemented to address the problems of examining epistasis and improving phenotype definition.

An inherent problem in the effort to link improved phenotype definition with complex genetic concepts such as epistasis is the difficulty inherent in detecting meaningful signals within the noisy vastness of the human genome. Various approaches have been developed to reduce the high dimensionality of genomic data, such as the use of gene ontology or pathway analysis.[16] Another approach for addressing high dimensionality is the use of statistical techniques such as machine learning algorithms (MLA). The utility of these approaches has been demonstrated recently in studies that have used techniques such as recursive partitioning to model epistatic relationships.[17] We have developed a strategy for using similar algorithms to identify potentially meaningful relationships between schizophrenia vulnerability genes. Our approach employs a MLA, stochastic gradient boosting using regression trees,[18] and performs a series of steps to winnow through a large genomics data set and to identify genes and SNPs that may interact with one another to produce brain tissue loss in schizophrenia.

In this report we also address the problem of phenotype definition by basing it on the organ that produces the clinical picture of the illness: the brain. Numerous well-replicated studies have documented that patients with schizophrenia have measurable brain tissue reduction when compared with healthy volunteers.[19-33] The tissue loss affects multiple brain regions, particularly overall cerebral volume and frontal and temporal lobe gray matter (GM) and white matter (WM); tissue loss in these regions is accompanied by a corresponding increase in cerebrospinal fluid (CSF). These brain volume reductions are present at the onset of the illness, but they also continue to progress over time, with the greatest severity occurring during the first few years after onset.[33] These tissue changes are considered to be due to an aberration in a neurodevelopmental process for several reasons: the typical onset of illness during adolescence, when the brain is undergoing major developmental changes; and the absence of an neuropathologic signature suggesting a degenerative process such as neuronal loss or gliosis.

Because we have conducted a prospective longitudinal study of patients with schizophrenia ascertained at the time of onset ("first episode patients") and followed them with repeated Magnetic Resonance (MR) scans for a time period up to eighteen years, we have been able to develop a refinement in phenotype definition by basing it on the occurrence of brain changes during the early course of the illness, as measured by repeated MR scans.[33]

Therefore, in this report we describe a data analytic strategy for evaluating epistatic relationships between multiple genes and their alleles with a biologically-based phenotype of schizophrenia. This strategy is summarized in Figure 1 and consists of seven steps:

Figure 1

Steps Involved in Stochastic Gradient Boosting

Measure brain change over time in multiple brain regions, using structural Magnetic Resonance (sMR) volumetric measures in a sample of first episode schizophrenia patients

Identify candidate genes for analysis, based on their relevance to neurodevelopmental processes that regulate brain change over time or their association with schizophrenia, and genotype patients for these genes (See Table 1)

Table 1

Fourteen Genes Evaluated for Epistasis

Gene*		Location	Association withSchizophrenia(score)**	Biological Plausibility ofProtein Product	Number of SNPs
BDNF	B	Chr11:27.63-27.7 Mb	+ (0.34)	Member of the nerve growth factor family induced by cortical neurons & necessary for survival of striatal neurons in brain.	14
Brain Derived Neurotrophic Factor
NRG1	N	Chr8:31.52-32.4 Mb	++ (1.07)	Signaling molecule involved in neurodevelopment and synaptic plasticity. Ligand for ERBB receptors in the CNS.	62
Neuregulin 1
DISC1	D	Chr1:229.83-230.24 Mb	+++++ (3.20)	Regulates neurite outgrowth, cortical development, embryonic and adult neurogenesis.	130
Disrupted in Schizophrenia 1
RELN	R	Chr7:102.9-103.42 Mb	++ (1.10)	A cell matrix protein involved in cell-cell interactions important for cell positioning and neuronal migration during brain development.	165
Reelin
ERBB4	E	Chr2:211.96-213.11Mb	+++ (1.50)	A receptor tyrosine kinase activated by neuregulin which induces mitogenesis and neuronal differentiation	303
Avian Erythroblastic Leukemia Viral Oncogene Homolog 4
GDNF	G	Chr5:37.85-37.88Mb	+ (0.33)	Encodes for a neurotrophic factor involved in differentiation of dopaminergic neurons	16
Glial Cell Line-Derived Neurotrophic Factor
PDE4B	P	Chr1:66.03-66.61Mb	+++++ (3.00)	A phosphodiesterase that acts as a second messenger to regulate cellular responses to neurotransmitters	75
Phosphodiesterase 4B
DTNBP1	DT	Chr6:15.63-15.77Mb	++ (0.77)	Contributes to neuroplasticity via regulation of pre-synaptic vesicular transport; interacts with structural proteins and regulates receptor surface expression post-synaptically.	18
Dystrobrevin Binding Protein 1
COMT	C	Chr22:18.31-18.34Mb	+ (0.56)	Metabolizes catecholamines such as dopamine, epinephrine, and norepinephrine	23
Catechol-O-Methyltransferase
RGS4RG	RG	Chr1:161.31-161.31Mb	++ (1.07)	Brain expressed protein involved in negative regulation of G-protein activity, mRNA expression decreased in schizophrenia	7
Regulator of G-Protein Signaling 4
GRM3	GR	Chr7:86.11-86.33Mb	+ (0.46)	Excitatory L-glutamate neurotransmitter receptor.	34
Glutamate Receptor, Metabotropic, 3
DAOA	DA	Chr13:106.12-106.14Mb	++++ (2.25)	Regualtes NMDA receptor activity via multistep pathway that oxidizes D-serine	17
D-Amino Acid Oxidase Activator
KCNH2	K	Chr7:150.27-150.31Mb	++++ (2.00)	Brain expressed voltage-activated potassium channel, variants associated with differences in cognitive functioning	13
Potassium Channel, Subfamily H Member 2
AKT1	A	Chr14:104.31-104.33Mb	+++ (1.29)	Protein product is mediator of growth factor-induced neuronal survival during neurodevelopment	5
Protein Kinase B alpha

One or two letter abbreviations used in subsequent tables to denote the gene associated with a specific SNP

Association score generated from studies included in the SZGene database (http://www.schizophreniaforum.org/res/sczgene/default.asp). Each score is the ratio of the number of positive association studies divided by the number of negative association studies.

Relative + rankings determined by the association score:

if >2.4,

if >1.8 to 2.4,

if >1.2 to 1.8,

if >0.6 to 1.2, and

if ≤0.6.

Use stochastic gradient boosting (SGB) to reduce the number of genes to those that have the strongest relationship with multiple measures of brain change, as indicated by variable importance measures for genes and their associated SNPs (See Table 2 and Figure 2)

Table 2

Variable Importance* of Genes for Predicting Brain Volume Change Based on Initial Stochastic Gradient Boosting

BrainRegion	AKT1	BDNF	COMT	DISC1	DAOA	DTNBP1	ERBB4	GDNF	GRM3	KCNH2	NRG1	PDE4B	RELN	RGS4
Cerebral Tissue		4 - 1		34 - 7			19 - 4			5 - 1	10 - 2	12 - 2	10 - 2
Cerebral GM				19 - 3	4 - 1		34 - 8		4 - 1			9 - 2	17 - 4
Frontal GM			5 - 1	9 - 2			45 - 9			9 - 2	4 - 1	6 - 1	9 - 2
Temporal GM	5 - 1		5 - 1	18 - 4			28 - 6	6 - 1	4 - 1			9 - 2	10 - 2
Parietal GM			9 - 2	25 - 5			16 - 4		4 - 1		4 - 1		23 - 5
Cerebral WM			4 - 1	10 - 2			36 - 8		4 - 1		14 - 3	15 - 3	6 - 1
Frontal WM		6 - 1		13 - 3			23 - 5				6 - 1	15 - 3	20 - 4
Temporal WM		4 - 1		4 - 1			42 - 8				8 - 2	8 - 2	9 - 2	5 - 1
Parietal WM		4 - 1	5 - 1	10 - 2			39 - 8				5 - 1	6 - 1	18 - 4
VBR				29 - 6			17 - 3		7 - 1		15 - 3	4 - 1	13 - 3
Surface CSF			8 - 1	8 - 2	4 - 1		36 - 7				5 - 1	20 - 4	17 - 3
Caudate				18 - 4	4 - 1		38 - 7		5 - 1		15 - 3		12 - 3
Putamen				7 - 1		5 - 1	41 - 9			4 - 1	8 - 1	13 - 3	10 - 2
Thalamus		7 - 1	5 - 1	9 - 2	4 - 1		32 - 7				4 - 1	5 - 1	28 - 5

The variable importance (VI) measure for each gene, as a predictor of each brain volume measure, is calculated using the variable importance of the top 20 SNPs for predicting each brain volume measure. It is a percent of variable importance, calculated by adding up the VI scores for the top 20 SNPs and dividing each SNP VI score by the total of the 20. A genewise VI score is then calculated by adding up the Vis for each SNP to get a gene summary. The numbers in each row show the genewise VI, followed by the number of SNPs contributing to that VI. The VI adds up to 100% within each row (apart from rounding errors).

Figure 2

Variable Importance of Genes for Predicting Brain Volume Change Based on Initial Stochastic Gradient Boosting

Repeat SGB for the genes that survive the initial SGB to identify the top 20 SNPs with the highest VI; reduce data to those genes/SNPs that may have a significant epistatic relationship by determining which ones share an influence on the same specific brain regions

Reduce data for those genes/SNPs by using a GLM analysis to find the genes/SNPs that have a significant gene/SNP by gene/SNP interaction

Further reduce data by requiring that gene/SNP interactions affect at least two brain regions in a biologically plausible pattern (See Table 3)

Table 3

Epistatic Relationships Between Genes in Specific Brain Regions

Two WayInteractions	CerebralTissue		SurfaceGray		ParietalGray		CerebralWhite		FrontalWhite		TemporalWhite		Surface CSF		VBR		Caudate		Putamen
	F	Pr>F	F	Pr>F	F	Pr>F	F	Pr>F	F	Pr>F	F	Pr>F	F	Pr>F	F	Pr>F	F	Pr>F	F	Pr>F
PDE4B-RELN
P_rs12730529							9.27	0.0002	9.47	0.0001			8.19	0.0004
R_rs499953							3.46	0.03	3.31	0.04			6.22	0.003
P_rs12730*R_rs499953							7.94	<.0001	7.03	<.0001			4.45	0.002
P_rs544024	2.82	0.04											6.24	0.0005
R_rs2229860	0.25	0.62											2.14	0.15
P_rs54402*R_rs222986	5.35	0.02											12.90	0.0005
P_rs11576970							8.00	0.005	9.44	0.003
R_rs499953							4.19	0.02	4.57	0.01
P_rs11576*R_rs499953							4.48	0.01	5.05	0.008
ERBB4-PDE4B
E_rs7565622							4.84	0.009	3.41	0.04
P_rs12730529							5.23	0.006	4.48	0.01
E_rs75656*P_rs127305							5.46	0.0004	4.65	0.001
DISC1-RELN
D_rs3738401	6.87	0.01													7.18	0.008
R_rs580884	1.99	0.14													5.84	0.004
D_rs37384*R_rs580884	3.86	0.02													6.32	0.002
D_rs3738401	1.87	0.16													4.64	0.01
R_rs580884	0.60	0.55													2.97	0.05
D_rs37384*R_rs580884	4.61	0.004													3.05	0.03
DISC1-ERBB4
D_rs821589	9.57	0.0024					6.49	0.0120
E_rs6435632	8.49	0.0003					5.26	0.0063
D_rs82158*E_rs643563	5.44	0.0054					4.38	0.0143
D_rs821589							3.24	0.0739			5.98	0.0158
E_rs2118891							1.13	0.2890			4.02	0.0470
D_rs82158*E_rs211889							8.19	0.0049			14.21	0.0002
D_rs11578905													5.25	0.0064	9.17	0.0002
E_rs1357139													4.52	0.0126	11.27	<.0001
D_rs11578*E_rs135713													4.05	0.0040	6.99	<.0001
DISC1-PDE4B
D_rs11578905			4.84	0.0094	8.01	0.0005
P_rs2455012			0.99	0.3208	1.51	0.2213
D_rs11578*P_rs245501			3.53	0.0319	3.54	0.0319
NRG1-RELN
N_rs17664708																	4.05	0.0463	3.43	0.0660
R_rs2237628																	0.42	0.6578	1.39	0.2538
N_rs17664*R_rs223762																	3.12	0.0476	6.26	0.0025

Examine relationships between specific genotypes and severity of brain tissue change over time in order to identify potential epistatic relationships that affect disease progression in the brain

Materials and methods

Study sample

Subjects for this study were recruited through the Iowa Longitudinal Study (ILS). The ILS was initiated nineteen years ago and includes a total cohort of 542 first episode schizophrenia patients. These patients were recruited from consecutive admissions to the University of Iowa Psychiatry inpatient service at a rate of 25–30/year; recruitment ended in 2007. They have been followed at six month intervals after initial intake, with assessment of clinical symptoms, psychosocial function, and treatment received. Intensive assessments (sMR and cognitive testing) are done at intake and at two, five, nine, twelve, fifteen, and eighteen years.[33] In this report we focus on those subjects for whom we have both genomic data and adequate sMR data to provide a relatively definitive determination of whether progressive brain change occurs over a time interval that is up to 12 years after intake. These comprise a total of 144 patients for whom we have a minimum of 2 scans and a maximum of six. (See Table 1 in Supplementary Materials for demographics and clinical characteristics.)

sMR data acquisition and analysis

Acquisition and analysis of sMR data have been previously described.[33,39-44] Details are available in supplementary materials.

Genotyping

To genotype study participants, DNA was prepared by high-salt extraction from whole blood and assayed using Illumina Infininum II array BeadChips which were designed, manufactured and completed by Illumina (San Diego, CA). The locus success rate and genotype call rate on the 26,122-SNP BeadChips were 97.2% and 99.8% respectively. We used a customized Illumina chip provided by Ortho-McNeill Janssen that contains SNPs for 1204 genes; the SNPs are primarily tag SNPs. The genes were selected for inclusion on the chip based on their association with major mental illnesses, their role in neurotransmission, their role in regulating major metabolic or developmental pathways, and their role in drug metabolism. For this study we selected fourteen genes for examination of statistical epistasis in relation to brain change measures. Because we were interested in the specific examination of the neurodevelopmental mechanisms influencing brain change in schizophrenia, we selected the genes based on the following criteria: established association with schizophrenia; influence on neurodevelopmental pathways; neurotransmitters implicated in schizophrenia. The rationale for the selection of these fourteen genes (810 SNPs) is summarized in Table 1.

Statistical analysis

To create a measure of percent change for use in our analyses, the MRI volumes were divided by the intracranial volume to correct for variability in head size. We have previously found that the greatest brain changes occur early in the course of the illness.[33] Therefore, our measures are a percent change score based on the earliest available scan and its subsequent one. The value for the first scan was subtracted from the value for the second scan; that value was then divided by the number of years between the scans; the result was then divided by the value for the first scan, and finally that number was multiplied by 100 to provide an index of percent change.

The goal of our statistical approach was to identify a logical and empirically-driven method for detecting the signal (potential epistatic relationships between genes that may affect brain volume changes) within our vast array of genomic data. (See Figure 1)

Our first step was designed to determine which of the fourteen candidate genes had the strongest relationship with brain tissue change over time. We used regression trees to identify the relative importance of genes (SNPs) for predicting the amount of change in the size of each of the brain regions. This analysis employed an ensemble approach with stochastic gradient boosting (SGB) (TreeNet).[18] Variable importance (VI) scores, generated by the SGB analysis, were used to identify the SNPs and genes that were retained during the initial steps of variable reduction. Variable importance scores show the relative contribution of each of the variables (SNPs) to predicting the outcome (a measure of brain volume change).

VI scores are relative; the largest value is assigned a value of 100 and the remaining variables are scaled accordingly. The variable importance measure for each gene, as a predictor of each brain volume measure, is calculated using the variable importance of the top 20 SNPs for predicting each brain volume measure. It is a percent of variable importance, calculated by summing the VI scores for the top 20 SNPs and dividing each SNP VI score by the total of the 20. A genewise VI score is then calculated by adding up the VIs for each SNP to get a gene summary.

The steps involved are summarized as follows.

Begin with the entire data set and sample rows and columns of the data matrix (SNPs and brain measures for each individual); randomly sample a subset of individuals/SNPs

Build a small tree (maximum of 8 nodes) and compute variable importance for SNPs

Randomly resample individuals and SNPS and repeat, building an ensemble of 1000 trees

Aggregate VI across the ensemble

From SNPs identified in the second SGB, we used general linear models (GLM) to further assess SNP*SNP interaction effects on brain volume changes over time. Analyses were performed using SAS (version 9.2; SAS, Cary, North Carolina). Percent brain volume change measure was entered as the dependent variable. Genotype of SNP pairs and SNP*SNP interaction terms were entered as independent measures in the statistical models. For all analyses, intracranial volume at initial MR scan, gender, imaging protocol (MR5 versus MR6), age at initial MR scan, and amount of neuroleptic exposure measured using dose-years were included as covariates. Because our technique is exploratory, we did not correct for multiple comparisons; we report only findings with P<.02.

Results

Table 2 and Figure 2 show the VI of the fourteen genes for predicting severity of brain volume changes, broken down into fourteen specific brain regions. Three genes stand out by virtue of having a high variable importance for all the brain regions that we examined: DISC1, ERBB4, and RELN. Two others, NRG1 and PDE4B, had high variable importance for twelve regions. The remaining genes had lesser (e.g., BDNF, COMT, and GRM3) or very little importance (e.g., AKT1, DAOA, GDNF, DTNBP1, KCNH2, and RGS4). We therefore retained the top five genes for further examination.

We then repeated SGB using these five genes and their associated SNPs (N=735). We again identified the top 20 SNPs for each brain region based on VI. The analysis further reduced the “candidate SNPs” to a total of 267 (Table 2, Supplementary Materials). A prominent role continued to be played by the SNPs belonging to DISC1 (58 discrete SNPs), ERBB4 (87 discrete SNPs), and RELN (61 discrete SNPs). NRG1 and PDE4B contribute fewer (23 and 38 respectively). We then reduced these data by identifying groups of SNPs that may have an epistatic relationship because they share an association with a given specific brain region. These potential interactions were then tested for statistical significance by being placed in a general linear model, using relevant covariates, and determining whether there was a significant SNP*SNP interaction. Eighteen SNP pairs were found to pass this filter.

The final filter in the search consisted of examining these pairs in relation to the brain regions with which they were associated. We limited this final group to those SNP pairs that had a significant relationship with more than one brain region and that was biologically plausible (e.g., reduced cerebral volume accompanied by increased CSF, as an indicator of a generalized tissue loss; reduced WM in two regions as indicator of tissue specificity). Eleven interactions were selected for analysis using these criteria. They are shown in Table 3. They include all five genes retained up to this point and a total of 17 SNPs: 5 for PDE4B, 4 for RELN, 4 for ERBB4, 3 for DISC, and one for NRG1.

Five of these SNPs have been previously identified as increasing schizophrenia risk or are in linkage disequilibrium with known schizophrenia risk alleles or alleles affecting a cognitive function relevant to schizophrenia.[45-52] Many of them involve RELN. Two of the SNPs involve an interaction between RELN and PDE4B alleles.

One is an interaction between RELN rs2229860 and PDE4B rs54402. This interaction is associated with a decrease in cerebral tissue and an increase in surface CSF. Rs2229860 is an identified schizophrenia risk allele on the RELN gene; it is a rare A/G mutation that results in a proline to arginine change at position 1703.[46] A second interaction is between RELN rs580884 and DISC1 rs3738401. This interaction is associated with decreases in cerebral tissue and an increase in VBR. Rs3738401 is an A/G polymorphism that is a missense mutation; the presence of the minor A allele leads to an arginine to glycine change at position 844. The significant epistatic relationship appears to derive from an AA/TG combination, which occurred in 8 subjects. When these two genotypes co-occur, they have the largest amount of cerebral tissue loss (−.769%) and the largest amount of VBR increase (14.65%). Both are statistically significant (p<.01). A third interaction is between RELN rs499953 and PDE4B rs11576970. This interaction is associated with decreases in frontal and cerebral white matter. The PDE4B SNP is located in or near block 5, 2.8 kilobases away from rs7412571, with which it is in linkage disequilibrium; rs7412571 is significantly associated with schizophrenia.[53] The epistatic relationship derives from a TC/GG combination, resulting in a significant decrease in cerebral (−6.10%, p<.01) and frontal (−6.62%, p<.03) WM.

A fourth interaction is between DISC1 rs11578905 and PDE4B rs2455012. This interaction is associated with a decrease in surface and parietal GM. The PDE4B SNP is in LD with two identified alleles, rs1354064 and rs2503222, both of which increase risk for schizophrenia.[48,53] A final interaction involving a previously identified allele is between NRG1 and RELN. This interaction is associated with an increase in the volume of the caudate and putamen. The RELN SNP, rs2237628, is in LD with rs2711870, which has been implicated with problems in shifting response set on a card sorting task, an impairment in executive function that has repeatedly been shown to be abnormal in schizophrenia.[54] The interaction of the RELN and the NRG1 heterozygotes leads to an increase of 2.22% in the caudate and 3.90% in the putamen.

Discussion

We have described a method for identifying potential epistatic relationships between genes that may contribute to schizophrenia and affect disease progression in the brain in patients suffering from schizophrenia. The method combines a machine learning algorithm, the ensemble method using stochastic gradient boosting, with traditional general linear model statistics. The method was used to identify genes/SNPs that were interacting with one another and predicting a continuous outcome measure that is a biologically meaningful phenotype (“intermediate phenotype”) for schizophrenia: changes in brain structure occurring after the onset of the illness. The method began with fourteen genes and 810 SNPs and reduced them to five genes and seventeen SNPs, identifying eleven interactions. Five of the SNPs involved in these interactions replicate previous research, in that these five SNPs have previously been identified as schizophrenia vulnerability markers or implicate cognitive processes relevant to schizophrenia. This ability to replicate previous work suggests that our method has potential for detecting a meaningful signal within the human genome.

RELN was the gene that was found to have the highest number of interactions. It had a total of six, involving PDE4B, DISC1, and NRG1. RELN’s relationship with schizophrenia is well-supported by linkage or association studies,[45,46] and it has also been shown to be associated with treatment resistance.[57] RELN expression is markedly reduced (by 50%) in mRNA and protein in post mortem brain tissue from schizophrenia patients.[58] Epigenetic abnormalities have also been found; hypermethylation in the promotor region may be an explanation for the reduced expression of RELN in the brain.[55,56] The known functions of RELN are highly consistent with current views that schizophrenia may be due to a disruption in neurodevelopmental processes. RELN encodes a large glycoprotein that affects both early and late neurodevelopment. It is responsible for neuronal migration, cell-cell interactions, and positioning of proliferating neurons early in development, and it also modulates neuronal and synaptic plasticity throughout life.[57]

DISC1 was also found to have significant epistatic relationships with other genes, including PDE4B and ERBB4 in addition to RELN. DISC1 is a well-established schizophrenia vulnerability gene.[59-62] It has multiple SNPs that have been extensively studied, including one identified in the present study (rs3738401).[63-69] The prevailing model for the role of DISC1 in schizophrenia emphasizes haploinsufficiency. DISC1 has decreased expression in mRNA and protein in both post mortem tissue and lymphoblastoid cell lines from patients with schizophrenia.[66] DISC1 functions are consistent with both the neurodevelopmental hypothesis and current thinking about the mechanisms of schizophrenia. It is a hub protein whose interactors modulate expression of neurodevelopmental, synaptogenic, and sensory perception genes.[61-62] Its established interactome includes PDE4B, PDE4D, NDE1, and NDEL1; however, pathway analysis suggests that it may in fact be much larger.[70] The DISC1 SNP identified in this analysis, rs3738401, is among those identified as a psychoactive drug target for the treatment of psychosis and/or schizophrenia through an Ingenuity pathway analysis.[70]

PDE4B was identified as having an epistatic relationship with DISC1 in our analyses, as well as with RELN and ERBB4. PDE4B has also been identified as a risk factor for schizophrenia.[48-49,53] It is a large gene (580 kb) with 17 exons; at least five isoforms occur as a result of alternative splicing. PDE4B2 and PDE4B4 isoforms have been found to be reduced in post mortem tissue from schizophrenia patients.[48] It is orthologous to Drosophila dunce, a gene involved in learning and memory; dunce mutants have altered axonal grown cone motility, synaptic plasticity, and neuronal function.[47] It also influences myelination, consistent with the findings in this study of significant WM loss associated with PDE4B genotypes. Its interaction with DISC1 occurs when the N-terminus of DISC1 binds to the UCR2 regulatory domain of PDE4B.[71] This interactome plays a critical role in CNS signaling and homeostasis. DISC1 sequesters PDE4B in an inactive state and releases it when cyclicAMP is elevated; it thereby regulates intracellular signaling by inactivating cAMP.

NRG1 was found to have an epistatic relationship with RELN in our study. It has been implicated in multiple studies as a potential susceptibility gene for schizophrenia.[72-76] NRG1 is a large gene (1200 kb) that encodes at least 15 isoforms through alternative splicing.[75] It is responsible for a host of functions that affect brain development and cell-cell communications and that are relevant to schizophrenia; they include axon guidance, myelination, glial differentiation, synaptogenesis, synaptic plasticity, and neurotransmission.[75] mRNA expression studies implicate variation in isoform expression, due to alternative splicing, as a the most likely molecular mechanism for the genetic association of NRG1 with schizophrenia.[77]

ERBB4 would be predicted to have an epistatic relationship with NRG1, based on the fact that ERBB4 encodes a receptor for NRG1,[75,78] but none was found in our analyses. It did have relationships with PDE4B and DISC 1, which are consistent with the hub protein nature of DISC1. ERBB4 is part of a family of genes for protein tyrosine kinases.[75] We found an interaction with a known RELN SNP that is related to cognitive function. This relationship is limited exclusively to two basal ganglia structures, the caudate and putamen. Unlike other components of brain parenchyma, these two structures have increased in size. The epistatic relationship may provide some clue as to the mechanisms of antipsychotic action and, perhaps, toxicity.

This study suggests that exploratory analyses using MLAs may be a fruitful approach for examining the relationships between multiple genes that may be implicated in schizophrenia. It is noteworthy that five of the SNPs identified through this approach have been previously identified as risk alleles, in LD with risk alleles, or relevant to functions that are disrupted in schizophrenia. This replication of prior findings is of particular interest because the prior studies have generally been based on case-control comparisons, while the current study is a case only study, and brain measures rather than diagnosis were the outcome variable. This study also suggests that examination of genomic factors in schizophrenia may be enhanced by using a novel phenotypic outcome measure: brain change over time as measured by MR imaging.

This study has several limitations. It is a case-only sample; an examination of the genetic associations with brain change in a case-control setting would strengthen these results. It is based on a relatively small sample and is primarily useful for hypothesis generation. Because it is an exploratory study, we did not correct for multiple comparisons. The findings need replication in an independent sample. Finally, from a genetic perspective, we know that copy number variations (CNVs) and epigenetic factors such as methylation play an important role in conferring liability to disease. We have not examined any of these factors in this study.