De novo mutations identified by whole-genome sequencing implicate chromatin modifications in obsessive-compulsive disorder.

Obsessive-compulsive disorder (OCD) is a chronic anxiety disorder with a substantial genetic basis and a broadly undiscovered etiology. Recent studies of de novo mutation (DNM) exome-sequencing studies for OCD have reinforced the hypothesis that rare variation contributes to the risk. We performed, to our knowledge, the first whole-genome sequencing on 53 parent-offspring families with offspring affected with OCD to investigate all rare de novo variants and insertions/deletions. We observed higher mutation rates in promoter-anchored chromatin loops (empirical P = 0.0015) and regions with high frequencies of histone marks (empirical P = 0.0001). Mutations affecting coding regions were significantly enriched within coexpression modules of genes involved in chromatin modification during human brain development. Four genes—SETD5, KDM3B, ASXL3, and FBL—had strong aggregated evidence and functionally converged on transcription’s epigenetic regulation, suggesting an important OCD risk mechanism. Our data characterized different genome-wide DNMs and highlighted the contribution of chromatin modification in the etiology of OCD.

INTRODUCTION

Obsessive-compulsive disorder (OCD) is a severe, disabling neuropsychiatric disorder characterized by repetitive, compulsive thoughts, and/or compulsive movements or rituals, with an estimated lifetime prevalence of 1 to 3% (–). OCD is typically diagnosed on the basis of observed behaviors, duration of symptoms, and impairment of function, rather than based on a biological understanding of the disease etiology. This has notably hindered progress in developing more precise diagnoses and better therapeutics to improve patient outcomes. Epidemiological studies have demonstrated the important contributions of genetic factors to the etiology of OCD. On the basis of concordance rates in monozygotic and dizygotic twin studies, the heritability of OCD is approximately 50% (). Although some putative chromosomal regions have been highlighted by traditional genetic study (), no causative genes for OCD have been confidently identified ().

More recently, genome-wide association studies (GWASs) (, ) have identified single-nucleotide polymorphisms (SNPs) associated with OCD with roles in glutamate signaling and synaptic functions. However, the SNP-based heritability of OCD is estimated to be 0.22 (), suggesting that the remaining contribution to the disease risk cannot be explained by weakly associated common SNPs (). The apparent contradiction between the profound heritability of OCD based on epidemiological studies and the limited polygenic risks estimated by GWAS, the so-called “missing heritability,” has given rise to the hypothesis of a “de novo paradigm,” which suggests the potential role of rare variants in the development of OCD. More recently, two whole-exome sequencing (WES) studies of parent-offspring OCD trios by Cappi et al. (, ) have shown compelling evidence for the role of de novo mutations (DNMs) in OCD, which adds to the accumulating proof for the hypothesis that de novo variants have a significant role in the genetic architecture of psychiatric disorders (). The two OCD exome studies (, ) (the number of trio = 20 and 184, respectively) estimated that about 22% of patients with OCD carried a damaging coding DNM that could estimate the occurrence of the disease. In the meantime, they highlighted genes, such as CHD8 and SCUBE1, which have a potential role in the pathology of OCD.

However, most of the large-scale investigations on OCD susceptibility have been focused on attributing disease risk to the protein-coding regions of the genome by targeted-sequencing () or exome-sequencing approaches (, ). Apparent highly penetrant pathogenic variants in intergenic, noncoding RNA (ncRNA), and large structural variant (SV) regions (, ) are also known, suggesting that whole-genome sequencing (WGS), which detects more classes and sizes of mutations than WES, could be considered as the preferred genomic platform to study OCD. Furthermore, the application of WGS in other neuropsychiatric disorder studies, such as autism disorder spectrum (ASD) (–) and intellectual disability (ID) (), has already demonstrated its unparalleled power for understanding the overall burden of DNMs at the genome scale and their contributions to diseases.

Thus, following the previous exome studies in OCD, we applied WGS and bioinformatics analysis to investigate a cohort of OCD parent-child trios to identify genome-scale de novo variants, including single-nucleotide variants (SNVs), small insertions and deletions (indels), and SVs. In 53 parent-offspring families, we found strong evidence for a high occurrence of mutation at promoter-interacting loops. We also observed that DNMs preferentially occurred on intolerant genes and affected genes regulating chromatin modification. Furthermore, we identified three high-confidence chromatin modifiers—SETD5 (SET-containing-domain 5), KDM3B (lysine demethylase 3B), and ASXL3 (transcriptional regulator 3)—as OCD candidate risk genes based on multiple lines of evidence. Last, our integrated analysis found that genes affected by OCD and Tourette disorder (TD) DNMs (–) had similar biological functions but affected distinct brain regions and cell-type expression patterns, which, in part, explains the complex comorbidity relationship between the two neurodevelopmental disorders.

RESULTS

Global identification of de novo variants and their distributions in OCD

We performed WGS of 160 DNA samples from a cohort of 53 families with 54 OCD probands and their unaffected parents without major psychoses (Materials and Methods and table S1). After the quality control, 53 probands from 51 trios and one quad, as well as their unaffected parents, remained for subsequent analysis with one trio family removed (the total number of participants is 157). On average, 89.1% of genomic regions were covered by 20 or more reads at the individual level. The average sequencing depth was 29.92× (fig. S1). Among 53 parent-offspring families, we identified 4143 de novo SNVs and short indels initially. After the Sanger-sequencing validation with ~96% positive rate (see Materials and Methods and table S2), we retained a list of 4062 DNMs (table S2) with a mutation rate of 1.34 × 10−8 (95% confidence interval: 1.26 × 10−8 to 1.41 × 10−8) per base pair. The number of DNMs in each proband ranged from 51 to 117 in the genome (zero to four in the exome). Thirty-five probands carried one or more exonic DNMs. The per-individual number of DNMs was 0.93 for SNVs and 0.13 for indels across coding regions, 62.5 for SNVs, and 14.7 for indels across the genome, similar to those reported in other family-based WES and WGS studies of psychiatric diseases (–).

