Axon guidance pathway genes are associated with schizophrenia risk.

In the present study, we analyzed schizophrenia (SCZ)-related genome-wide association studies (GWAS) to identify genes and pathways associated with SCZ. We identified 1,098 common genes (1,098/9,468) and 20 shared KEGG pathways (both P<0.01) by integrating candidate genes from the European and American SCZ-related GWAS. The pathways related to axon guidance, long term potentiation and arrhythmogenic right ventricular cardiomyopathy (ARVC) were highly significant (P<10-3). Moreover, 15 axon guidance pathway-related genes were associated with SCZ. The association between axon guidance pathway genes and SCZ was validated by a two-stage case-control study on Shandong migrants in northeastern China. Moreover, individuals with the rs9944880 TT polymorphism in the deleted in colorectal cancer (DCC) gene were associated with SCZ. These findings indicate that the axon guidance pathway genes and the rs9944880 SNP in DCC gene are associated with SCZ pathogenesis.

Introduction

Nearly 1% of the world population is affected by schizophrenia (SCZ), which is a severe neuropsychiatric disorder that is characterized by cognitive dysfunction (1–3). The etiology of SCZ pathogenesis is complex and lacks effective therapeutic program. The role of genetics in SCZ has been widely explored. Studies on twins and families have indicated that genetic factors play a key role in the onset of SCZ. The concordance rate of monozygotic twins is ~50%, which is higher than fraternal twins and other siblings (4). In a family without SCZ patients, the morbidity of first relatives is lower than the lineal consanguinity of SCZ patients (5). Estimated heritability of SCZ ranges from 65 to 80% (6,7). Recently, genome-wide association studies (GWAS) has been used to identify molecular and genetic mechanisms underlying SCZ (8) and several hundred loci that are related to SCZ risk have been identified (9–20). However, analyses of GWAS results are not robust and account for only a small fraction of critical genes.

Although many genes are associated with SCZ susceptibility, the mechanistic details are yet to be elucidated. The in silico methods are the best tools available to functionally annotate genetic factors that contribute to SCZ. Among these, the pathway-based method is the most reliable way to functionally analyze GWAS data. This method assumes that genes belonging to the same pathway are most likely to confer the risk for a complex disease and is used to identify the enriched functional categories (21,22). For example, Wang et al used the pathway-based approach to identify 102 genes in sterol transport and metabolism that were associated with high-density lipoprotein cholesterol levels (23). The pathway-based analysis provides biological insights that are not readily available from a single marker analysis and is more effective in replicating different independent studies. Therefore, in this study, we performed Gene Set Enrichment Analysis (GSEA) of two independent SCZ-related GWAS to identify the pathways and genetic factors that contribute to SCZ. We further validated our findings in a two-step case-control study involving Shandong migrants from northeastern China.

Materials and methods

Ethics statement

All procedures in this study were conducted as approved by the Ethics Committee of Qiqihar Medical University (Qiqihar, China) and in accordance with the Helsinki Declaration. All participants or their guardians gave written informed consent prior to inclusion.

Study population

The northeastern part of China is a border area with frequent population migration, especially from the Shandong province. In this study, we included two independent cohorts of migrants that self-reported as originating from Shandong. The initial discovery cohort consisted of 600 unrelated Han SCZ patients (286 females and 314 males; mean age, 37.3 years) and 600 control subjects (278 females and 322 males; mean age, 38.2 years). For validation, we included an independent sample set consisting of 1,731 SCZ patients (829 females and 902 males; mean age, 34.6 years) and 1,849 controls (936 females and 913 males; mean age, 35.0 years). The patients were recruited from the Second Affiliated Hospital of Qiqihar Medical University, Harbin First Specific Hospital, Mudanjiang Mental Hospital, Wealth Mental Hospital and Beian Mental Hospital in Heilongjiang province. The patients were clinically diagnosed according to the Chinese Classification of Mental Disorders version 3 (CCMD-3) and the Diagnostic and Statistical Manual of Mental Disorders-IV (DSM-IV). All control individuals were enrolled from the Third Affiliated Hospital of Qiqihar Medical University and matched for age and sex to patients. Two experienced psychiatrists interviewed the healthy controls individually using the Structured Clinical Interview for DSM-IV-TR Axis I Disorders-Nonpatient Edition (SCID-NP), and controls were clinically determined to be free of psychiatric disorders and individual or family history of mental illness.

Patient sample preparation and genotyping

