Genome-wide association study in obsessive-compulsive disorder: results from the OCGAS.

Obsessive-compulsive disorder (OCD) is a psychiatric condition characterized by intrusive thoughts and urges and repetitive, intentional behaviors that cause significant distress and impair functioning. The OCD Collaborative Genetics Association Study (OCGAS) is comprised of comprehensively assessed OCD patients with an early age of OCD onset. After application of a stringent quality control protocol, a total of 1065 families (containing 1406 patients with OCD), combined with population-based samples (resulting in a total sample of 5061 individuals), were studied. An integrative analyses pipeline was utilized, involving association testing at single-nucleotide polymorphism (SNP) and gene levels (via a hybrid approach that allowed for combined analyses of the family- and population-based data). The smallest P-value was observed for a marker on chromosome 9 (near PTPRD, P=4.13 × 10(-)(7)). Pre-synaptic PTPRD promotes the differentiation of glutamatergic synapses and interacts with SLITRK3. Together, both proteins selectively regulate the development of inhibitory GABAergic synapses. Although no SNPs were identified as associated with OCD at genome-wide significance level, follow-up analyses of genome-wide association study (GWAS) signals from a previously published OCD study identified significant enrichment (P=0.0176). Secondary analyses of high-confidence interaction partners of DLGAP1 and GRIK2 (both showing evidence for association in our follow-up and the original GWAS study) revealed a trend of association (P=0.075) for a set of genes such as NEUROD6, SV2A, GRIA4, SLC1A2 and PTPRD. Analyses at the gene level revealed association of IQCK and C16orf88 (both P<1 × 10(-)(6), experiment-wide significant), as well as OFCC1 (P=6.29 × 10(-)(5)). The suggestive findings in this study await replication in larger samples.

Introduction

Obsessive-compulsive disorder (OCD) is a psychiatric condition characterized by persistent, intrusive, senseless thoughts and urges (obsessions) and repetitive, intentional behaviors (compulsions). Affected individuals tend to recognize that their thoughts and behaviors are excessive and unreasonable, and often struggle to resist them. The lifetime prevalence of OCD is estimated to be between 1–3%, based on national and international population-based surveys [1, 2]. Patients experience a chronic or episodic course with exacerbations that can substantially impair social, occupational, and academic functioning; according to the World Health Organization, OCD is among the most disabling medical conditions worldwide [3].

Since the early twentieth century, clinicians have suspected that heredity plays an important role in OCD susceptibility. Consistent with this, increased OCD prevalence was identified among first-degree relatives of case probands (11.7%), compared to relatives of controls (2.7%; the Hopkins OCD family study) [4]. In the same study, the prevalence of OCD in siblings of early onset probands was 17.9% (λsib =7.8) and thus within the range found in other psychiatric disorders for which a heritable component was reported, such as bipolar disorder and panic disorder. Family studies of child and adolescent OCD probands report prevalence rates of 7% to 15% in first-degree relatives, and 20–30% of these probands had one or more first-degree relative with OCD [5]. These findings are consistent with previous reports of an increased familial loading in probands with early age at onset. In addition, a recent population-based study suggested that (while there is no significant effect of a shared environment on OCD risk) the contribution of genetic factors to OCD risk is as high as 50% [6].

Two genetic linkage studies of OCD have found a suggestive linkage peak on chromosomal region 9p24 [7, 8]. The OCD Collaborative Genetics Study (OCGS) also found suggestive genetic linkage peaks on chromosomal regions 3q27-28, 6q, 7p, 1q, and 15q, but none were genome-wide significant [9].

There have been numerous candidate gene association studies in OCD [10]. A positional candidate of particular interest is the neuronal and epithelial glutamate transporter gene (SLC1A1). It is located in the region of the replicated linkage peak on chromosome 9p24, and there is considerable evidence from imaging, animal, and treatment studies that abnormal glutamatergic transmission may be involved in the pathophysiology of OCD [11, 12]. Several studies have reported positive associations of OCD with this gene, but with different SNPs in each [13-19]. A recent meta-analysis of these studies did not resolve the issue and was not able to report experiment-wise significance for SNPs annotated to SLC1A1 [20].

