Genome-Wide Association Study of Treatment Refractory Schizophrenia in Han Chinese

We report the first genome-wide association study of a joint analysis using 795 Han Chinese individuals with treatment-refractory schizophrenia (TRS) and 806 controls. Three loci showed suggestive significant association with TRS were identified. These loci include: rs10218843 (P = 3.04×10−7) and rs11265461 (P = 1.94×10−7) are adjacent to signaling lymphocytic activation molecule family member 1 (SLAMF1); rs4699030 (P = 1.94×10−6) and rs230529 (P = 1.74×10−7) are located in the gene nuclear factor of kappa light polypeptide gene enhancer in B-cells 1 (NFKB1); and rs13049286 (P = 3.05×10−5) and rs3827219 (P = 1.66×10−5) fall in receptor-interacting serine/threonine-protein kinase 4 (RIPK4). One isolated single nucleotide polymorphism (SNP), rs739617 (P = 3.87×10−5) was also identified to be associated with TRS. The -94delATTG allele (rs28362691) located in the promoter region of NFKB1 was identified by resequencing and was found to associate with TRS (P = 4.85×10−6). The promoter assay demonstrated that the -94delATTG allele had a significant lower promoter activity than the -94insATTG allele in the SH-SY5Y cells. This study suggests that rs28362691 in NFKB1 might be involved in the development of TRS.


Introduction
Schizophrenia is a severe psychiatric disorder with a prevalence estimated to be approximately 1% [1] in the world and 0.6% in Taiwan [2]. It is the third-leading cause of disability among individuals age between 15 and 44 [3]. Its clinical manifestations are characterized by distortion of reality, delusions, hallucinations, altered emotional reactivity, disorganized behavior, social isolation and cognitive impairment. The etiology of schizophrenia is not well understood but it has been postulated as a complex disease with an estimated heritability as high as 80% [4] [5].
Antipsychotic medication is the major treatment for schizophrenia. However, one fifth to one third of schizophrenic patients do not respond to antipsychotic treatments [28,29]. These patients with treatment refractory schizophrenia (TRS) have persistent psychotic symptoms combining with poor social/work function in spite of administering at least two trials of sufficient antipsychotic doses and adequate treatment duration [28]. Comparing with those patients with adequate responses to antipsychotic treatments, patients with TRS had significantly lower levels of catecholamine in cerebrospinal fluid or plasma [30], increased cortical atrophy [31,32], and a lower level of plasma tryptophan concentrations [33]. Therefore, TRS may be a distinct and homogenous subgroup of schizophrenia.
To identify the genetic variants susceptible for schizophrenia, this study performed the first GWAS focusing on TRS in a Han-Chinese population. We identified several novel genetic loci which were not associated with schizophrenia. Our findings may pave a new way to elucidate the underlying molecular mechanism of schizophrenia and to improve the treatment for TRS.

Demographic information
Demographic data from 522 TRS patients and 806 controls is listed in Table

Data quality
The average call rate was 99.860.3% for each subject genotyped in this study. Gender determined from the GWAS result for all the subjects were consistent with recorded data. 694,436 (79.99%) of the 868,114 SNP in the autosomes passed the quality control filter and had an average call rate of 99.860.4% (Supplementary Table S1). The results of principal component analysis showed no significance for population stratification between TRS patients and controls, (P.0.05, and Fst statistics between populations ,0.001) (Supplementary Figure S1). Furthermore, genomic control with a variance inflation factor l = 1.042 (trend test), estimated posterior to the regular GWAS, also indicated no substantial population stratification. These SNPs were then taken for further GWAS analysis.