Genomic DNA was isolated from peripheral white blood cells using a commercial Axygen Blood Genomic DNA Miniprep kit (Axygen Biosciences, Hangzhou, China). The DNA samples were verified by 0.8% agarose gel electrophoresis and only those preparations with no signs of degradation were used to construct DNA pools. The DNA concentrations were quantified in a NanoDrop™ 2000 spectrophotometer (Thermo Fisher Scientific, Inc., Waltham, MA, USA) and each sample was diluted to 50 ng/µl. The DNA pooling was performed thrice.

For allelotyping, six DNA pools from the discovery sample set were individually genotyped on Affymetrix GeneChip Human Mapping 6.0 chips (Affymetrix: Thermo Fisher Scientific, Inc., Santa Clara, CA, USA). The processing of the genomic DNA and the array hybridization were performed according to the manufacturer's instructions. The hybridization intensities of two probes for each allele were derived from the raw scanning files.

For the validation of the candidate loci, DNA samples from the validation panel consisting of 1,731 SCZ patients and 1,849 controls were genotyped using ABI 7900HT Real-Time PCR system (Applied Biosystems: Thermo Fisher Scientific, Inc., Foster City, CA, USA).

Statistical analysis

The age and sex distributions between SCZ and control subjects in the two stages were tested by Student's t-test and Chi-square test, respectively. P<0.05 was considered statistically significant. The relative allele signal (RAS) score for each SNP was estimated as the ratio of signal intensity from each allele as analyzed by the R package SNP MaP. The differences in allele frequency between SCZ and controls were evaluated by the Z combined test in the pooled dataset. Association analysis was performed with PLINK v1.07 software (http://pngu.mgh.harvard.edu/~purcell/plink/) using additive, recessive and dominant models in the validation stage. We assumed that the prevalence of SCZ is 1% and used GAS Power Calculator to perform the power calculation. The sample size used in the validation stage could detect a risk allele with 73% power at a false-positive rate of 5%.

Integrative analysis of published GWAS data