A genome-wide association study was recently completed by the International OCD Foundation Genetic Collaborative (IOCDF-GC) [21] with a combined study sample from 22 sites, comprising 1 465 cases and 5 557 ancestry-matched controls, as well as 400 complete trios. In case-control analyses, the lowest P-values were found for two SNPs in perfect LD (rs11081062 and rs11663827; r2=1, D′=1 in reference samples from the 1000 genomes project [22]; P=2.49×10−6 and P=3.44×10−6) Both were located on chromosome 18 within the discs, large (Drosophila) homolog-associated protein 1 gene (DLGAP1), a member of the neuronal postsynaptic density complex. In the trio analysis, a SNP near the BTB (POZ) domain containing 3 gene (BTBD3; rs6131295, chromosome 20), exceeded the genome-wide significance threshold with a P-value=3.84×10−8. However, when trios were meta-analyzed together with the case-control samples, the P-value for this variant was 3.62×10−5, losing genome-wide significance. Although no SNPs were found to be associated with OCD at a genome-wide significant level in the combined trio–case–control sample, a significant enrichment of methylation QTLs (P<0.001) and frontal lobe expression quantitative trait loci (eQTLs) (P=0.001) was observed within the top-ranked SNPs (P<0.01) from the trio–case–control analysis, potentially pointing to a broad role in gene expression in the brain, and possibly in the etiology of OCD.

In this paper, we report the results of the most recent genome-wide association study of OCD, the OCD Collaborative Genetics Association Study (OCGAS). Investigators at 8 research centers in the United States, including Brown University, Columbia University, Johns Hopkins University, Massachusetts General Hospital, Harvard School of Public Health, the National Institute of Mental Health, University of California Los Angeles, and the University of Southern California, conducted this collaborative study. The aim of the study was to identify common variants associated with OCD using an integrative analyses pipeline that involved association testing at both SNP and gene levels. The approach included OCD-affected patients (in a family-based analytic setting) with an early onset of the disorder (eighteen years and younger) who had been comprehensively assessed via an identical diagnostic approach at all recruitment centers.

Materials and methods

Sample

A total of 1 065 families were included in this study (comprising 1 406 patients with OCD and 2 895 individuals in total); 621 families were recruited and assessed specifically for this study at one of the five participating recruitment sites or the National Institute of Mental Health; 444 families had previously been evaluated in one of the earlier studies at Hopkins University or by one of the collaborating sites. The sample comprised of 460 complete trios (including an affected proband and both parents); 155 pedigrees with a proband and an unaffected sibling and 450 families with another structure (complex family structure). An additional 192 probands without an additional family member present in the study (singletons) were included. A breakdown of the families by site is found in Supplementary Table S1.

For study inclusion, probands were required to meet DSM-IV criteria for OCD [23] with onset of obsessions and/or compulsions before the age of 18 years (mean = 9.4 years; SD=6.35). Subjects disease, schizophrenia, severe mental retardation that does not permit an evaluation to characterize the psychiatric disorder, Tourette disorder (TS), or OCD occurring exclusively in the context of depression (secondary OCD) were excluded. In addition, individuals were removed from the sample if they were previously diagnosed with brain pathology including brain tumors, Huntington’s Disease, Parkinson’s Disease, or Alzheimer’s Disease. Each case was evaluated by a PhD-level clinical psychologist using the Structured Clinical Interview for DSM-IV (SCID) modified and extended to include additional symptom and diagnostic information as indicated in the supplementary material. Final diagnostic status was assigned based on the consensus of two psychiatrists or psychologists reviewing the case independently. Both parents of the proband were also recruited whenever possible. When parents were unavailable for participation, unaffected siblings were recruited. Genotyping was performed at the Johns Hopkins SNP Center using Illumina’s HumanOmniExpress bead chips (Illumina Inc., San Diego, CA, USA). More details on the diagnostic assessment and the genotyping process are provided in the supplementary materials and methods.

To increase the power of the study to detect significant association, we also included 1 984 unrelated controls (genotyped with Illumina’s HumanOMNI1-QUAD bead chip) from a previously published study on Parkinson’s disease (dbGaP accession number phs000196.v2.p1) [24, 25]. Individuals with a self-reported or diagnosed neuropsychiatric disorder at the time of enrollment were excluded from the present study.

Quality control of genotyping data

We followed a stringent quality control (QC) protocol that was designed to minimize occurrence of false positive signals, which included checking the relatedness of samples (i.e. verifying the relationship reported from the participating clinical centers) and reported sex (based on deviations from expected heterozygosity rates based on x-chromosomal markers in the analysis; PLINK standard parameters were used). In addition, evidence for genotyping errors at the sample / marker level was evaluated by searching for an excess of “Mendelian inconsistencies” in the data (as indicated by substantial deviation from the empirical distributions in both measurements). Multidimensional scaling (MDS) analyses were performed on singleton OCD cases and unselected controls, as implemented into Plink [26]. Samples were removed when they significantly deviated in the first two MDS dimensions (> 4 standard deviation from the mean). Although inherently robust against population stratification, substantial heterogeneity in the population-based study cohorts would have lead to decreased power for our analytical approach (see below). More details on the QC process are provided in the supplementary materials and methods section. Use of the described filters resulted in a final analysis dataset comprising 5 061 individuals from the OCGAS samples and additional controls, with available genotypes for 549 123 autosomal markers.