Next, we calculated the relative mutation rates (number of DNMs of a particular type divided by the total number of DNMs in the study) globally for five general mutation types: exonic, intronic, untranslated region (UTR), ncRNA, and intergenic. Similar to the results reported previously in a large WGS study in ASD (), differences between OCD probands and Simons Simplex Collection (SSC) controls () in most types did not reach statistical significance after the multiple testing correction except for ncRNA (P = 0.02). Furthermore, we defined a damaging DNM category as one causing de novo loss-of-function (LoF; frameshift indel, stop-gain, or splice-site mutation) or pathogenic missense [predicted by SIFT () and PolyPhen-2 () with a damaging effect, and the PHRED-scaled score of combined annotation dependent depletion (CADD) () ≥ 20]. We then evaluated the rate of damaging mutations in OCD against SSC controls () by calculating the relative exonic mutation rates against the synonymous mutation rates, since the usage of synonymous mutations as an internal control should be resistant to potential artifacts caused by the comparison of data from different studies. We observed a significantly higher mutation rate resulting in damaging DNMs in OCD [P = 0.02, odds ratio (OR) = 2.32]. Moreover, no significant difference was observed in guanine-cytosine content (GC content) (Wilcoxon’s test P = 0.1) or gene length (Wilcoxon’s test P = 0.9) between DNM genes in cases and those in controls, indicating that the observed significance in the rate of OCD damaging mutation was not introduced by biases, such as GC content or gene length.

DNMs in OCD are enriched in promoter-interacting loops and regions with high frequencies of histone marks

Next, we investigated whether OCD mutations occurred more frequently in functional genomic regions by integrating genomic annotations from multiple different data sources, including PsychEncode (), ENSEMBL (), and Roadmap epigenomics (). We observed no global enrichment of DNMs directly hitting promoters, enhancers, or full-length genes with noncoding regions after correcting for multiple testing, which was also similar to a previous large WGS study of schizophrenia (SCZ) trios (). However, we observed that our OCD mutations were significantly enriched in the promoter-anchored chromatin loops (empirical Pcorrected < 0.001 and OR = 1.23; Fig. 1A), which are promoter-interacting distal regulatory elements. To rule out the potential bias by study design and sequencing technical issues, we collected two additional control datasets (, ) and found that the control DNM distribution in promoter-anchored chromatin loops was similar in different studies (fig. S2; Cochran test P = 0.91). We compared our OCD data to the three control data, estimated SE by subsampling test and leave-one-chromosome-out (LOCO) (Materials and Methods), and obtained consistent enrichment with no cross-study heterogeneity (subsampling: meta-analysis OR = 1.25, P < 2.2 × 10−16, Cochran test P = 0.94; LOCO: meta-analysis OR = 1.25, P < 2.2 × 10−16, Cochran test P = 0.15).

Fig. 1.

Rate of DNMs across genomic regions.

(A and B) Distribution of the per-individual rate of mutations within (A) different functional genomic regions and (B) different chromatin states. The P values represent the significance of differences in controls or patients with other diseases and were corrected for multiple testing using a false discovery rate. The X axes corresponded to ORs estimated by the Fisher’s tests, and error bars represent 95% confidence intervals. (C) Schematic and Integrative Genomics Viewer plots are showing two dSVs validated by quantitative polymerase chain reaction (qPCR). (D) Comparison of intolerance of genes affected by the same four types of DNMs. SSC control, Simons Simplex Collection siblings; pLI, the probability that a gene is intolerant to a LoF mutation; dSV, de novo SV.

Chromatin states, which are the combination of multiple epigenomic marks in a spatial context, can capture information describing different classes of genomic elements such as promoters, enhancers, transcribed, repressed, and repetitive regions. Thus, we extended our analysis of DNM burden to investigate their occurrence at epigenomic marks in cases versus controls by mapping mutations to the 15 core chromatin states (). We observed a significantly higher mutation occurrence within regions of zinc finger genes and repeats (chromatin state 8; empirical Pcorrected < 0.0001 and OR = 9.65; Fig. 1B), which is characterized by high frequencies of histone marks H3K9me3, H4K20me3, and H3K36me3 and relatively low frequencies of other marks. Similar to above, the DNM distributions among chromatin states were highly concordant between different control studies [Pearson correlation coefficient (R) > 0.96; fig. S2], and the enrichment result was consistent across these studies (subsampling: meta-analysis OR = 7.24, P = 1.23 × 10−10, Cochran test P = 0.87; LOCO: meta-analysis OR = 7.46, P < 2.2 × 10−16, Cochran test P = 0.07).

Small de novo structural variation affecting FBL gene in patient with OCD

Aside from identifying point mutations and indel, we also explored whether de novo SVs (dSVs) had a role in OCD. After machine learning detection, manual inspection, and population frequency filtration (Materials and Methods), we identified seven potential dSVs not presented in the GnomAD population. Two of them were quantitative polymerase chain reaction (qPCR) verified as true de novo variants on the basis of the sample availability. One of them was identified as a heterozygous de novo deletion (chr19:40331154–40331255) of the second intron of gene FBL on a male patient (ID: WOC3_114_1_PT) with early-onset OCD (table S3). This dSV skipped most of the second intron and two bases of the 5′ of exon3 and inserted one base to the 3′ of exon2, leading to a frameshift-like alteration on FBL (Fig. 1C). FBL encodes Fibrillarin, a nuclear protein that methylates glutamine-104 of histone H2A (H2AQ104me) and modulates chromatin structure and activity (). The other confirmed dSV was on another male patient (ID: WOC4_34), also with early-onset OCD. The dSV was a heterozygous deletion of a 703–base pair (bp) fragment (chr12:64800102–64800404), which affected the part of the first intron of XPOT (Fig. 1C and table S3) and was predicted to have little effect on the gene-based function annotations from ANNOVAR ().

De novo coding mutations preferentially hit intolerant genes and affected genes regulating chromatin modification

Next, we explored the functional impact of the identified DNMs. Since genes differ in their tolerance to mutations, as measured by the probability of being LoF intolerant (pLI) (), we then investigated whether this intolerance differed across genes affected by different mutation classes (Fig. 1D). We observed that OCD-associated genes were significantly more mutation-intolerant than SSC controls with regard to coding DNMs (P = 0.003), damaging DNMs (P = 0.03), or missense DNMs (P = 0.002), but not for synonymous DNMs (P = 0.25).

Moreover, we reasoned that transcriptomic data from the developmental human brain would improve our understanding of OCD pathophysiology, as the cerebral cortex has been consistently implicated in OCD (). We examined the expression data from the Brainspan () and constructed prenatal and postnatal coexpression networks by weighted gene coexpression network analysis (WGCNA) () to access three sets of gene lists: set 1 contained protein-altering DNMs (Table 1), set 2 contained genes with promoters hit by DNMs, and set 3 was a combination of set 1 and set 2.