The P-values for association between the two SCZ-related GWAS were derived from dbGaP database (http://www.ncbi.nlm.nih.gov/gap/) and downloaded for integrative analysis. The first study (dbGaP accession no. pha002857) consisted of 1,351 SCZ patients and 1,378 control subjects of European ancestry, whereas the second study (dbGaP accession no. pha002859) included 1,195 SCZ patients and 954 control subjects of American ancestry. Raw data of both studies were generated from Affymetrix GeneChip Human Mapping 6.0 chips. The detailed information of the two studies is shown in Table I. Moreover, SCZ meta-analysis data were obtained from the Psychiatric Genomics Consortium (https://www.med.unc.edu/pgc).

Table I.

Basic characteristics of European and American SCZ GWAS datasets.

			Genes
Study datasets	Study population	Sample size (case/control)	SNPs passing QC	With SNPs	P<10−5	P<10−2
pha002857	European	1,351/1,378	700,487	44,462	12	4,703
pha002859	American	1,195/954	809,854	44,840	4	5,863
Overlap		–	537,598	44,410	0	1,098

SCZ, schizophrenia; GWAS, genome-wide association studies.

Gene-SNP mapping

We downloaded the gene annotation file from the UCSC human genome annotation database (http://hgdownload.cse.ucsc.edu). The physical position of SNPs was obtained from the Affymetrix chip annotation file. The positional information of the genes and the SNPs was based on the human genome assembly GRCh37 (Feb. 2009). Based on the positional information, the SNPs identified in all three studies were mapped to the gene region and 20 kb upstream and downstream of gene boundaries to cover the coding sequence and most of the regulatory features with the PLINK v1.07 software (http://zzz.bwh.harvard.edu/plink/index.shtml). The SNP set P-value for a gene was calculated from the minimum P-value of the gene.

GSEA and pathway enrichment analysis

We performed functional analysis for all candidate genes using the GSEA method. The number of minimum candidate genes included in a pathway was set to 2 and it was considered to be significant with adjusted P<0.01. To test the repeatability of both GWAS data, the enrichment of one dataset was analyzed with the other dataset as the background. The pathway enrichment analysis was performed using the GEne SeT AnaLysis Toolkit (WebGestalt) (24). The pathway files were compiled from the KEGG pathway (http://www.genome.jp/kegg/pathway.html). P-value was calculated based on the hypergeometric distribution.

Results

Integrative analysis of the European and American SCZ GWAS datasets

We mapped SNPs to their corresponding genes to identify genetic factors related to SCZ in the European and American GWAS datasets. We used P<0.01 as the cutoff and identified 4,703 candidate genes in the European samples and 5,863 candidate genes in the American samples. Then, we integrated the two candidate gene lists and identified 1,098 common genes (~11.59%, 1,098/9,468), which was 1.79 times of the expected value of each individual pair-wise comparison (614 genes, χ2=471.4, df=1, P<10−4). This implied considerable replication between the two GWAS datasets for individual genes (Table I).

Moreover, we performed GSEA to confirm the replication between the two studies. We compared 4,703 candidate genes from the European study against the entire gene set from the American study, which were ranked according to their P-values. The candidate genes of the European study were enriched at the top of the gene list from the American study (Fig. 1A, ES=0.36390138, NES=1.7276831, nominal P<10−3). The GSEA of 5,863 candidate genes from the American study also showed analogous enrichment at the top of the European gene list (Fig. 1B, ES=0.381695544, NES=1.6385169, nominal P<10−3). These results suggested that the two SCZ GWAS datasets were highly overlapping and contained SCZ-related genes with similar biological function.

Figure 1.

Ranks of significant genes in the American and European GWAS. GSEA of candidate genes (A) in the European GWAS associated with SCZ in the American GWAS and (B) in the American GWAS associated with SCZ in the European GWAS. (C) The common SCZ-related pathways in the European and American GWAS. GWAS, genome-wide association studies; GSEA, Gene Set Enrichment Analysis; SCZ, schizophrenia.

Analysis of SCZ-related pathways

For better insight into the molecular mechanism of SCZ, we performed KEGG pathway enrichment analysis for candidate genes identified by both studies with the WebGestalt software (24). We identified 29 enriched KEGG pathways for the European study and 34 for the American study (adjusted P<0.01). As shown in Fig. 1C, 20 pathways (~83.33%, 20/34) were shared by the two studies.

Next, we compared the 20 common pathways to identify susceptibility genes in SCZ meta-analysis data from the Psychiatric Genomics Consortium. As shown in Fig. 2C, we validated six pathways, namely, axon guidance, vascular smooth muscle contraction, calcium signaling pathway, long term potentiation, MAPK signaling pathway and arrhythmogenic right ventricular cardiomyopathy (ARVC). In both studies, axon guidance, long term potentiation and ARVC were highly significant (P<10−3). In particular, the axon guidance represents a key stage in the formation of neuronal network and the pathway deserves our further study. Therefore, we focused on the axon guidance pathway in this study. The European, American and SCZ meta-analysis identified 23, 39 and 109 genes related to axon guidance pathway, respectively. Integrating the three datasets, we identified 15 axon guidance pathway genes that were shared among the three studies (Table II, Fig. 3).

Figure 2.

Ranks of significant genes in the pooled GWAS dataset. (A) Axon guidance pathway genes enriched in the pooled GWAS. (B) Fifteen significant axon guidance pathway-related genes identified in this study. (C) Common SCZ susceptibility pathways in the SCZ meta-analysis, European and American datasets. GWAS, genome-wide association studies; SCZ, schizophrenia.

Table II.

The 15 common axon guidance genes.

	-log10 (P-value)
Gene	European	American	SCZ meta-analysis
DCC	2.8091083	3.7437634	7.001612
DPYSL2	2.352617	2.5673513	3.7548857
EFNA5	3.3699791	2.5774102	4.387744
EPHA4	2.031097	2.1434524	3.0377238
EPHA5	2.103915	2.219971	5.9553895
EPHA6	2.111259	2.9742846	4.5623975
FYN	2.79588	2.2160964	5.4071918
LRRC4C	2.3809066	2.3208447	7.528907
PPP3CA	2.2432117	2.5166981	3.2800066
ROBO2	2.7068586	3.1870866	7.847146
SEMA6D	2.2464936	2.6424463	5.4734087
SLIT1	3.4405723	2.8052082	3.5275486
SLIT3	2.6044989	2.230475	6.288213
SRGAP3	3.5915904	2.2222834	4.070334
UNC5C	3.3293831	2.4387794	4.3944135

SCZ, schizophrenia; DCC, deleted in colorectal cancer; DPYSL2, dihydropyrimidinase-like 2; EFNA5, ephrin-A5; EPHA4, EPH receptor A4; EPHA5, EPH receptor A5; EPHA6, EPH receptor A6; FYN, FYN proto-oncogene, Src family tyrosine kinase; LRRC4C, leucine rich repeat containing 4C; PPP3CA, protein phosphatase 3 catalytic subunit α; ROBO2, roundabout guidance receptor 2; SEMA6D, semaphorin 6D; SLIT1, slit guidance ligand 1; SLIT3, slit guidance ligand 3; SRGAP3, SLIT-ROBO Rho GTPase activating protein 3; UNC5C, unc-5 netrin receptor C.

Figure 3.

The common axon guidance-related genes from European, American and SCZ meta-analysis datasets. SCZ, schizophrenia.

Validation of SCZ susceptibility genes in clinical samples

Next, we performed pooled GWAS on the clinical samples to verify the SCZ susceptibility genes involved in axon guidance in the Han population. GSEA pathway analysis identified the 15 susceptibility genes in the pooled GWAS data (Fig. 2A). The 15 SCZ susceptibility genes were deleted in colorectal cancer (DCC) (P=9.96e-08), dihydropyrimidinase-like 2 (DPYSL2) (P=1.75e-04), ephrin-A5 (EFNA5) (P=4.09e-05), EPH receptor A4 (EPHA4) (P=0.00092), EPH receptor A5 (EPHA5) (P=1.11e-06), EPH receptor A6 (EPHA6) (P=2.74e-05), FYN proto-oncogene, Src family tyrosine kinase (FYN) (P=3.92e-06), leucine rich repeat containing 4C (LRRC4C) (P=2.96e-08), protein phosphatase 3 catalytic subunit α (PPP3CA) (P=0.00052), roundabout guidance receptor 2 (ROBO2) (P=1.42e-08), semaphorin 6D (SEMA6D) (P=3.37e-06), slit guidance ligand 1 (SLIT1) (P=0.0003), slit guidance ligand 3 (SLIT3) (P=5.15e-07), SLIT-ROBO Rho GTPase activating protein 3 (SRGAP3) (P=8.51e-05) and unc-5 netrin receptor C (UNC5C) (P=4.04e-05). The pathway enrichment of these 15 axon guidance genes was significant with P<10−3 (Fig. 2B). DCC, LRRC4C and ROBO2 were the top genes in the list with P<10−7.

We further investigated these 15 genes in the pooled GWAS and SCZ meta-analysis datasets, and identified 7,237 common loci related to the 15 genes. After Bonferroni correction, we obtained 14 significant loci in four genes, DCC, FYN, ROBO2 and SEMA6D (Table III). Among the 14 candidate loci, two SNPs, rs9944880 in DCC gene and rs7630091 in ROBO2 gene were significant in both the meta-analysis and pooled GWAS datasets (P<10−4).

Table III.

The SCZ-associated loci after Bonferroni correction.

Gene	Chromosome	SNP	Meta P-value	Pooled P-value
DCC	18	rs9944880	4.19E-06	5.55109E-05
FYN	6	rs7764009	1.86E-04	0.375738765
ROBO2	3	rs7636341	1.42E-05	0.447750198
ROBO2	3	rs1516457	1.63E-05	0.163967215
ROBO2	3	rs11128505	1.98E-05	0.013580237
ROBO2	3	rs7633954	2.65E-05	0.526997858
ROBO2	3	rs7630091	4.14E-05	3.95512E-05
ROBO2	3	rs3849481	5.69E-05	0.086994752
ROBO2	3	rs6802433	8.72E-05	0.657963545
SEMA6D	15	rs1559677	1.02E-05	0.00947373
SEMA6D	15	rs281299	1.24E-05	0.019215293
SEMA6D	15	rs1353580	1.24E-05	0.006132857
SEMA6D	15	rs13313462	1.44E-05	0.163967215
SEMA6D	15	rs1948572	1.51E-05	0.899343456

Meta P-value was obtained from SCZ meta-analysis data of the Psychiatric Genomics Consortium. Pooled P-value was calculated from the pooled GWAS data. SCZ, schizophrenia; DCC, deleted in colorectal cancer; FYN, FYN proto-oncogene, Src family tyrosine kinase; ROBO2, roundabout guidance receptor 2; SEMA6D, semaphorin 6D; GWAS, genome-wide association studies.

We validated rs9944880 and rs7630091 in an independent pooled dataset that included 1,731 SCZ patients and 1,849 health controls. The rs9944880 located in DCC gene was associated with SCZ patients (Table IV). The T allele was associated with SCZ risk under additive model (P=0.004, OR=1.19, 95% CI=1.06–1.33). Individuals with rs9944880TT genotype were more likely to contract SCZ under recessive model (P=0.0001, OR=1.67, 95% CI=1.29–2.17 for TT vs. C+CT). No significant association was found in dominant model (P=0.218, OR=1.10, 95% CI=0.95–1.26). However, the SNP rs7630091 in ROBO2 was not associated with SCZ risk in our validation cohort under the three models.

Table IV.

Verification of gene loci in SCZ patients from Shandong migrants.

		Genotype count			Allele frequency
SNP	Group	AA	AB	BB	A	B	Padd	Pdom	Prec
rs9944880		CC	CT	TT	C	T	0.004	0.218	0.0001
	SCZ	1,167	410	154	0.793	0.207
	Controls	1,283	464	102	0.819	0.181
rs7630091		GG	AG	AA	G	A	0.680	0.693	0.134
	SCZ	913	602	216	0.701	0.299
	Controls	962	687	200	0.706	0.294

For the SNP rs9944880, the symbol A indicates allele C and B indicates allele T. For the SNP rs7630091, the symbol A indicates allele G and B indicates allele A. All SNPs were analyzed under additive (add), dominant (dom) and recessive (rec) models. SCZ, schizophrenia.

Discussion

SCZ is a complex polygenic disease, which is associated with many genes and non-coding loci (25). Traditional GWAS analysis has mainly focused on the most significant SNPs that are associated with the pathogenesis of SCZ (26). Therefore, the less significant genes have not yet been analyzed. Therefore, we developed a comprehensive approach by combining a sub-significance strategy with pathway analysis. When we used stringent common genome-wide significance level of P<10−5, only 12 and 4 candidate genes were identified in the European and American studies, respectively. Moreover, none of these genes were common to either group. This suggested poor repeatability between the two studies or a substantial loss of information due to stringent criterion. Therefore, when the threshold was set as P<10−2, we identified 4,703 and 5,863 candidate genes in the European and American studies, respectively, and 1,098 genes (12%) were common in the two gene lists. Wang et al performed a similar analysis with low significance standard to identify novel susceptibility genes in the coronary artery disease (27). Moreover, the independent GSEA for European and American studies showed 47% (20/43) common pathways between the two studies. Taken together, these results suggest the feasibility of adopting a less stringent significance level combined with pathway analysis to identify molecular mechanisms of polygenic diseases.

Based on this strategy, we found high degree of repeatability in the genes identified by the two studies. Thus, we identified 20 candidate pathways associated with SCZ in the two independent SCZ GWAS databases (European and American) from NCBI. We also verified the top three significantly enriched pathways (P<10−4) through pathway analysis of SCZ meta-analysis dataset.

The relationship of the axon guidance pathway with SCZ has not been well studied. The aging of human cerebral cortex depends on the proliferation and differentiation of neural progenitor cells, and proper migration and positioning behavior of neurons (28). Disruption in any of these processes can result in aberrant development of the cerebral cortex and is central to many mental illnesses such as SCZ. However, further investigations are necessary to prove this aspect beyond doubt.

In our study, we identified two core genes, ROBO2 and DCC, in the axon guidance pathway. Axons integrate directional information multiple guidance cues and their receptors during development. ROBO2 and DCC are axon guidance receptors that can function individually or in combination with other guidance receptors to regulate downstream effectors (29,30).

DCC is a membrane receptor for Netrin-1, which is mainly expressed in various types of axon beams of the developing central nervous system (CNS) (31). DCC regulates axon guidance, projections, branching and a number of important physiological processes of the CNS neurons after binding to Netrin-1 (32–34).

ROBO2 is a membrane protein that belongs to a class of conserved and development-related transmembrane receptor family roundabout, which promotes axon guidance of embryonic CNS, dendritic branching, axonal growth and neuronal cell migration (35,36). ROBO2 cooperates with Slit proteins to guide the major forebrain projections and mainly regulates axon guidance (37).

DCC and ROBO2 as well as their respective ligands, Netrin-1 and Slits, guide the callosum projection neurons across the midline (38–42). Glial wedge and indusium griseum represent the main glial cell populations located near the center, which secrete nerve guide factors such as Netrin-1 and Slits. These neural guidance factors interact with their corresponding membrane receptors, DCC and ROBO2, on the neuronal axons projected by the corpus callosum. These processes ensure the axons project correctly and finally form corpus callosum (38–42).

In this study, we demonstrated that rs9944880 polymorphism in DCC gene on chromosome 18q21.2 was associated with SCZ. From a mechanistic point of view, this polymorphism probably altered the axon guidance pathway resulting in SCZ. In conclusion, our study showed that aberrant axon guidance pathway was central to SCZ pathogenesis by combining pathway analysis with less stringent significance criterion in analyzing gene sets from two GWAS. This study needs to be confirmed in large cohorts of SCZ patients.