Statistical analyses

For the single-marker association analyses, we used a recently published method for combining family-based and population-based data. More precisely, family-based association testing of within-family information is combined with population-based analyses of between-family information and the association analyses of unrelated study subjects [27]. This hybrid approach is inherently robust to population stratification and potentially increases statistical power compared to a classic meta-analysis design [27]. As part of our analytical pipeline we therefore computed P-values for all autosomal markers and both the within and the between family information using PBAT [28]. The two P-values were subsequently combined using a weighted Z-score statistic as implemented via METAL [29].

Gene-based statistics were derived as fixed Z scores, as implemented in FORGE [30]. More information on the approach is provided elsewhere [31]. Here we used the single-marker results of our combined OCGAS GWAS. Information about the correlation pattern in the data was provided through usage of HapMap phase 2 samples [32]. A maximum of 1 000 000 permutations were used per gene (adaptive approach) and analyses included an additional +/− 20kb sequence information based on positions obtained from ENSEMBL v70 [33]. We used the gene-based results in two ways: First we used them to agnostically search for genes that are associated with OCD. Second, we used them to follow-up on gene-set based results from the IOCDF-GC study, which reported an enrichment of association signals for two gene sets that comprised high confidence targets of two miRNA families [21]. In addition, we used information from a global interactome for Homo sapiens in order to identify high confidence interactors of DLGAP1 and GRIK2 (confidence threshold > 0.95) [34]. No sub-network reduction was applied and only genes representing first neighbors of DLGAP1 and/or GRIK2 in the global interactome were considered for this analysis. We used gene-based results in order to identify an enrichment of association with OCD in interaction partners of DLGAP1 and the ionotropic kainate 2 glutamate receptor (GRIK2). This analysis was motivated by the assumption that biologically closely related genes of these previously described OCD risk genes would make reasonable candidates for hypothesis driven downstream analyses. Based on an assumed similarity between these genes with regards to an involvement in common biological processes and provided that, on a broader level, these biological processes themselves are associated with the phenotype under study it seems reasonable to hypothesize that focusing on biologically closely related genes (see above) helps to identify new disease genes (through reduction of some of the multiple testing burden in standard GWAS). More information on the usage of the interactome data and the visualization of the resulting networks are found in the supplementary material and methods. In brief, we would like to emphasize that our analyses represent a simple, hypothesis driven approach for inclusion of the human interactome data, rather than an exhaustive network science approach (making use among others of topological features of the interactome other then the status of direct interaction). While the latter has been successfully used with neuropsychiatric traits in the past, we felt that this evolving field still needs some further improvements and therefore decided to refrain from these kinds of analyses.

Results

Single SNP association

Analysis of the 549 123 autosomal markers revealed a nominally significant result at the level of α=0.05 for 27 283 markers. Among these, 456 markers reached a P-value of less than 0.001 and 54 markers were identified with a P<0.0001. As expected due to the stringent study design, no evidence for population stratification was observed (see also Supplementary Figure S1). Figure 1 shows the ‘Manhattan’ plot of the association results across the autosomal genome and Table 1 reports the association between single SNPs and OCD with P<1×10−5. No marker tested for association met the standard of genome-wide significance. The smallest P-value for our study was observed with a marker on chromosome 9 (rs4401971) at a P-value of 4.13×10−7. This SNP is 1.28 Mb from the 5′ end of the protein tyrosine phosphatase, receptor type, D gene (PTPRD, isoform 1; according to RefSeq). The second most significantly associated SNP (rs6876547, P=1.76×10−6) is located in a region of cadherin clusters and it is of note, that a second, independent SNP (rs6452234, P=1.13×10−5, r2<0.2, distance ~650kb) is located in the same region. The nearest flanking protein coding genes are the cadherin 9, type 2 (CDH9, 1 308 kb) and the cadherin 10, type 2 (CDH10, 927 kb) genes. A ‘regional association plot’ of these two top regions is provided in Figure 2.

Figure 1

Manhattan plot for OCGAS GWAS

Shown are the result for the hybrid analysis of the within and between family component of the OCGAS GWAS. A thin blue line indicates level of suggestive evidence for association (1 × 10−5) and a thin red line indicates genome-wide significance (5 × 10−8).

