Monozygotic twin patients with treatment-resistant schizophrenia were recruited at Osaka University Hospital. The details of the patients were previously described (4). The twins were 59-year-old Japanese females, both diagnosed with treatment-resistant schizophrenia. Each subject was diagnosed and assessed by at least two trained psychiatrists according to the Diagnostic and Statistical Manual of Mental Disorders, fourth edition (DSM-IV) criteria based on a structured clinical interview. Written informed consent was obtained from subjects after the procedures were fully explained.

The generation of iPS cells from lymphoblastoid B-cell lines and in vitro differentiation of iPS cells into excitatory neurons were performed as previously described (4). In brief, immortalized lymphoblastoid B cells were electroporated with plasmids encoding hOCT3/4, hSOX2, hKLF4, hL-MYC, hLIN28 and dominant negative mTP53. After electroporation, cells were seeded onto SNL feeder cells and grown to form colonies. Subsequently, colonies similar to human ES cells were clonally isolated, morphologically selected and evaluated for the expression of pluripotency markers. For in vitro differentiation, iPS cells were transferred under feeder-free conditions onto Matrigel (Corning)-coated dishes and cultured in the chemically defined medium Essential 8 (Thermo Fisher Scientific) according to the manufacturer's instructions. The in vitro differentiation of iPS cells into neurons through neurogenin2 (Ngn2) overexpression was performed essentially according to Zhang et al. (22). On day 8, 1 μM clozapine (Sigma-Aldrich) or mock treatment was added to the culture medium. After 7 days, Ngn2-induced neurons were harvested for DNA isolation. The gene expression levels in iPS neurons measured by RNA-seq and related statistics were obtained from our previous study (4). In our previous study, total RNA was isolated from three clones obtained from each patient and equal amounts of total RNAs from each patient were combined and sequenced using the Illumina HiSeq2000 system (BGI, Beijing, China) (4). The reads were aligned to the human reference genome hg19. Identifying of DEGs between two samples was performed based on a Poisson distribution (BGI, Beijing, China). The calculated p-values were adjusted using FDR q-value to correct for multiple testing. Genes were identified as differentially expressed between the mock and clozapine treatments or between responder and non-responder at a significance level of FDR q-value < 0.05.

Genomic DNA was extracted from neuronal nuclei using a Blood & Cell Culture DNA kit (Qiagen). We used the same three clones of each twin as those analyzed by our previous RNA-seq analysis (4) for methylation analysis. Bisulfite conversion of 500 ng of genomic DNA was performed with the EZ DNA methylation kit (Zymo Research). DNA methylation levels were assessed according to the manufacturer's instructions using Infinium® HumanMethylation450 BeadChips (Illumina Inc., San Diego, CA, USA), which enable the examination of DNA methylation status at 485,577 CpG sites, and the resulting data were analyzed using the methylation analysis module within the BeadStudio software (Illumina Inc.). The data were mapped to the hg19 genome. The regions described as TSS200 and TSS1500 in this annotation file were defined as the promoter regions. For methylation analysis, IDAT files were processed using the R software package minfi (23). The raw signal data were normalized by the preprocessQuantile function. The methylation status of each CpG site was represented as a β-value that ranged from 0 (completely unmethylated) to 1 (fully methylated). We analyzed 341,629 sites that satisfied the following criteria: (1) β-values with detection p-values < 0.01; (2) CpGs with probes having <3 beads; (3) no probe single nucleotide polymorphisms (SNPs) with minor allele frequencies (MAFs) ≥ 1% in the HapMap-JPT population; (4) no probe cross-reactivity; and (5) no SNPs at CpG sites and single-base extension sites. A list of probes and their corresponding MAF values in the Japanese population was derived from Okamura et al. (24). We excluded sites with SNPs and in cross-reactive regions based on a list derived from Chen et al. (25). We performed dmpFinder function in the minfi package to identify positions that were differentially methylated. Methylation sites were identified as differentially methylated between the mock and clozapine treatments or between responder and non-responder at a significance level of p-value < 0.01. The calculated p-values were adjusted using FDR q-value to correct for multiple testing, although there were no methylation sites differentially methylated between the mock and clozapine treatments at a significance level of FDR q-value < 0.05.

The same direction means that both the mRNA expression level and the DNA methylation level are higher or lower in clozapine compared to mock treatment. The opposite direction indicates highly express and lower methylate in clozapine compared to mock treatment, or vice versa.