Association analysis
Data analysis was first performed for the 522 TRS patients and 806 controls (Figure 1). Preliminary results revealed 19 SNPs with suggestive significant associations with TRS (Supplementary Table  S2, 10 28 ,P,10 25 ). Fourteen markers were retained after cross platform validation with the Sequenom platform and showed a concordance rate of over 98% (Supplementary Table S2).
Four major clusters with more than one SNPs located within 500 kb of each other were identified from the 14 validated SNPs (Table 2, Figure 2). The first locus, comprising rs10218843 (P = 6.73610 26 ) and rs11265461 (P = 5.90610 26 ), and is located approximately 10 kb downstream of signaling lymphocytic activation molecule family member 1 (SLAMF1) on chromosome 1; the second locus contains two SNPs, rs4699030 (P = 8.41610 27 ) and rs230529 (P = 1.07610 26 ), is located in the introns of nuclear factor of kappa light polypeptide gene enhancer in B-cells 1 (NFKB1) on chromosome 4; three SNPs, rs739617 (P = 1.46610 25 ), rs17158926 (P = 3.99610 25 ) and rs17158930 (P = 3.08610 25 ) are clustered in the introns of dedicator of cytokinesis 4 (DOCK4) on chromosome 7; and the last locus which consists of two SNPs, rs13049286 (P = 1.23610 25 ) and rs3827219 (P = 1.23610 25 ), is located in receptor-interacting serine/threonine-protein kinase 4 (RIPK4) on chromosome 21. These loci are located in the regions with high LD (except for chromosome 21) (Supplementary Figure S2). Multipoint/Haplotype analysis also showed that these clusters were associated with TRS, the cluster on chromosome 4 had the highest P value with global score P = 2610 25 (Supplementary Table S3).
In addition to the above SNPs in clusters, five other SNPs also showed suggestive significant association with TRS. These SNPs are: rs461409 (P = 2.      Figure S3). Except for rs13049286 and rs3827219 with odds ratio (OR) of approximately 3, all other SNPs identified in this study showed modest effects with OR between 1.06-1.81 (Table 2).

Joint Analysis in with additional TRS patients
The 14 SNPs showing suggestive significance were then genotyped in an independent cohort of 273 TRS patients. An average call rate of 99.3760.21% was achieved for each subject. Joint analysis was then carried out in the 795 cases and 806 controls. Of the 14 SNPs showing suggestive association in the initial analysis, 7SNPs remain suggestively associated with TRS after joint analysis ( Table 2). These SNPs are: rs10218843 (P joint = 3.04610 27 ) and rs11265461 (P joint = 1.94610 27 ), which both are located in SLAMF1; rs230529 (P joint = 1.47610 27 ) and rs4699030 (P joint = 1.94610 26 ) are located in NFKB1; rs739617 (P joint = 3.87610 25 ) is in DOCK4; rs13049268 (P joint = 3.05610 25 ) and rs3827219 (P joint = 1.66610 25 ) which are located in RIPK4.

Testing in schizophrenic patients
The top SNPs showing suggestive significance were then tested in an independent cohort of 1982 schizophrenic patients whose responses to antipsychotic treatments were not determined and additional 2000 controls. An average call rate of 99.3760.21% was achieved for each subject. However, none of these SNPs were significantly associated with this group of patients (Supplementary Table S4), suggesting that these SNPs were specifically associated with TRS and not a broad phenotype of schizophrenia.

Re-sequencing of NFKB1
Because the lowest P values in both single and multi-point analysis were observed for the SNPs located on NFKB1 on chromosome 4, we next aimed to identify variants with functional consequence in NFKB1. Re-sequencing was performed on the exons, intron-exon boundaries, and a 2-kb region covering the promoter of NFKB1 in a discovery cohort of 94 TRS patients and 94 controls. Twenty-three genetic polymorphisms including 11 novel variants and 2 non-synonymous changes (R231L and R534H) were identified in NFKB1 (Table 3). The rs28362491 SNP with an ATTG deletion in the promoter region of NFKB1 (94delATTG) was reported to affect nuclear protein binding and gene transcription in colonic epithelial cells [34]. In a test with 520 TRS cases and 806 controls, rs28362491 was associated with TRS (P = 6.69610 25 ). rs28362491 is in linkage disequilibrium with rs230529 and rs4699030 (r 2 = 0.741 and 0.714, respectively) ( Figure 3).

Functional analysis of rs28362491
Since rs28362491 has been reported to affect nuclear protein binding and gene transcription in colonic epithelial cells, it is possible that the deletion also alters the efficiency of transcription in neuronal cells. The promoter assay showed that the construct containing the -94delATTG promoter displayed significantly