Table 1

Results for OCGAS GWAS

SNPs listed are strongest associated GWAS variants (P < 0.0001) in the hybrid analysis of the within and between family component. Only one SNP (“lead SNP” = SNP with lowest p-value) per clump is given. Clumps are defined using PLINK (--clump-p1 0.0001, --clump-p2 0.05, --clump-r2 0.20, --clump-kb 5000, LD based on additional controls that have also been used for the hybrid analysis). Headings are as follows: MARKERINFO (# - Rank place for clump for which this SNP is the lead SNP; CHR – Chromosome; SNP – marker name for lead SNP in clump; BP – position; MA – minor allele in additional controls; MAF – frequency of this SNP in additional controls;), OCGAS (INFO – number of informative families in the within family component analysis; DIR – direction of effect with respect to the MA; P – p-value for the lead SNP in the clump; N – number of additional SNPs in the clump meeting the aforementioned criteria), and ANNOTATION (GENE – gene in which lead SNP is physically located; LEFT GENE and RIGHT GENE – nearest flanking genes).

MARKERINFO						OCGAS				ANNOTATION
#	CHR	SNP	BP	MA	MAF	INFO	DIR	P	N	GENE	LEFT GENE	RIGHT GENE
1	9	rs4401971	11890045	A	0.41	619	− −	4.13 × 10−7	31	NA	LOC646114	LOC100049717
2	5	rs6876547	25572301	G	0.19	471	++	1.76 × 10−6	6	NA	CDH10	MSNL1
3	7	rs1343795	54313171	A	0.08	229	− −	9.69 × 10−6	6	NA	FLJ45974	LOC222005
4	5	rs6452234	24922789	A	0.21	429	− −	1.13 × 10−5	10	NA	CDH10	MSNL1
5	14	rs1014951	24039559	T	0.1	311	++	1.49 × 10−5	6	JPH4	LOC100131731	DHRS2
6	8	rs7462051	80266553	A	0.18	462	++	1.69 × 10−5	21	NA	IL7	LOC100128963
7	16	rs1544352	19713882	C	0.16	445	++	1.94 × 10−5	19	NA	C16orf62	C16orf88
8	4	rs1532154	157279663	C	0.21	427	− −	2.22 × 10−5	3	NA	FTHP2	hCG_1814936
9	11	rs509876	79811268	G	0.13	406	++	2.29 × 10−5	4	NA	LOC646112	LOC729790
10	20	rs10392	37550935	A	0.19	467	− −	2.87 × 10−5	18	PPP1R16	ACTR5	FAM83D
11	9	rs2821204	11683901	G	0.2	492	++	2.89 × 10−5	29	NA	LOC646114	LOC100049717
12	10	rs1088258	97109019	A	0.33	577	− −	3.56 × 10−5	22	SORBS1	PDLIM1	LOC643981
13	8	rs2278144	25634077	T	0.13	382	− −	3.94 × 10−5	2	NA	CDCA2	EBF2
14	15	rs8026755	37247722	A	0.33	594	− −	4.03 × 10−5	8	MEIS2	LOC145845	LOC390576
15	3	rs838209	176352763	C	0.34	616	− −	4.14 × 10−5	5	NA	LOC730168	TBL1XR1
16	5	rs1773562	168653502	C	0.1	312	++	5.03 × 10−5	4	SLIT3	LOC728095	CCDC99
17	8	rs7003102	126303493	G	0.14	416	++	5.14 × 10−5	6	NSMCE2	KIAA0196	TRIB1
18	8	rs1254734	26833198	A	0.2	486	− −	5.38 × 10−5	11	NA	LOC100127897	LOC100132229
19	2	rs7593878	44358504	C	0.15	378	− −	5.44 × 10−5	6	NA	LRPPRC	PPM1B
20	2	rs1686740	182189405	T	0.11	345	++	6.10 × 10−5	5	NA	LOC729026	LOC100127923
21	13	rs9541148	35095884	C	0.37	624	− −	6.73 × 10−5	5	NA	LOC100130499	LOC100129452
22	18	rs1671253	50040561	T	0.29	555	− −	7.04 × 10−5	19	DCC	LOC100132995	LOC100133176
23	6	rs973714	133578733	A	0.34	615	++	7.44 × 10−5	12	EYA4	LOC285735	TCF21
24	10	rs3902042	17065067	T	0.39	641	++	7.71 × 10−5	12	CUBN	RSU1	TRDMT1
25	20	rs8120171	50712059	T	0.19	414	− −	7.85 × 10−5	1	ZFP64	SALL4	ERP28P
26	4	rs1685101	75216365	A	0.24	467	− −	7.89 × 10−5	15	NA	EPGN	EREG
28	3	rs9845643	179531062	G	0.2	457	− −	8.58 × 10−5	7	PEX5L	USP13	LOC647249
29	9	rs1126590	92077132	C	0.39	634	− −	8.60 × 10−5	12	SEMA4D	PA2G4P6	LOC100128670
30	20	rs1699787	16533162	C	0.23	518	++	8.68 × 10−5	12	KIF16B	LOC100131642	RPLP0P1
31	9	rs2183738	71211200	G	0.14	366	++	8.73 × 10−5	8	NA	C9orf71	LOC347097
32	1	rs1209698	185201061	C	0.34	612	++	8.87 × 10−5	32	C1orf26	LOC100129295	IVNS1ABP
33	8	rs7005206	130620813	A	0.15	377	− −	9.13 × 10−5	7	NA	LOC100129525	LOC100130376

Figure 2

Regional association plots for top regions in OCGAS GWAS

The most associated marker in the region (see table 2 and S2, purple dot) is centered in a genomic window of 5 Mb (hg19). P-values for the OCGAS GWAS are given. The linkage disequilibrium (LD) strength (r2; data from the 1000 genomes project European samples) between the sentinel single nucleotide polymorphism and its flanking markers is demonstrated by the coloring of the dots for the neighboring markers (ranging from red = high to blue = low). The recombination rate (cM/Mb; second y axis) is plotted in blue. Plots are given for the a) chromosome 9 region harboring PTPRD (clump 1 in table 2), b) the chromosome 5 region (clump 2 in table 2), and c) the chromosome 16 region harboring the two genes that were identified in the gene-based analysis, IQCK and C16orf88.