Table 1.

List of all protein-altering DNMs identified in the present study.

VUS, the variant with uncertain significance. NA, not available.

Proband ID	Position	Mutation	Gene symbol	Mutation type	Protein change	SIFT	Polyphen2 HDIV	CADD	Classification
WOC5_160_1	18:31324188	GC > G	ASXL3	Frameshiftdeletion	p.F1460Lfs*5	NA	NA	35	Damaging
WOC3_114_1_PT	6:29912028	AGG > AG	HLA-A	Frameshiftdeletion	p.D251Tfs*45	NA	NA	32	Damaging
WOC5_137_1	1:109242401	CAAATAAA>CAAA	PRPF38B	Frameshiftdeletion	p.K324Efs*59	NA	NA	34	Damaging
WOC5_26	2:114482963	CA > C	SLC35F5	Frameshiftdeletion	p.L408Rfs*11	NA	NA	35	Damaging
WOC5_160_1	19:16611793	CT > CTT	C19orf44	Frameshiftinsertion	p.K65Efs*15	NA	NA	14.66	Damaging
WOC5_6	19:2939254	AGTGAGGGAATGACACCACCCTTACCCAAGGAGGCA>A	ZNF77	Splice-sitemutation(frameshiftdeletion)	p.A41Lfs*5	NA	NA	16.62	Damaging
WOC5_160_1	3:49713641	C > T	APEH	Nonsense	p.Q199X	NA	NA	35	Damaging
WOC4_28	13:32937534	T > G	BRCA2	Nonsense	p.L2732X	NA	NA	39	Damaging
WOC5_102_1	14:20779873	G > A	CCNB1IP1	Nonsense	p.R224X	NA	NA	37	Damaging
WOC4_34	4:47427852	C > G	GABRB1	Nonsense	p.Y414X	NA	NA	27	Damaging
WOC4_165_1	6:39286844	C > T	KCNK16	Nonsense	p.W93X	NA	NA	43	Damaging
WOC4_170_1	22:25320164	C > T	SGSM1	Nonsense	p.R1008X	NA	NA	46	Damaging
WOC84	1:118616533	G > A	SPAG17	Missense	p.R777C	D	D	29.5	Damaging
WOC4_170_1	X:51640699	C > T	MAGED1	Missense	p.R515C	D	D	28.9	Damaging
WOC5_157_1	5:137722015	A > G	KDM3B	Missense	p.E362G	D	D	28.6	Damaging
WOC5_096_1	17:19284412	A > G	MAPK7	Missense	p.Y158C	D	D	27.5	Damaging
WOC5_62	1:36479519	C > T	AGO3	Missense	p.R192W	D	D	27.4	Damaging
WOC5_178_1	3:142741802	C > G	U2SURP	Missense	p.P376A	D	D	25.9	Damaging
WOC4_130_1_WY	3:9476069	C > T	SETD5	Missense	p.R77C	D	D	24.6	Damaging
WOC5_093_1	17:18195985	G > C	TOP3A	Missense	p.P324A	D	D	25.2	Damaging
WOL3_046_1_ZXY	19:58438675	C > T	ZNF418	Missense	p.G207R	D	D	21.7	Damaging
WOC5_093_1	1:206648328	C > T	IKBKE	Missense	p.R32C	D	D	28.8	Damaging
WOC5_099_1	4:142949945	G > T	INPP4B	Missense	p.A922D	D	D	25.7	Damaging
WOC5_149_1	8:131072874	G > A	ASAP1	Missense	p.T1048M	D	D	24.5	Damaging
WOC5_17	2:71817395	C > T	DYSF	Missense	p.A1152V	D	P	27.9	VUS
WOC5_164_1	12:96617430	A > T	ELK3	Missense	p.N29I	D	P	25	VUS
WOC5_4F	13:76381720	G > T	LMO7	Missense	p.R201L	D	P	23.2	VUS
WOC5_097_1	3:38595773	C > T	SCN5A	Missense	p.V1586M	D	P	23.1	VUS
WOC5_130_1	20:33594249	T > C	TRPC4AP	Missense	p.T606A	T	P	25.9	VUS
WOC5_161_1	X:123519751	G > A	TENM1	Missense	p.T1944I	D	B	23.6	VUS
WOC5_155_1	11:62416110	A > G	INTS5	Missense	p.V481A	D	B	23.5	VUS
WOC5_078_1_X	15:41810235	C > T	RPAP1	Missense	p.S1314N	T	B	22	VUS
WOC5_60	16:68225447	A > G	NFATC3	Missense	p.T959A	D	B	21.7	VUS
WOC5_091_1	6:32123547	G > A	PPT2	Missense	p.M140I	D	B	20.8	VUS
WOC5_076_1_WSH	14:91928489	T > C	PPP4R3A	Missense	p.M451V	T	B	20.4	VUS
WOC5_15	17:40149144	T > C	DNAJC7	Missense	p.S38G	T	B	19.65	VUS
WOC5_4F	3:102157378	C > T	ZPLD1	Missense	p.P32L	T	B	18.89	VUS
WOC5_094_1	19:47207624	G > A	PRKD2	Missense	p.P74S	T	B	16.23	VUS
WOC5_138_1	19:36223863	C > T	KMT2B	Missense	p.P2138L	D	P	15.46	VUS
WOC5_116_1	19:57036758	G > A	ZNF471	Missense	p.S367N	D	B	15.06	VUS
WOC5_152_1	12:113759149	C > T	SLC8B1	Missense	p.R54H	T	B	14.72	VUS
WOC5_39	2:217012900	C > T	XRCC5	Missense	p.T524I	T	B	14.16	VUS
WOC5_39	9:125551863	A > C	OR5C1	Missense	p.I218L	T	B	12.96	VUS
WOC5_084_1_PS	6:151671185	T > A	AKAP12	Missense	p.D455E	T	B	6.24	VUS
WOC5_150_1	18:76753591	G > C	SALL3	Missense	p.V534L	T	B	0.008	VUS
WOC5_091_1	X:114468473	CCCGCCGCCGCCGCCGCC>CCCGCCGCCGCCGCC	LRCH2	Nonframeshiftdeletion	p.G44del	NA	NA	16.51	VUS