Genes associated with neuronal and synaptic functions were defined as those genes annotated with Gene Ontology terms including the following words: “neuro,” “synapse,” or “synaptic.” We used the information from the Gene Ontology project (http://geneontology.org/) and R software package GO.db. We obtained 2,756 genes after the procedures. Furthermore, we selected 2,399 genes overlapped with 15,466 genes that we treated in the RNA-seq and DNA methylation analyses.

Gene functional enrichment analysis was performed using Metascape software (26). Metascape reports a term that an input gene list overrepresented as an enrichment cluster and prevents redundancies in terms across different ontology sources.

To examine the extent of overlap between the DEGs and the genes related to neuronal and synaptic functions, we calculated the p-value by hypergeometric distribution test and the expected numbers as follows:

where x is the number of genes that overlapped between the DEGs and the genes related to neuronal and synaptic functions, and m and n are the numbers of the DEGs and the genes related to neuronal and synaptic functions (2,399 genes), respectively. N is the total number of genes (15,466 genes that we treated in the RNA-seq and DNA methylation analyses).

We show an overview of the present study in Figure 1. We performed genome-wide DNA methylation profiling using differentiated neurons from iPS cell clones obtained from clozapine responders and non-responders. We analyzed 341,629 of 485,577 methylation sites based on robustness criteria (see section Methods). To identify positions that were differentially methylated following clozapine treatment, we compared the DNA methylation levels of each methylation site between the mock and clozapine-treated neurons in each twin. We extracted the changed methylation sites from each comparison at a significance level of p-value < 0.01. As a result, we identified 2,229 and 1,384 methylation sites in the clozapine responder and non-responder, respectively (Supplementary Tables 1, 2). Among these changed methylation sites, only 21 sites were overlapped between the twins. After we excluded the 21 overlapped sites from the changed methylation sites in each twin, these patient-specific methylation sites were involved 1,789 and 1,145 genes in the clozapine responder and non-responder, respectively. We used these methylation data in the subsequent analyses.

We recently identified DEGs following clozapine treatment in iPS neurons from the same twins by performing RNA-seq analysis (4). We found 714 DEGs specific to the clozapine responder and 277 DEGs specific to the clozapine non-responder (FDR q-value < 0.05). Among 714 DEGs identified in the clozapine responder, 84 genes were located near the responder-specific methylation sites (Supplementary Table 3). Of the 84 DEGs, 17 DEGs, including MECP2 (methyl-CpG-binding protein 2), which has been implicated in schizophrenia and autism as well as Rett syndrome (27, 28), differed in the opposite direction of the change in the methylation sites in the promoter regions (Table 1; Supplementary Figure 1), and 32 DEGs showed the same direction of change for the mRNA expression level and the DNA methylation level in the corresponding gene body (Table 2). Of the DEGs associated with the promoters and gene bodies, 2 genes [STOX2 (storkhead box 2) and TBC1D16 (TBC1 domain family member 16)] were common. Among 277 DEGs identified in the clozapine non-responder, 15 genes were located near the non-responder-specific methylation sites (Supplementary Table 4). Of these 15 DEGs, 2 DEGs differed in the opposite direction of the change in the methylation sites in the promoter regions (Table 1), while 3 DEGs showed the same direction of change in the mRNA expression level and the DNA methylation level in the corresponding gene body (Table 2).

Then, we examined the proportion of genes associated with neuronal and synaptic functions (Table 3). Among the genes near the changed methylation sites, 300 of the 1,789 genes (16.8%) in the clozapine responder and 190 of the 1,145 genes (16.6%) in the clozapine non-responder were associated with neuronal and synaptic functions. Among the genes with altered methylation and expression, 30 of the 84 genes (35.7%) in the clozapine responder and 1 of the 15 genes (6.7%) in the clozapine non-responder were associated with neuronal and synaptic functions. The number of the genes with neuronal and synaptic functions in the clozapine responder was significantly 2.3 times higher than the expected value (p-value = 2.96E-06; #observed genes / #expected genes = 2.30), whereas the number of the genes in non-responder was lower than the expected value (p-value = 0.220; #observed genes / #expected genes = 0.430). When we performed gene functional enrichment analysis to investigate the functions of these 84 DEGs, we found axon guidance (FDR q-value = 1.59E-07) and neuron projection morphogenesis (FDR q-value = 1.21E-05) in the gene functional enrichment clusters (Supplementary Figure 2). In contrast, we did not find any statistically significantly enriched function in the clozapine non-responder. Among the genes whose expression changes were consistent with the methylation changes in the gene body or promoter regions, 19 of the 47 genes (40.4%) in the clozapine responder were significantly associated with neuronal and synaptic functions (p-value = 2.53E-05; #observed genes / #expected genes = 2.61), while there were no associations in the clozapine non-responder (p-value = 0.430; #observed genes / #expected genes = 0.00).

We finally examined the differences of the basal methylation states of the responder and non-responder. We identified 3,065 and 2,980 methylation sites that changed between the twins in the mock and clozapine treatments at a significance level of p-value < 0.01, respectively. Among these changed methylation sites, 308 sites were overlapped between the two treatments. Each changed methylation sites were involved 2,305 and 2,241 genes in the mock and clozapine treatments, respectively. The RNA-seq analysis found 1,705 DEGs in the mock treatment and 806 DEGs in the clozapine treatment (FDR q-value < 0.05). Among 1,705 DEGs identified in mock treatment, 133 genes were located near the 3,065 changed methylation sites. Among 806 DEGs identified in the clozapine treatment, 88 genes were located near the 2,980 changed methylation sites. The 133 DEGs in mock treatment were overlapped with only 5 genes among the 84 DEGs response to clozapine in the clozapine responder. Additionally, there were no genes overlapped between 133 DEGs in mock treatment and 15 DEGs in the clozapine non-responder.