We used the single-maker association signals to follow up the results from the recent IOCDF-GC GWAS [21]. For the markers listed in Table 1 of that study, we found evidence for association with a marker (rs2205748, same effect allele and direction of effect) on chromosome 6 near the genes GRIK2 (distance to gene 1.94 Mb) and the HECT domain and ankyrin repeat containing E3 ubiquitin protein ligase 1 gene (HACE1, 713 kb) with an uncorrected P-value of 0.0493 in our study (P=8.52×10−6 in the original IOCDF-GC trio-case-control analysis). One additional region (harboring the gene for the cAMP-specific phosphodiesterase 4D, PDE4D) showed a nominal significant signal; however, the effect (based on the same effect allele) was in the opposite direction (compared with the IOCDF-GC finding). Though no other marker among the tested SNPs from the IOCDF-GC GWAS reached a significant result, it is of note that, for the region around DLGAP1 an independent SNP in our analyses (rs3866988) yielded a P-value of 2.67×10−4 and thus indicating this region as a susceptibility factor for OCD (r2=0.023, D′=0.316 between rs3866988 and rs11081062/ rs11663827, based on 1 000 genomes data). It is also of note, that 12 out of 15 markers (counting the DLGAP1 markers in perfect LD only once) show the same direction of effects for the effect alleles in both the IOCDF-GC and OCGAS analyses (binominal test with probability 0.5 = sign test: P=0.0176). Detailed results for the follow-up of the IOCDF-GC results are found in Table 2.

Table 2

Follow up of IOCDF-GC GWAS results

SNPs listed are strongest associated GWAS variants in trio, case-control and combined trio-case-control samples reported in table 1 of Stewart et al. 2013. Only the results from the trio-case-control sample are listed in the table (IOCDF –GC: P) along the results from the current GWAS (OCGAS: P). Directions of effects (IOCDF-GC: DIR and OCGAS: DIR) are given based on the A1 alleles in the IOCF analyses (IOCDF-GC: A1). For the IOCDF-GC GWAS directions are given for all study samples (EU, AJ, SA, trios). In case of the OCGAS GWAS the directions refer to the combined results (within and between family component). For SNPs not present in the OCGAS results a proxy SNP was used that was identified using the SNAP tool (OCGAS: PROXY; max distance to query SNP 500kb, rsq min 0.8, reference data: 1000 genomes project pilot 1 release).