We first observed that module M16 in the prenatal network contained significantly more set1 genes than expected (two-sided Fisher’s exact test P = 0.00012). M16 was enriched in functions including histone modification, neuron projection morphogenesis, and synapse organization, according to Metascape (Fig. 2A and table S4). When exploring the connectivity of DNM genes in M16, we observed three DNM genes acting as hubs, i.e., ASXL3 [degree = 161, kME (a measure of intramodular connectivity) = 0.78], SETD5 (degree = 45, kME = 0.75), and KDM3B (degree = 246, kME = 0.87), in M16 at an average degree of 41. Next, we observed that module M11 in the postnatal network contained significantly more set1 genes than expected (P = 0.007). M11 was mainly enriched in epigenetic regulation–related functions, such as covalent chromatin modification, histone lysine methylation, and chromatin remodeling (Fig. 2B and table S4). When investigating the connectivity of OCD DNM genes within M11, we also found two OCD genes, SETD5 (degree = 26, kME = 0.92) and KDM3B (degree = 27, kME = 0.84), that act as hubs in M11 with an average degree of 9. However, we did not observe similar results for set2 or set3 genes in both networks.

Fig. 2.

The biological significance of coding DNMs.

(A and B) Enrichment of DNM-affected genes within coexpression modules in (A) prenatal and (B) postnatal human brains (left) and Gene Ontology (GO) biological process (BP) enrichment results according to Metascape for corresponding modules (right). Heatmaps indicate the P values of the degree of enrichment. (C) Integrative analysis of the potential pathogenicity of coding DNM-affected genes. Thirteen independent sources of evidence were used to rank DNM-affected genes. (D) Schematic of the chromatin-regulating activities of SETD5, KDM3B, and ASXL3. (E) Western blot result (left) and cross-species conservation at the discovered mutation site (right). (F) qPCR validation of mutation of (E). CHD8, chromodomain helicase DNA binding protein 8; FMRP, fragile X mental retardation protein; RBFOX, RNA binding fox-1 homolog 1; G2C, the genes to cognition program.

Chromatin modifiers as potential contributors to OCD pathology

As noted, our analysis of mutation enrichment, dSV location, and coding DNM coexpression consistently implicated chromatin regions and chromatin modification–related pathways in OCD. Thus, we hypothesized that a subset of DNM genes related to the highlighted pathways might be more relevant to OCD. Thus, we prioritized all 45 DNM genes (Table 1) according to 13 independent sources of evidence (median Jaccard similarity index = 0.0244) (Fig. 2C; Supplementary Materials and Methods for details). We found three genes—i.e., SETD5, KDM3B, and ASXL3—with the most supporting evidence (≥6) and were all mutation-intolerant (pLI ≥ 0.9) chromatin modifiers that acted as network hubs. SETD5 functions as a histone methyltransferase and monomethylates Lys9 of histone H3, whereas KDM3B is a histone demethylase that explicitly demethylates Lys9 of histone H3. Also, ASXL3 is a part of the polycomb repressive deubiquitinase (PR-DUB) complex that deubiquitinates histone H2A lysine119 (Fig. 2D). Furthermore, as for FBL, a gene affected by de novo structure variant of patient WOC3_114_1_PT, a previous study has already confirmed its ability to methylated histone H2A glutamine-104 and regulated the nucleolar activity (Fig. 2D).

To confirm whether these mutations affected genes’ functions in chromatin modifications as expected, we selected SETD5 and ASXL3 genes to perform wild-type (WT) versus mutant experiments in human embryonic kidney (HEK) 293T cells for the hypothesis verification (Materials and Methods). We first constructed expression plasmids with WT/mutant SETD5 and ASXL3, fused with a FLAG-tag for overexpression experiments in HEK293T cells. Next, we transfected HEK293T cells with plasmids expressing the empty pcDNA3.1 vector (as a control), WT SETD5 or ASXL3, and p.R77C SETD5 or p.F1460Lfs*5 ASXL3. We collected protein lysates from each cell group and examined histone H3K9me1 or H2AK119ub levels via Western blotting after 48 hours. We found that the p.R77C SETD5 mutant’s transfection could not retain the monomethylation level of the histone mark H3K9 compared to the transfection of empty pcDNA3.1 vector or WT SETD5, indicating that p.R77C mutation impaired the enzymatic activity of SETD5 (Fig. 2E). Meanwhile, we observed that transfection of WT ASXL3 decreased the H2AK119 ubiquitination level compared with control with pcDNA, confirming the histone deubiquitination function of WT ASXL3. In contrast, the H2AK119 ubiquitination level of cells with ASXL3 mutant of p.F1460Lfs*5 resembled the pcDNA control, suggesting that this mutant lacked the normal histone deubiquitination function (Fig. 2E).

Chromatin modifiers might involve in OCD through neurotransmitters

OCD has been hypothesized to be caused by alterations in neurotransmitter pathways, such as the glutamate system, dopamine system, and serotonin system. Therefore, we hypothesized that chromatin modifiers identified in the present study—such as SETD5, KDM3B, ASXL3, and FBL—may regulate neurotransmitters’ expression by epigenetic modification. To investigate the expression regulation between chromatin modifiers and neurotransmitters in OCD, we used brain gene expression data from a study conducted by Jaffe et al. (), which contains the expression profile of the dorsolateral prefrontal cortex of postmortem brains isolated from OCD cases and nonpsychiatric controls. We first calculated the Pearson’s Rs between three chromatin modifiers and neurotransmitters in OCD cases and in controls. We then measured the coexpression regulation (ΔCo-exp) by calculating the absolute coexpression difference between patients with OCD and controls (|ΔCo-exp| = |ROCD – Rcontrol|) (Fig. 3A and table S5). We found that the overall coexpression between KDM3B and dopamine genes was significantly altered between OCD cases and controls (Wilcoxon’s test P = 0.025). In contrast, the coexpression between ASXL3 and glutamate and between ASXL3 (Wilcoxon’s test P = 0.047) and serotonin were marginally disrupted (Wilcoxon’s test P = 0.045; Fig. 3A).

Fig. 3.

Disruption of coexpression between key de novo genes.