Discussion
This study presented the data of the first GWAS on TRS in a Han-Chinese population. Therefore, the use of genomic control did not substantially change the results of this GWAS (Supplementary Table S5). The sample size is comparatively smaller than in previous GWAS for schizophrenia due to the smaller number of TRS patients available Since schizophrenia is an etiologically heterogeneous disorder, narrowing schizophrenia down to TRS may represent a more discrete and genetically homogeneous group to identify genes involved in the etiology of schizophrenia.
Among the genetic loci with suggestive significance identified by this study, three regions with more than one significant SNP in each region stood out after joint analysis. rs4699030 and rs230529 are located in the introns of NFKB1. This gene encodes for two functional different proteins [35], one for the 105 kD (p105) protein, and the other for a 50 kD (p50) protein. P105 contains seven copies of ankyrin-like sequence in the carboxyl terminal region which is similar to those in I-kappa B kinase (IkB), therefore, p105 may also have functions similar to IkB. P105 could associate with either c-Rel or RelA in the cytoplasm to inhibit Rel protein-specific transcription [36,37]. P50 is one subunit of NFkappa B with repression activity. NF-kappa B is well distributed and a highly conserved dimeric transcriptional factor which regulates more than 200 genes [38]. Different dimeric combinations of NF-kappa B are found in different tissues and respond differently to regulate gene expressions. P50 assembles either with other NF-kappa B subunits, such as RelA, RelB, c-Rel, or with itself as a homodimer to repress NF-kappa B-dependent gene transcriptions [39]. The heterodimer of p50 and RelA subunit is the most abundant form of NF-kappa B [40].
Both rs10218843 and rs11265461 are located near SLAMF1 (also known as CD150), which is a signaling lymphocyte activation molecule and a member of the CD2 family belonging to the immunoglobulin superfamily of receptors. SLAM is a costimulatory molecule involving T cell activation and also as a receptor for measles virus [41], a bacterial sensor [42], and responsible for the NKT lineage development [43].The activation of SLAM has been shown to associate with numerous distinct downstream activities, including augmenting T cell mediated cytotoxicity through a sequence of signaling transduction, including NF-kappa B activation, Stat1 phosphorylation and Tbet induction [44]; increasing recruitment of protein kinase C (PKC)-theta and the activation of NF-kappa B p50, both of which are involved in enhancing T helper 2 cytokines production and natural killer T cell development [45,46]. Two SNPs, rs13049286 and rs3827219 are both located in RIPK4. The expression of RIPK4 is well distributed, including the brain, and found to interact with PKC-delta [47]. This gene is important in sensing cellular stress, such as infection, inflammation, cellular differentiation programs and DNA damage. It also mediates downstream signaling, in particular the activation of NFkappa B and the induction of apoptosis [48].
The identification of NFKB1, RIPK4, and SLAMF1 in this study suggests that the NF-kappa B pathway plays an important role in the pathogenesis of TRS. NF-kappa is also found to be related with schizophrenia since the genetic variants in RELA gene, which encoded the protein RelA as one subunit of NF-kappa B, were reported to be associated with schizophrenia and the patients startle responses in a Japanese population [49]. NF-kappa B is a  key transcriptional factor in the regulation of the expression of many inflammatory factors, such as cytokines, chemokines, and adhesion molecules [50]. Several studies have demonstrated that schizophrenic patients' cerebrospinal fluid and plasma had abnormal levels of cytokines [51,52,53], the aberrations were especially more pronounced in TRS [54,55,56]. Song et al. found the elevated level of cytokine in first-episode schizophrenic patients was associated with the activation of NF-kappa B [53]. Thus, these suggest abnormal inflammatory response could lead to TRS. One common variant rs28362691 (-94ins/delATTG) identified from resequencing NFKB1 was found to associate with TRS. This SNP is located in 19 base pairs upstream of a functional kB binding site [57]. The promoter assay showed that the NFKB1 promoter with the -94delATTG allele had a lower promoter activity in SH-SY5Y cells in comparison with the -94insATTG allele. This implies that the -94delATTG allele may result in lower expression of NFKB1. Changes in NFKB1 expression could alter the level of p105 and induce divergent dimeric combinations of NF-kappa B, which might cause disturbances in cytokine regulations, and lead to a failure of antipsychotic treatment. However, the association between the abnormal levels of cytokines and NF-kappa B in patients with TRS remains to be established. Two novel non-synonymous polymorphisms (R231L and R534H) were also identified from resequencing. However, these polymorphisms have extremely low frequency in the Han population and their effects on NFKB1 function remain to be elucidated.
Other genetic loci identified in association with TRS in this study suggests genes involved in neuronal development (DOCK4 [58]). A recent study conducted in the Jewish population also identified a SNP (rs2074127) in DOCK4 associated with schizophrenia [59] however, this SNP is not in LD with the DOCK4 SNPs reported in this study (Supplementary Figure S4). Its role in the development of TRS remains to be elucidated.
We also compared our results with previous genetic studies showing associations with TRS (such as variants in CYP3A4 [60], CYP3A5 [60], DRD3 [60], HTR2A [61], HTR3A [62]), none of the variants or their nearby SNPs had a P value lower than 10 25 in this study. It implicated that these above genes involving in metabolic enzymes and receptors of antipsychotics were not associated with TRS. Furthermore, we also compared our data with previous GWAS on schizophrenia. We examined the P values of the loci previously reported to be associated with schizophrenia in our data. We also checked the vicinity (200 kb) of the reported loci in our data. Only rs1602565 on chromosome 11 showed nominal association (P = 3.17610 24 ) in this study (Supplementary  Table S6 and Supplementary Figure S5). These data suggested that the genetic loci identified in this study were specifically associated with TRS.
None of the loci identified in this study reached genome-wide significance, this could due to tour sample size lack the statistical power to detect common variants of susceptibility. However, we focused on TRS within schizophrenia, which may represent a more homogeneous group. An independent TRS group was also not available for replication study. Future replication studies in additional population of TRS are required.
In conclusion, we report the first GWAS on TRS in the Han Chinese population. Our data suggest that the NF-kappa B pathway may play an important role in the pathogenesis of TRS. We have also provided the functional effect of the -94delATTG allele showing the possible mechanism of NFKB1 in TRS. Further studies are required to confirm the association of the risk alleles identified in this study across different populations, to identify the causative genetic variants, and to elucidate the underlying molecular mechanisms, which may help to improve treatments for refractory schizophrenia.