MARKERINFO			IOCDF-GC			OCGAS			GENES (kb)
SNP	CHR	BP (HG19)	P	A1/A2	DIR	P	DIR	PROXY (rsq)
rs11898020	2	144282078	2.65E-04	A/G	+−+ −	0.6357	−	N/A	ARHGAP15 (intronic)
rs10165908	2	158315629	0.0169	C/T	+++−	0.5816	+	rs2198581 (0.83)	CYTIP (15), ACVR1C (68)
rs1838733	5	58533392	3.82E-05	T/C		0.0482	+	rs6858946 (0.95)	PDE4D (intronic)
rs26728	5	106946056	1.01E-04	T/C	−++ −	0.3289	−	rs154191 (0.85)	EFNA5 (intronic)
rs4868342	5	173504522	3.20E-05	C/T	++++	0.8425	−	N/A	HMP19 (intronic)
rs9499708	6	104445367	2.96E-06	C/T	− −+ −	0.1287	−	N/A	GRIK2 (1927), HACE1 (731)
rs2205748	6	104462555	8.52E-06	G/A	++++	0.0493	+	N/A	GRIK2 (1944), HACE1 (713)
rs182320	6	130073291	2.25E-05	C/T	++++	0.5801	+	N/A	ARHGAP18 (42), C6orf191 (79)
rs6531002	8	12722703	0.0067	T/C	−++ −	0.4897	−	rs7823534 (0.96)	LONRF1 (110), KIAA1456 (80)
rs11611761	12	33025612	0.115	A/C	− −+ −	0.2676	−	rs17543624 (0.87)	PKP2 (intronic)
rs297941	12	50319086	4.99E-07	G/A		0.8846	−	N/A	FAIM2 (21), AQP2 (25)
rs9652236	13	72688774	5.14E-06	T/G	++++	0.4017	+	N/A	DACH1 (247), MZT1 (594)
rs11081062	18	3662879	2.92E-05	T/C	+0++	0.6274	+	N/A	DLGAP1 (intronic)
rs11663827	18	3663631	2.31E-05	A/G	++++	0.6844	+	N/A	DLGAP1 (intronic)
rs485186	19	49207206	9.94E-06	G/A	++++	0.5067	+	N/A	FUT2 (coding-synon)
rs6131295	20	11996267	3.63E-05	G/A	+	0.5663	−	N/A	BTBD3 (89), SPTLC3 (993)

Gene-based findings

Gene-based analyses for 21 567 genes (protein-coding genes and miRNAs) resulted in an experiment-wise significant result for two genes, the IQ motif containing K gene (IQCK) and the chromosome 16 open reading frame 88 gene (C16orf88; Pcorr<0.0215; P<1×10−6). Both genes are located in the same chromosomal region (16p12.3) and share their leading SNP, i.e. the SNP with the lowest P-value that is annotated to the gene (rs1544352; P=1.94×10−5). The only other gene that shows a P-value < 1×10−4 is the orofacial cleft candidate1 gene (OFCC1; P=6.29×10−5). A more detailed list of association results for the gene-based analyzes are given in Table 3 and a ‘regional association plot’ for the OFCC1 locus is provided in Figure 2.

Table 3

Gene-based results

Genes listed are strongest associated protein-coding genes in the FORGE analysis (P < 0.001). Annotation information is taken from ENSEMBL v70 (GENEINFO: CHR, START, STOP), gene names (GENEINFO: GENE) are following HGNC nomenclature. The results (ASSOCIATION: P – gene-based p-value; P_CORR – experiment-wise corrected gene-based p-value; N – number of SNPs annotated to the gene; minP – p-value of SNP with lowest p-value that was annotated to the gene) are based on the Z_FIX statistic as implemented into FORGE. Annotation used a +/− 20kb (from GENEINFO: START and STOP) flanking region for each gene under study.

GENEINFO				ASSOCIATION
CHR	START	STOP	GENE	P	P_CORR	N	minP
16	19714902	19729557	C16orf88	9.99E-07	0.019	7	1.94E-05
16	19727778	19868907	IQCK	9.99E-07	0.019	20	1.94E-05
6	9596343	10211841	OFCC1	6.29E-05	1	152	2.86E-04
14	24037244	24048009	JPH4	1.16E-04	1	5	1.49E-05
15	79252289	79383115	RASGRF1	1.18E-04	1	53	6.25E-04
19	2100988	2164464	AP3D1	1.96E-04	1	17	1.93E-03
9	116148597	116163613	ALAD	2.56E-04	1	23	9.83E-03
1	185087220	185126204	TRMT1L	2.93E-04	1	9	1.73E-04
4	75230860	75254468	EREG	3.06E-04	1	19	7.89E-05
20	37590942	37668366	DHX35	3.55E-04	1	7	6.22E-05
6	88117701	88221352	C6orf165	3.55E-04	1	15	1.84E-03
15	41186628	41196173	VPS18	4.36E-04	1	6	1.02E-03
1	185126212	185260897	SWT1	5.00E-04	1	22	8.87E-05
6	88180341	88222054	SLC35A1	6.11E-04	1	10	1.84E-03
1	244571796	244615436	ADSS	6.30E-04	1	14	2.49E-04
9	97418353	97480105	C9orf118	6.32E-04	1	16	4.96E-04
9	71939488	72007371	FAM189A2	6.47E-04	1	18	9.37E-04
3	149191761	149221068	TM4SF4	6.67E-04	1	23	6.00E-04
10	97071528	97321171	SORBS1	6.84E-04	1	86	1.34E-04
16	19566562	19718115	C16orf62	7.54E-04	1	37	1.94E-05
19	19287712	19303425	MEF2BNB	8.09E-04	1	11	8.85E-04
10	90562654	90580303	LIPM	8.46E-04	1	21	2.73E-03
14	24025216	24029480	THTPA	8.90E-04	1	8	1.49E-05
1	204042243	204096863	SOX13	9.17E-04	1	37	1.57E-03
1	231359509	231376933	C1orf131	9.33E-04	1	15	4.05E-03
6	133561736	133853258	EYA4	1.00E-03	1	84	7.44E-05