(A) Coexpression between FBL, SETD5, KDM3B, and ASXL3 and neurotransmitter (glutamate, serotonin, and dopamine)–related genes in the prefrontal cortex of patients with OCD. Each dot represents a gene involved in one neurotransmitter system, and the Y axis shows the absolute differences in the coexpression (ΔCo-exp) between patients and controls. P values were calculated using Fisher’s tests and indicate whether ΔCo-exp was significantly larger than expected. (B) Disruptions caused by KDM3B and SETD5 had opposite trends. For each neurotransmitter system, the ΔR values of SETD5 and KDM3B were plotted against each other, and their correlations were calculated.

By examining the neurotransmitters with the most altered coexpression between controls and OCD cases, we observed that the tyrosine hydroxylase (TH) gene, which is involved in the conversion of tyrosine to dopamine, had the highest disruption by KDM3B (ΔKDM3B = 0.89 from R = 0.22 in control to R = −0.66 in OCD, P = 0.0065; Fig. 3A and table S5). The MAOA gene, which is involved in serotonin and dopamine metabolism, was also found to be disrupted in patients with OCD by KDM3B (ΔKDM3B = 0.66 from R = −0.07 in control to R = 0.59 in OCD, P = 0.06). The coexpression between SETD5 and the serotonin gene HTR1D, also known as 5-HT, was significantly changed (ΔSETD5 = 0.75 from R = 0.35 in control to R = −0.4 in OCD, P = 0.002). The coexpression between FBL and dopamine gene SV2C, which regulated secretion in neural and endocrine cells, was significantly changed (ΔFBL = 0.76 from R = −0.23 in control to R = 0.53 in OCD).

Since SETD5 methylates the histone mark H3K9, and KDM3B demethylates the histone mark H3K9 (Fig. 2D), we hypothesized that SETD5 and KDM3B might regulate the expression of the same neurotransmitters differently. As expected, we observed negative Rs between ΔCo-exp and ΔCo-exp in the glutamate, serotonin, and dopamine systems with R = −0.35, −0.31, and − 0.56, respectively (Fig. 3B). In contrast, ASXL3, which promotes deubiquitination of the histone mark H2A119, a different histone mark, was not correlated with either SETD5 or KDM3B in any of the neurotransmitter regulatory systems.

Integrated analysis of DNMs identified convergent and divergent patterns between OCD and TD

OCD and TD exhibit many similarities, including genetics, phenotypes, and epidemiology (). Thus, to further understand the complex genetics of OCD, we compared and contrasted OCD and TD through the analysis of rare DNMs (Supplementary Materials and Methods). By directly comparing the overlapping genes between OCD and TD, combined from the current study and published exome data (, , –), we observed more gene overlaps between these disorders than expected (all exonic DNMs: P = 1.21 × 10−38; LoF DNMs: P = 9.19 × 10−18; missense DNMs: P = 1.01 × 10−34) (Fig. 4A). Next, we investigated whether the enrichment of OCD and TD genes across 15 disease-related gene sets (Supplementary Materials and Methods) were also similar. We observed a strong positive correlation (Pearson R = 0.54, P = 0.01) between enrichments of the two disorders by comparing the log10-transformed enrichment ORs (Fig. 4B), indicating that a number of functional biological pathways are shared between the two disorders.

Fig. 4.

Comparison of genes affected by coding mutations in OCD and TDA.

(A) Venn plots showing the relationships between TD DNM–affected genes (curated from three studies) and OCD DNM–affected genes [both DNMs from Cappi et al. (, ) and the current study]. (B) Correlations between gene-set enrichments of TD DNM and OCD DNM–affected genes, which were not introduced by similarity or dependency between functional gene sets (average Jaccard similarity index = 0.025). Each dot represents one of the 15 manually curated gene sets (Supplementary Materials and Methods). The X axis represents the OR of TD DNM gene enrichment within this gene set, and the Y axis represents that of the OCD DNM–affected genes. Rs and P values were calculated using Pearson’s correlation analysis. (C and D) Distinct expression patterns of TD and OCD genes among (C) brain regions and (D) brain cell types. Red lines indicate the nominal P value threshold (0.05). TC, temporal cortex; STR, striatum; PFC, prefrontal cortex; PC, parietal cortex; OC, occipital cerebral wall; MD, mediodorsal nucleus of the thalamus; HIP, hippocampus; DTH, dorsal thalamus; CB, cerebellar cortex; AMY, amygdala.

Since a large proportion (75%) of DNM genes in OCD are not shared by TD, we hypothesized that this difference might elucidate specific developmental timing, brain regions, or neuronal cells likely to be involved in OCD pathogenesis. Using BrainSpan data () for developmental periods and brain regions, we observed that the OCD mutation genes were enriched in the dorsal thalamus [EWCE () P = 0.034], whereas the TD mutation genes were enriched in the occipital cortex (P = 0.039; Fig. 4C). On the basis of a study conducted by Dronc () for adult brain cell–type data, the OCD DNM genes were found to be enriched in astrocytes (P = 0.023), whereas the TD mutation genes were found to be enriched in GABAergic neurons (P = 0.006; Fig. 3D). No enrichment differences were observed according to brain developmental periods. These results indicate that there are clear, distinct developmental patterns that separate these two disorders.

DISCUSSION

Motivated by the growing interest in identifying ultra-rare, potentially highly penetrant genetic variants underlying the pathogenesis of neuropsychiatric disorders (, , ), in the present study, we described the first genome-wide DNM profiling of families of trios with OCD by WGS. We identified 4147 genome-wide DNMs, 57 of which affected exonic regions in probands. By exploring the global burden of DNMs in OCD probands, we observed a higher-than-expected proportion of mutations occurring at promoter-anchored chromatin loops and regions with a high frequency of histone marks, such as H3K9me3, H4K20me3, and H3K36me3 (chromatin state 8). When examining the properties of genes directly affected by mutations in OCD, we observed that DNM genes were significantly more mutation-intolerant than SSC controls, consistent with the hypothesis that DNMs might be a significant risk factor for OCD (, ).

Using expression data from multiple regions of the developing human brain, we discovered that the DNM genes identified in our OCD whole-genome screening showed highly correlated gene expression patterns in prenatal and postnatal development separately. The prenatal phase coexpression network included one single enriched module, M16, related to histone modifications and neuronal and synapse organization, whereas the postnatal phase included one single enriched module, M11, highlighting mainly chromatin modifications. These enriched modules from two sequential phases both identified chromatin regulation as a related top pathway, emphasizing its importance in OCD. These results also suggest that additional functional activities are involved in the development of OCD during prenatal development because M16 has other complex neuronal functions, such as neuron projection morphogenesis, synapse organization, and cell-cell adhesion, in addition to epigenetic regulation.