Subjects
This study was approved by the Institute Review Board of Taipei Veterans General Hospital Kaohsiung, Kai-Suan Psychiatric Hospital, Jianan Mental Hospital, Bali Psychiatric Center, Tsyr-Huey Mental Hospital, Yuli Veterans Hospital, National Taiwan University Hospital and Academia Sinica. Written informed consent was obtained from all the study participants.
A total of 522 unrelated patients with TRS, including 289 males (55.4%) and 233 females (44.6%), were recruited from Yuli Veterans Hospital, Taipei Veterans General Hospital, Kaohsiung Kai-Suan Psychiatric Hospital, Tsyr-Huey Mental Hospital, Jianan Mental Hospital, and Bali Psychiatric Center. DNA samples from additional 273 TRS patients were obtained from National Taiwan University Hospital and were used in joint analysis. In addition, DNA samples from 1982 schizophrenic patients were obtained from Yuli Veterans Hospital, Taipei Veterans General Hospital, and National Health Research Institutes. However, the responses to antipsychotic treatments for 1982 schizophrenic patients were not determined.
All patients were diagnosed according to the criteria of DSM-IV for schizophrenia. TRS was defined using a modified Conley and Kelly's protocol [28]. Briefly, schizophrenic patients with the following criteria were identified as TRS: No improvement in clinical impression (defined as 3 or more in the global improvement subscale of clinical global impression (CGI-I)) after at least two six-week trials of antipsychotic therapy at a dose equal to or higher than the equivalent daily dose of 600 mg of chlorpromazine for typical antipsychotics, or for second-generation antipsychotics (risperidone: 6 mg/day; olanzapine: 20 mg/ day; quetiapine: 800 mg/day; ziprasidone: 160 mg/day; amisulpride: 800 mg/day; zotepine: 300 mg/day), as well as for patients who were administered the last-line antipsychotic pharmacotherapy, clozapine (50-300 mg/day). All patients with TRS showed more than 5 years of persistent illness (defined as 4 or more in the severity of illness subscale of clinical global impression (CGI-S)). Informed consent was obtained from all participants. Only the Han-Chinese population, which accounts for 98% of the population in Taiwan, was recruited for this study.
The control (N = 2806) was randomly selected from the Han-Chinese Cell and Genome Bank in Taiwan described previously [63], in which more than 3,300 controls were collected and randomly selected through registry.