Query of the interactome [34] for high confidence interactors of DLGAP1 and GRIK2 identified 169 interactors of DLGAP1 and 161 interactors for GRIK2 (246 interactors for both of them together). Gene-based analyses revealed that 21 (out of 246) genes showed a nominal significant P-value (one-sided P=0.075). This included 16 genes for DLGAP1 (one-sided P=0.069) and 14 genes for GRIK2 (one-sided P=0.135). It is of note that GRIK2 itself was identified as interactor of DLGAP1 and showed a nominal significant gene-based P-value of 0.03. Among the list of other nominal significant genes are the neuronal differentiation 6 gene (NEUROD6, P=0.010), the synaptic vesicle glycoprotein 2A gene (SV2A, P=0.026), the ionotropic, AMPA 4 glutamate receptor gene (GRIA4, P=0.039), and the solute carrier family 1 (glial high affinity glutamate transporter), member 2 gene (SLC1A2, P=0.035). A full list of these results is given in Supplementary Table S4.

We also attempted to replicate a reported enrichment of association signals for two gene sets that comprised high confidence targets of two miRNA families (see the Stewart et al OCD GWAS). Analyses of the miRNA families of predicted miR-130ac/301ab/301b/301b-3p/454/721/4295/3666 (miRNA set 1) and miRNA-219-5p/508/508-3p/4782-3p (miRNA set 2) targets showed nominal significance for 3 out of 45 high confidence targets in miRNA set 1 and 14 out of 145 high confidence targets in miRNA set 2. The number of nominal significant genes for both sets are higher then expected by chance; however, for both sets, no significant enrichment for OCD associated genes was found in either set (one-sided P=0.323 and one-sided P=0.087, respectively). It is of note that one of the high confident targets in miRNA set 2 (the epiregulin gene, EREG) did demonstrate an experiment-wise significance level (correcting for the number of tests performed in miRNA set2; P=3.06×10−4; Pcorr=0.044).

Discussion

Here we report the results of a GWAS including a sample of 1 406 patients with OCD, that was predominantly family-based, but which also included a case-control subsample to increase power (resulting in a total sample of 5 061 individuals). The study identified interesting candidate genes for OCD, but failed to detect any genome-wide significant findings. This is similar to what has been observed for other psychiatric phenotypes such as schizophrenia, for which genome-wide significant findings were only achieved after either starting with more samples in the discovery step [35, 36] or making use of large follow-up samples, that were orders of magnitude larger then the original discovery samples [37, 38]. We suspect that, with additional samples, the findings of genetic studies for OCD will be more robust, and our currently suggestive findings may reach genome-wide significance.

Of particular interest is the signal near the PTPRD gene (rs4401971), the most significant found in this study (P=4.13×10−7). PTPRD is a member of the receptor protein tyrosine phosphatase family, which comprises transmembrane signaling molecules that regulate a variety of cellular processes including cell growth and differentiation [39]. Pre-synaptic PTPRD promotes the differentiation of glutamatergic synapses [40-43] and interacts with Slit and NTRK-like family member 3 (SLITRK3) which acts as a postsynaptic adhesion molecule. Together both proteins selectively regulate inhibitory GABAergic synapse development [44]. This is interesting because molecules in the same family, SLITRK5 and SLITRK1, have been shown to be associated with TS and OCD, the former in a mouse model and the latter in a TS genetic study [45, 46].