Chromatin-modifier genes are associated with the risk of various neurodevelopmental disorders (, ). In the present study, we observed three genes with top evidence-based rankings in OCD (SETD5, KDM3B, and ASXL3) and a gene disrupted by a damaging dSV (FBL). These genes are mainly chromatin modifiers and are affected by de novo damaging mutations despite their extremely high intolerance to protein-altering variants. SETD5 and KDM3B both regulate H3K9 histone methylation, which is critical in heterochromatin maintenance (). This is consistent with our observation that patients with OCD harbor more mutations in chromatin regions (chromatin state 8) packed with a high frequency of H3K9 histone marks (). LoF mutations in SETD5 are already implicated heavily in ASD () and ID (). Mouse models of SETD5 haploinsufficiency have demonstrated impaired cognition and memory as well as inflexibility in behaviors (), suggesting possible similarities to OCD clinical phenotypes. KDM3B is part of an important group of histone lysine methylases (KMTs) and histone lysine demethylases (KDMs), which are involved in gene regulation and expression. De novo and inherited pathogenic variants in KDM3B can cause a syndrome characterized by ID, short stature, and facial dysmorphism (). The ASXL3 gene (OMIM® number: 615485) is associated with the Bainbridge-Ropers syndrome, a developmental disorder characterized by delayed psychomotor development and severe ID (). The gene encodes a polycomb group protein that is part of the PR-DUB complex, functioning to deubiquitinate H2AUb1 (). Thus, previous genetic studies have established their important connections to various other mental disorders through different epigenetic regulatory mechanisms. FBL acts as a protein methyltransferase by mediating methylation of “Gln105” of histone H2A (H2AQ104me), a modification that impairs the binding of the FACT complex (). It is expressed in both mRNA and proteins in the human brain (), with a pLI mutation intolerance score of >0.95. Although there is no direct link between FBL and psychiatric disorders, previous study has linked SNPs of FBL to phenotypes, such as nicotine metabolite ratio (), suggesting that it might be a potential candidate for studying disorders with addictive behaviors.

An exome study by Cappi et al. () identified a chromatin remodeler gene, CHD8 (encoding chromodomain helicase DNA binding protein 8), as a risk gene for OCD and multiple other psychiatric diseases (, ). Although we did not find DNMs in CHD8 in our study because of the polygenic nature of complex diseases and our sample size, both SETD5 and KDM3B have been found to be binding targets of CHD8 (). This, therefore, confirms the concept of the involvement of chromatin regulation in OCD pathogenesis.

The pathophysiology of OCD is associated with abnormalities in the cortico-striatal-thalamic-cortical (CSTC) circuitry and dysregulation of glutamate, serotonin, and dopamine systems within this network. Recently, an integrative OCD model that incorporates circuitry, neurochemistry, and genetic/epigenetic elements has been proposed by Pauls et al. (), suggesting that individuals with OCD carrying genetic risks may be vulnerable to the impact of environmental factors that may trigger expression modifications of glutamate, serotonin, and dopamine system–related genes through epigenetic mechanisms. Our analysis showed that the coexpression patterns between three system genes and chromatin modifiers were significantly altered in the prefrontal cortex in patients with OCD. These results suggest that chromatin modifications involving SETD5, KDM3B, ASXL3, and FBL may be upstream regulators of neurotransmitter system expression, which controls necessary neurocognitive functions. Disruption of any part of this cascade may lead to abnormal obsessive phenotypes.

Despite only having 53 parent-offspring families in the current study, one of the most consistent observations in our findings is the enrichment of chromatin modification function, indicating highly significant, robust evidence for epigenetic regulation in OCD and suggests potential epigenomic vulnerability in OCD. Consequently, although synaptic proteins are commonly believed to be involved in OCD pathogenesis (), our findings highlight the existence of other contributing mechanisms extending beyond the synapse. Although further large-scale studies are required to prove the pathogenic role of epigenetic regulation by DNMs in OCD, given the moderate sample size and statistical significance in the present study, our observations are credible, considering the fact that similar findings have been reported for ASD and SCZ, along with the fact that CHD8 was identified as a top candidate risk in a recent OCD exome study and is also a chromatin remodeler gene.

Last, our findings provide evidence that OCD shares a strong genetic etiology with TD, although the transcriptional patterns differ by brain region and neuronal cell types in these two neuropsychiatric disorders. Similar enrichment in neurodevelopmental genes and prenatal periods was observed for TD and OCD genes, consistent with similar abnormalities in neurotransmitters’ metabolism () and symptoms. Our findings further support that both diseases are developmental disorders (). However, OCD genes show a uniquely high expression in the dorsal thalamus, consistent with the CSTC circuit in the OCD pathology model. This convergence and divergence between TS and OCD may be interpreted from the aspect of polygenic risk. Given the extensive genetic and phenotypic heterogeneity underlying both diseases and that only DNMs have been examined in the present study, our findings likely represent only a small fraction of the convergent and/or divergent aspects of the two disorders. Nonetheless, our analytic approach provides an important step in understanding the underlying comorbidities between OCD and TD.

MATERIALS AND METHODS

Subjects

We collected 53 unrelated parent-offspring families consisting of 160 individuals based on the availability of genomic DNA from whole-blood and the completeness of phenotype information. The parent-offspring families were composed of 52 trios (each family with one OCD-affected offspring, 156 samples in total) and one family with two OCD-affected offspring (four samples in total). Detailed information for each participant is provided in table S1. All the OCD affected patients met the following conditions: (i) The patients were diagnosed as having OCD according to Diagnostic and Statistical Manual of Mental Disorders, Fourth Edition (DSM-IV) () criteria. (ii) between the age of 18 and 65 years, and (iii)with a Yale-Brown Obsessive-Compulsive Scale total score cutoff of 16. The patients were excluded if they (i) included DSM-IV criteria for other disorders other than OCD, (ii) had moderate to severe suicidal ideation, and (iii) were pregnant or lactating females. Patients diagnosed with SCZ, schizoaffective disorder, ASD, pervasive developmental disorder not otherwise specified, or ID were excluded from the present study. All of the parents were screened for mental disorders by a structured interview using Mini International Neuropsychiatric Interview (), and the families with the parents of any DSM-IV Axis I psychiatric diagnosis were excluded. All participants’ informed consent was obtained, as approved by the Institutional Review Boards of Shanghai Mental Health Center.