Genotyping and Quality Control
Genomic DNA was isolated from peripheral blood using PUREGENE DNA purification system (Gentra Systems, Minneapolis, MN). Whole-genome scan was conducted using Affyme-trixH Genome-wide Human SNP Array 6.0 (Affymetrix, Santa Clara, CA, USA) on 522 TRS patients and 806 controls and genotyping was performed by National Genotyping Center at Academia Sinica. Genotype calling was determined by Affymetrix Power Tool 1.10.2 (Affymetrix) using default parameters.
Quality control of genotype data was performed by examining several summary statistics. First, individual's gender was double checked by calculating the ratio of loci with heterozygous calls on the X chromosome; After calculating total successful call rate and the minor allele frequency (MAF) of cases and controls for each SNP, SNPs were excluded if one of the following conditions applied: (1) only one allele appeared in both cases and controls; (2) the total call rate was less than 98%; (3) the total MAF was less than 5% and the total call rate was less than 99%; (4) significant deviation from Hardy-Weinberg equilibrium in the control group (P,10 24 ).

Population stratification
Detection of possible population stratification that might influence association analysis was carried out using principal component analysis (supplementary Methods S1, Supplementary  Table S7) with genotype data for 100,000 SNPs located at equally spacing across the human genome. Plink (Supplementary Methods S2) was performed to examine if the subjects were related with each other. Variance inflation factor for genomic control was estimated based on all qualified SNPs (Supplementary Methods S3).

Data Analysis
Genotyping data analysis was carried out by comparing the frequencies of allele and genotype between cases and controls using the following single-point methods: genotype, allele-type, trend test (Cochran-Armitage test), dominant, and recessive models. Empirical P-values were also obtained with 10 8 simulations. SNPs with P-values less than a = 10 28 , a cut-off for the multiple-comparison adjusted by Bonferroni correction, were considered to be significantly associated with the traits. SNPs with P-values between 10 28 and 10 25 were considered to have suggestive significant association. Quantile-quantile (Q-Q) plots were then used to examine P-value distributions (Supplementary Figure S6).
Multi-point/haplotype analysis was performed using the Haplotype Score Test [64] implemented in haplo.stat, a suite of S-PLUS/R routines for the analysis of indirectly measured haplotypes, for regions with more than two SNPs having genetic evidences (P-value,10 25 ). Regions were tested with independent 10 6 simulations.

Validation and joint analysis
Autosomal SNPs with P-value,10 25 from GWAS in 522 TRS cases and 806 controls were further validated using MALDI-TOF mass spectrometry (SEQUENOM MassARRAY, Sequenom, San Diego, CA, USA). The SNPs retained after cross-platform validation were then genotyped in the additional 273 TRS cases. Joint analysis was then carried out with all the 795 TRS cases and 806 control.

Direct Sequencing
Selected candidate genes were re-sequenced in a discovery cohort consisted of 94 TRS patients and 94 controls. Exons, 200 bp of exon/intron junctions, and a 2-kb region covering the promoter of the selected genes were sequenced using Applied Biosystems 3730 (CA, USA). Contig assembly and SNP identification were determined using Sequencher 4.5 Demo (Gene Codes Cooperation, Ann Arbor, MI, USA). All PCR products were bidirectionally sequenced.

Plasmid Construction for Luciferase Reporter Assay
To assay for the NFKB1 promoter activity, the NFKB1 promoter encompassing the -94ATTG polymorphism (from 21000 to 21) from patients with homozygous -94ATTGATTG (W) or -94ATTG (D) were first cloned into the pGEM-T Easy vector (Promega, Madison, WI, USA) with the forward primer: 59-CCCGGGCCTGATTACTGATGTTTTAAAGCT-39 and the reverse primer: 59-CTCGAGTTCCTGGCTGGAAATTCCCA-CTGA-39. Both the W and D fragments were then released from the pGEM-T Easy vector and subcloned into the upstream region of the firefly luciferase gene of the pGL3-basic vector (Promega). All constructs were subjected to sequencing to confirm the orientation and integrity.

Transient Transfection and Luciferase Assay
A total of 1610 5 SH-SY5Y cells were seeded in each well of a 24-well plate and transfected with 450 ng of each reporter construct along with 50 ng of pRL-TK vector (Promega) containing the Renilla luciferase gene as an indicator for normalization of transfection efficiency. Transfections were performed by using FuGENE HD (Roche Applied Science, Indianapolis, IN, USA) according to manufacturer's instructions. Cells were incubated for 24 hours and analyzed for luciferase activity with the Dual-Luciferase Assay System (Promega). Firefly luminescence was normalized to Renilla luminescence and reported as relative luciferase activity. All experiments were performed in triplicate and repeated at least three times.