Mice deficient in PTPRD show impairment in learning and memory tasks, especially spatial learning, and exhibit enhanced long-term potentiation (LTP), a form of activity dependent plasticity [47]. This is relevant to OCD, in that memory deficits have been reported for this condition [48]. PTPRD also has been shown to be associated with restless leg syndrome (RLS) [49] and hemizygous deletions were detected in four unrelated attention deficit hyperactivity disorder (ADHD) probands [50].

Our second strongest association finding is located in a region of cadherin clusters. It is known that CDH9 has been reported associated with OCD [51]. Furthermore, six SNPs between cadherin 10 (CDH10) and cadherin 9 (CDH9) were significantly (as well as in a replication sample) associated with autism spectrum disorders (ASD) [52] and was a top (although not significant) association signal in anorexia nervosa [53].

Analysis of previously identified GWAS hit regions for the IOCDF-GC study [21] revealed that 12 out of the 15 strongest signals in that study showed associations with the same direction of effects (based on the same effect alleles) in both the IOCDF-GC and OCGAS analyses (sign test P=0.0176). This observation supports the hypothesis that future collaborative studies (meta- and mega-analyses) will be able to identify OCD genes at the level of genome-wide significance. Support for association with OCD was particularly found in two previously identified genetic regions (see above). The most interesting region harbors the DLGAP1 gene. In the IOCDF-GC study signals in this gene were the top signals from the case-control analysis. Although the specific SNP from the former study was not found significant in this study, there was a prominent signal for an independent nearby marker (P=2.67×10−4). This observation potentially serves as independent evidence for association of markers in this region. DLGAP1 is a member of the neuronal postsynaptic density complex and is in the same family as the DLGAP3 gene, which has been shown in a convincing way to be responsible for OCD-like behaviors in DLGAP3 knockout mice models [54]. A marker in the second genetic region (harboring GRIK2), also a top signal in the IOCDF-GC GWAS showed nominal significance in our study (P=0.045). Although failing to reach experiment-wise significance it is of note that DLGAP1 and GRIK2 interact (FUNCOUP [34]). Following up these signals in a gene-set analysis for high confidence interaction partners of DLGAP1 and GRIK2 (identified using FUNCOUP [34]) showed a trend for association and pointed to a potential role of a set of DLGAP1 and GRIK2 interactors in the etiology of OCD, involving genes such as NEUROD6, SV2A, GRIA4, and SLC1A2. Interestingly PTPRD was part of this gene set. It is of note that for both DLGAP1 and GRIK2 we identified more then 160 interactors each, thus indicating they are highly connected nodes (hubs) in the interactome. Earlier findings have suggested that only genes that are essential in (early) development (“essential genes”) tend to encode hub proteins, while the vast majority of disease genes are nonessential and show no tendency to encode hub proteins [55]. This might point to DLGAP1 and GRIK2 as being essential genes in early (neuro-)development. However, more studies are warranted to confirm this hypothesis.

Gene-based association analyses identify significant signals in two genes, IQCK and C16orf88; however, this finding is not amenable to interpretation at this time. The next strongest gene-based signal in OFCC1 is intriguing because a variant within the exon of this gene has been reported in seven affected individuals of family segregating OCD, TS, and ADHD [56].

This study employed careful and comprehensive phenotyping; MD and PhD-trained psychologists conducted all assessments and recruitment sites used the same assessment approach and all cases were reviewed at a single site (JHU) to ensure diagnostic replicability. Given the likely etiologic heterogeneity of OCD, the comprehensive phenotyping provides opportunities for future studies (through subtyping based e.g. on factor-analytical approaches).

A limitation of this study is that the power to detect significant association signals for small effect sizes is limited by the sample size. This has been clearly demonstrated in many GWAS efforts for other neuropsychiatric conditions (e.g., schizophrenia) in which studies of comparable size have failed to identify significant association signals. Yet, when sample sizes have increased (in some case 10 or more fold), then multiple significant signals were identified. Power was improved in this study by using a hybrid analytic approach that included both family and case-control samples in the analysis. In conclusion, while this OCD GWAS study did not identify a study-wide significant association finding, several of the strongest findings are particularly interesting. There are both plausible biologic hypotheses and prior genetic evidence, for either OCD or related conditions, for the two most significant association findings, PTPRD and CDH9/CDH10. The finding that the signals in this GWAS concur with respect to allele and direction significantly with the top signals in a second reported OCD GWAS, suggests that similar genetic underpinnings of OCD are identified in both studies. Moreover, while not a replication in the full sense, findings for the genes DLGAP1 and GRIK2 suggest that they are both relevant candidates genes for OCD. The suggestive findings in this study await replication in larger samples.