WGS, data processing, and DNM identification

Genomic DNA extracted from whole-blood– or lymphoblast-derived cell lines were assessed for quality by PicoGreen and gel electrophoresis and then sequenced by Novogene (Novogene Biosciences Inc., Beijing, China). DNA quantity was measured by Qubit 3.0, and at least 1 μg of nondegraded genomic DNA was used for genomic library preparation and WGS. We sequenced all trio samples, which have never been previously sequenced, on an Illumina HiSeq 4000 sequencing platform (150-bp paired-end reads). Using the processing pipeline, all de novo variants were identified as high quality and then were manually inspected by visualization of aligned reads using the Integrative Genomics Viewer (IGV) (). Variants were annotated by using wANNOVAR (). The average depth of coverage of our WGS data was 29.9×, with an average of 99.87% median alignment rate. After one trio was excluded for quality control, the remaining 53 trios were submitted for subsequent analyses.

Using the Sentieon DNAseq () processing pipeline, reads were aligned to GRCh37.63 human genome reference build using BWA-mem () v0.7.15-r1140 and sorted by Sentieon Dedup. Then, duplicate reads were marked and removed by Sentieon LocusCollector. The indel realignment and base quality score recalibration were performed using Sentieon Realigner and QualCal. Next, we performed the quality check on the BAM files using SAMtools (version 1.3.1) flagstat () and WGA metric by Picard () (version 2.5.0).

De novo SNVs and indels were called using a combination of two algorithms, SAMtools bcftools (version 1.3.1) () and the Bayesian framework TrioDeNovo () with the default settings on a per-family basis. Putative de novo variants were filtered on the basis of the following thresholds: (i) genotype quality ≥ 30; (ii) minimum sequencing depth is 15 read points in the proband and both parents; (iii) mapping quality ≥ 30; (iv) homozygous in father and mother with allele balance (AB) < 0.05; (v) heterozygous in a proband with AB between 0.3 and 0.7; (vi) de novo quality score from TrioDeNovo ≥7; (vii) no overlap with known regions annotated as segmental duplication; (viii) minor allele frequency ≤ 5.0 × 10−3 in 1000 Genomes project (), ExAC (), EVS (), and gnomAD () databases; and (ix) variant clustered within a distance of 50 bp has been removed. Last, all de novo variants were detected with aligned reads from families in silico BLAT search () and then were manually inspected by visualization of aligned reads using IGV ().

SNV and indel annotation

Variants were annotated by four different annotation groups as follows: (1) Variant type: Each variant was classified by mutation types using wANNOVAR (), including SNV and indel (<50 bp), such as frameshift, stop gain (nonsense), nonsynonymous SNV (missense), and synonymous SNV, etc. (2) Genomic location annotation: Genomic locations of variants were annotated by wANNOVAR () and assigned in the following priority: coding, intron, UTR, upstream, downstream, ncRNA, and intergenic. (3) Gene sets: Gene lists associated with brain development or neuropsychiatric disorders were selected as follows: (3.1) CHD8 target genes were defined as the union of lists from two chromatin immunoprecipitation sequencing studies (, ). (3.2) FMRP target genes were defined as the union of lists from Darnell et al. () and Ascano et al. (). (3.3) RBFOX splice targets were selected from Weyn-Vanhentenryck et al. (). (3.4) Human postsynaptic density proteins were extracted from the Genes2Cognition database (). (3.5) Constrained genes were defined as having a pLI score ≥ 0.9 in the ExAC database (). (3.6) Genes encoding chromatin modifiers were downloaded from Chen et al. (). (3.7) Neurodevelopmental disorder risk genes were obtained from Stessman et al. ().

Identification of dSV

Four types of methods have been used to detect SVs: Manta (), cn.MOPS (), DELLY (), and LUMPY (). We applied SV2 () to merge SVs within 50 bp to 10 Mb into a VCF file for each family, and SVs with “PASS” in both FILTER tag and DENOVO_FILTER tag are selected as the candidate dSVs. In addition, we perform different filter conditions for final dSVs as follows: For male patients, we screen for the homozygous SVs that are absent in unaffected parent on autosome and X chromosome or absent in unaffected father on Y chromosome. For a female patient, we screen for the heterozygous/homozygous SVs, which are absent in unaffected parents on autosome and X chromosomes. The potential dSV detected by SV2 was further filtered by two procedures. First, we inspected the IGV of each dSV in the corresponding trios with at least two experts and removed those without corresponding dosage change. Second, we browsed GnomAD database and removed all dSV that had ever appeared (i.e., more than 50% length of dSV was covered by a SV of same type from GnomAD) in the nondisease Han Chinese population. Last, we validated the existence of candidate dSV by qPCR.

Estimation of relative mutation rates and definition of subtypes of coding mutations

To calculate the genome- and exome-wide mutation rate, we first downloaded the ranges of all exonic and UTR regions of hg19 from the UCSC Genome Browser and calculated the total number of base pairs in coding regions with more than 20× coverage across the whole genome. Then, we estimated the average mutation rates in a 95% confidence interval for all 53 patients with OCD using the “t.test” function in R. To determine the differences in DNM distribution between patients with OCD and healthy individuals, we included DNMs from SSC siblings () used as controls. Since the detected number of DNM slightly differed between studies, we applied Fisher exact test to evaluate the ORs, which measure the relative rate of DNMs in healthy controls compared to those in patients with OCD. We then partitioned all coding mutations according to Annovar () annotations: We defined frameshift indel or start site/stop site/splice site mutations as LoF mutations and nonsynonymous SNVs as missense mutations.

Comparison of mutation severity, gene intolerance, and mutation distribution

We compared CADD-phred scores between patients with OCD and SSC siblings from An et al. () for both exonic and genomic mutations. For the evaluation of gene intolerance, we collected pLI () scores for genes affected by different kinds of coding DNMs (LoF, missense, damaging, and synonymous) in patients with OCD. P values were calculated by a two-sided Wilcoxon rank-sum test. Since CADD-phred corresponds to the percentile rank of mutation severity and does not follow a normal distribution, we applied the cumulative distribution curve to visualize the distribution of CADD-phred score and Wilcoxon rank-sum test to compare the predicted severity. To evaluate the influence of GC bias in different genome regions, we compared the GC content and gene length between DNMs in our OCD data and SSC controls. Data for GC content and gene length were collected from GenBank. P values were calculated by a two-sided Wilcoxon rank-sum test.

We compared the results for the OCD patients with those of other neuropsychiatric diseases (developmental delay, ID, SCZ, and ASD) obtained from the PsyMuKB database (). Since detailed information for single individuals in some of these studies was unavailable, we applied a straightforward Fisher’s test on mutation counts adjusted by synonymous mutation.

To investigate the functional genomic distribution of mutations, we first obtained four definitions of functional regions of the human genome from ENSEMBL () and PsychENCODE (): (i) Full-length gene: We obtained ranges of all coding genes from biomaRt with entries “start_position” and “end_position.” (ii) Promoter: For all genes obtained in (i), we found their longest transcripts and the corresponding transcription start sites (TSSs). We defined the promoter region as 2-kb upstream and 1-kb downstream of TSSs. (iii) Promoter Interaction Region: We downloaded from PsychENCODE resource page the promoter anchored loop file, which was obtained from HI-C data of the human prefrontal cortex. (iv) Enhancer: From PsychENCODE resource page, we downloaded the enhancer-gene interaction file generated using HI-C data of the prefrontal cortex and extracted all enhancer regions. We also downloaded the core 15 chromatin state annotations from the Roadmap project (), which was calculated by chromHMM ().

Permutation test

To verify that the mutation distribution enrichment was not driven by random effects, we applied permutation tests to the positive results (promoter loop and chromatin state 8). First, we randomly chose 53 controls and repeated the Fisher’s test. This procedure was repeated 10,000 times, and the distribution of OR was used to generate a P value for OR = 1. Last, we removed all mutations on one chromosome and repeated Fisher’s test on the remaining mutations. This procedure was repeated for each of the 22 autosomal chromosomes.

WGCNA and analysis of coexpression patterns of SETD5, KDM3B, FBL, and ASXL3

BrainSpan RNA sequencing data () were applied to the following coexpression network analysis. We split the expression profiles into two different period sets by prenatal and postnatal samples, including all brain regions. For each subset, only genes with Reads Per Kilobase per Million mapped reads (RPKM) of >0 in at least half of the samples and coefficient of variance of >0.3 were retained for WGCNA () analysis. The remaining data were log10-transformed. The soft threshold was set at 12 for the prenatal subset and 16 for a postnatal subset. The minimum module size was set at 20. We constructed the signed networks by blockwiseModules function in the WGCNA R package based on Pearson R and partitioned genes into modules with a tree cut height of 0.3. After the detection of coexpression modules in both subsets, we tested whether genes carrying OCD coding mutations were enriched in some of the modules by the Fisher’s exact test. The enriched modules were then used for the Gene Ontology enrichment and network analysis.

To explore the potential dysregulation between four key genes (SETD5, KDM3B, FBL, and ASXL3) and three neurotransmitter systems (glutamate, serotonin, and dopamine), we first calculated the Pearson’s R between the key genes and all other genes in the prefrontal cortex expression data from Jaffe et al. (). We used the two-sided Wilcoxon rank-sum test to check whether the overall coexpression between the chromatin modifiers and neurotransmitters was significant between patients with OCD and controls. To see whether any single neurotransmitter gene was significantly dysregulated, we compared the |ΔCo-exp| of each gene to the whole distribution of |ΔCo-exp| of corresponding key genes and generated a nominal P value.

Plasmids

The SETD5-FLAG and ASXL3-FLAG expression plasmids were purchased from Youze Biotechnology Co. Ltd. (Hunan, China), and the SETD5 R77C and ASXL3 F1460Lfs*5 mutated plasmids were introduced using the Q5 Site-Directed Mutagenesis Kit (New England BioLabs). All plasmids were confirmed via Sanger sequencing.

Western blot analysis and antibody

HEK293T cells were lysed in 1× lysis buffer on ice for 30 min and centrifuged at 12,000 rpm for 15 min. Cell lysates were separated by electrophoresis on 4 to 20% SDS–polyacrylamide gel electrophoresis gels and then transferred onto nitrocellulose membranes. The membranes were blocked with 5% milk in TBST (Tris-buffered saline, 0.1% Tween® 20) for 1 hour and incubated overnight at 4°C with the corresponding primary antibodies: rabbit anti-FLAG and anti–β-tubulin (Sigma-Aldrich), rabbit anti-H3K9me1 (Abcam), and rabbit anti-H2AK119ub (Cell Signaling Technology). Next, the experiment was followed by incubation with the horseradish peroxidase–labeled Goat Anti-Rabbit immunoglobulin G (Beyotime) for 1 hour at room temperature 20° to 25°C. Quantitation of immunoblots was then performed via densitometric analysis using the ImageJ software.

Comparison of OCD and TD DNM genes

We downloaded publicly available DNM data for OCD and TD from the psyMuKB database (). A total of two OCD studies (, ) and three TD studies (–) were included in the present dataset. All genes with at least one missense or LoF mutation were recorded. We combined our nonsynonymous mutation genes with published OCD DNM genes and compared them to all previous TD DNM genes.

To analyze the overlap of TD and OCD genes, we performed hypergeometric tests with a background gene list as all protein-coding genes. P values were calculated for two-tailed tests. To analyze the functional characteristics of TD and OCD genes, we first obtained a list of 15 gene sets, mainly about the central nervous system, from psyMuKB (). We tested enrichment TD/OCD genes in all these gene sets by Fisher’s exact test. Then, we performed a Pearson correlation analysis on −log10 (OR) of genes from two diseases. To analyze the spatiotemporal- and cell type–specific expression patterns of TD and OCD genes, we applied EWCE R package (), a bootstrap enrichment tool based on gene cell-type specificity matrix, to conduct enrichment analysis on two datasets: (i) Brainspan () for period and region analysis and (ii) Dronc data () for adult brain cell–type analysis. Specificity was defined in the corresponding paper (). We first calculated the specificity matrix of two expression datasets by “generate.celltype.data” function. Next, the enrichment of all TD/OCD genes was tested by “bootstrap.enrichment.test” function. We defined the background gene list as all genes annotated in the tested expression dataset.