Peripheral blood samples were obtained from Ningbo Kangning Hospital (Zhejiang, China) in 2011. Patients were clinically diagnosed by trained psychiatrists using the Diagnostic and Statistical Manual of Mental Disorders, 4th Ed. (30). Case 1 refers to a male patient with paranoid SCZ and case 2 refers to a male patient with undifferentiated SCZ. Patients were verified to be without any serious or unstable medical illness. A healthy male individual was used as a negative control. Table I lists the demographic and clinical characteristics of the subjects in the present study. The negative control male and the supervisors of SCZ patients provided their written informed consent for the study. The present study was approved by the Human Research Ethics Committees of Ningbo University and Ningbo Kangning Hospital (Zhejiang, China).

Methods for MBD-seq were similar to those described in a previous study (31). Genomic DNA was isolated from peripheral blood samples of the three individuals using Qiagen Puregene kit (Qiagen, Hiden, Germany). DNA was then sonicated into 200- to 400-bp fragments by 20 cycles of 30 sec sonication, using a Bioruptor sonicator (Diagenode, Liège, Belgium). Subsequently, the active motif methyl collector kit (Qiagen) was used to complete MBD2 enrichment. According to the manufacturer's protocol, ~1 mg of sonicated genomic DNA was incubated with MBD2-His-conjugated protein and magnetic beads. DNA was purified for each biological replicate, and then pooled post-enrichment. After enrichment, DNA purification columns (Qiagen) were used to purify both the methylated fraction and supernatant fractions, as reported previously (31). Finally, MBD2-enriched genomic DNA from each sample was submitted for high-throughput sequencing analysis. Libraries were sequenced twice using an Illumina Hiseq 2000 platform (Illumina, San Diego, CA, USA), obtaining two datasets of 51- and 101-bp single-end reads, which were used separately in order to identify the MBD2-bound fraction of the human genome, respectively.

MBD-seq reads were filtered using FASTX-Toolkit software with a cutoff of <Q20 in order to obtain high-quality reads. Remaining reads were then mapped to the human genome (hg19; ensembl) using bowtie and <2 mismatches (32). DNA methylation peaks for each sample were identified using MACS software (Miltenyi Biotech, Bergisch Gladbach, Germany) (33) with a threshold of P<10–5. Reads from the healthy sample were used as a negative control for the experiments.

Genomic regions of each peak were compared with the known genes in the hg19 reference genome, including coding and lncRNAs, using the RefSeq database (http://www.ncbi.nlm.nih.gov/refseq/). The regions were assigned to one of six classes associated with known genes using HOMER software (http://homer.salk.edu/homer/microarray/) (34) The six classes of genome regions were intergenic, 3′ untranslated regions (UTR), intron, exon, 5′UTR, promoter-transcript start site (promoter-TSS) and transcript termination site (TTS). Promoter-TSS were identified from −1000 bp to +100 bp, while TTS was identified from −100 bp to +1000 bp. Subsequently, intergenic peaks identified using the RefSeq database were compared with genomic regions of lncRNAs from NONCODE version 4.0 database (http://www.noncode.org/). Images of peak visualization were obtained using Drompa 2.3.1 version software (35).

Genomic locations of CpG islands were downloaded from UCSC (https://genome.ucsc.edu/) (36). Genomic regions of each peak were compared with those of CpG islands in order to identify peaks overlapping with or containing CpG islands. The number of CG sites was counted and compared with that in random genomic regions with the same length distribution.

GO terms and their corresponding depth levels in the GO tree were obtained from the GO database (http://geneontology.org/page/go-database). For each GO term, the P-value was calculated using a hyper-geometric test and adjusted with false rate discovery using the 'multtest' package in R statistics. Enriched GO terms with adjusted P-values <0.05 were then obtained. Enriched Kyoto Encyclopedia of Genes and Genomes (http://www.genome.jp/kegg/) pathways were obtained.

Sequencing was repeated twice and two datasets of MBD-seq were obtained for each SCZ sample. Sample reads were filtered using a cutoff of Q20, which resulted in 40–75% of the original reads, of which 65–90% matched with hg19 human genome reads (Table II). Using MACS software for peak calling, 1,236 and 1,552 peaks were obtained in the first dataset from patients with paranoid and undifferentiated SCZ, respectively. In the second dataset, 1,397 and 1,437 peaks were obtained from patients with paranoid and undifferentiated SCZ, respectively. The results of two datasets were compared by HOMER software with default parameters. Compared with the second dataset, there were 12.14 and 12.37% different peaks in the first dataset of paranoid and undifferentiated SCZ, respectively. Peak regions from the two datasets were then merged and peaks shared by the two datasets were selected using HOMER software. Finally, 757 and 871 peaks were obtained from patients with paranoid and undifferentiated SCZ, respectively. The average peak length was 288 and 431 bp in patients with paranoid SCZ and undifferentiated SCZ, respectively. As demonstrated in Fig. 1, chromosome distributions in the two SCZ subtypes were similar (pearson correlation test: r=0.84, P<0.01). G+C contents and CG percentages were similar in the patient with paranoid SCZ (G+C=49.0%; CG=1.9%) and that with undifferentiated SCZ (G+C=51.5%; CG=2.7%). By contrast, G+C and CG contents of random sequences with the same length distribution were markedly lower in the patients with paranoid SCZ (G+C=37.3%; CG=0.8%) and undifferentiated SCZ (G+C=38.4%; CG=0.9%), compared with those in the control patient. Three and 19 peaks overlapped with CpG islands in the paranoid and the undifferentiated SCZ sample. Among them, four peaks in the undifferentiated SCZ sample were located in the gene promotor regions, including tripartite motif containing 45 (TRIM45), major histocompatibility complex, class II, DP α 1 (HLA-DPA1), thioesterase super-family member 6 (THEM6) and chromosome 6 open reading frame 123 (C6orf123). As demonstrated in Fig 2, the intergenic region was the most enriched (paranoid SCZ, 54.96%; undifferentiated SCZ, 44.20%), followed by the intron region (paranoid SCZ, 40.69%; undifferentiated SCZ, 43.97%).

There were 291 coding genes, 18 non-coding genes and 1 pseudogene in the paranoid SCZ sample, whilst 411 coding genes, 21 non-coding genes and 10 pseudogenes were observed in the undifferentiated SCZ sample. GO term enrichment analysis demonstrated that the functions of protein-coding genes in the paranoid SCZ sample were predominantly associated with the nervous system. The five most enriched GO terms were axon guidance, neuron development, neuron projection development, axonogenesis and neuron projection morphogenesis (P<0.01; Fig. 3A). A number of these GO terms are associated with SCZ. For example, contactin-associated protein-like 2 has been shown to mediate cell-cell interactions in the nervous system and is associated with SCZ (37,38) and synapsin III, a member of the synaptic vesicle-associated protein family, was shown to be a candidate gene for SCZ (39). Furthermore, hyper-methylated genes have been shown to be enriched in several important pathways, including the Wnt signaling pathway (40), calcium signaling pathway (41) and a number of pathways associated with cancer (P<0.05).

The functions of hypermethylated protein-coding genes observed in the undifferentiated SCZ sample were associated with the nervous system, including nervous system development, axonogenesis, neuron differentiation and cell morphogenesis (Fig. 3B). There were 32 hypermethylated genes shared by paranoid and undifferentiated SCZ data-sets (Table III). GO term enrichment analysis demonstrated that these shared hypermethylated genes were enriched in the following functions: Axon guidance, axonogenesis, neuron projection morphogenesis and neuron projection development. The results of the present study suggested that SCZ may be caused by aberrant hypermethylation of a number of genes associated with the nervous system.

There were 7 (0.92%) and 12 (1.38%) peaks observed in the promoter regions in the paranoid and undifferentiated SCZ samples, respectively. The 7 genes in the paranoid SCZ sample were zinc finger protein 641, RNA binding protein, fox-1 homolog (RBFOX1), solute carrier family 45 member 2, ubiquitin-conjugating enzyme E2Q family-like 1, cluster of differentiation 8a, retinoic acid early transcript 1K pseudogene and GRAM domain containing 3. RBFOX1 (Fig. 4) is a type of RNA binding protein that is capable of regulating alternative splicing. RBFOX1 binds with ataxin-2, which causes familial neurodegenerative diseases. Therefore, RBFOX1 is involved in neuron development (42). It was hypothesized that the hypermethylated RBFOX1 promoter may affect RBFOX1 expression, thereby regulating biological processes involved in SCZ.

For the undifferentiated SCZ sample, genes with hypermethylated promoters included claudin 5, prostaglandin reductase 1, HLA-DPA1, C6ORF123, THEM6, phosphodiesterase 8B (PDE8B), zinc finger protein 622, TRIM45, C11ORF16, polymerase (RNA) I polypeptide E, tensin like C1 domain containing phosphatase and potassium voltage-gated channel, delayed-rectifier, subfamily S, member 2. Among these, PDE8B (Fig. 4) has been reported to be associated with SCZ (43). A study has demonstrated that PDE8B expression is upregulated following treatment with olanzapine and down-regulated following treatment with lithium, which are two drugs that are used to treat schizoaffective disorders (44).

A total of 17 hyper-methylated lncRNAs were identified in the paranoid SCZ sample. Among them, 3 were shown to be associated with cancer (disrupted in renal carcinoma 3, cancer susceptibility candidate 2 and plasmacytoma variant translocation 1; data available on request). SCZ patients are more likely to suffer from cancer compared with the general population. Therefore, SCZ may lead to alterations in epigenetic modifications for a number of cancer-associated genes and, therefore, increase the risk of cancer (45). In the undifferentiated SCZ sample, 20 lncRNAs were identified, including 5 anti-sense RNAs to protein-coding genes, which were shown to be associated with SCZ. Sodium leak channel non-selective-anti-sense RNA 1 (NALCN-AS1; data available on request) is the anti-sense of NALCN that has been reported as a SCZ candidate gene (46–48). Microtubule-associated protein tau-AS1 (MAPT-AS1; data available on request) is the anti-sense lncRNA gene of MAPT that has been shown to be associated with neurodegenerative disorders and SCZ (49). Guanine nucleotide binding protein (G protein) α stimulating activity polypeptide 1-AS (GNAS-AS; data available on request) is the anti-sense lncRNA gene of GNAS1 that may influence a patients' susceptibility to SCZ (50). POU domain, class 4, transcription factor 1-AS1 (POU4F1-AS1; data available on request) is the anti-sense of POU4F1 that exhibits functions associated with the nervous system (51).

The development of DNA sequencing technology has enabled the identification of an increasing number of human lncRNAs. Although their functions are currently unknown, lncRNAs are likely to be transcribed and be involved in a number of biological processes. The present study has identified a number of peaks in the intergenic regions that may correspond to novel lncRNAs. Therefore, lncRNAs were obtained from NONCODE v4.0 and their genomic regions were compared with the peaks of the intergenic regions identified in the present study. The results identified 131 out of 416 (31.5%) and 114 out of 385 (29.6%) intergenic peaks in the paranoid and undifferentiated SCZ, respectively, corresponding to 434 and 320 lncRNAs from NONCODE database for paranoid and undifferentiated SCZ, respectively. The results of the present study suggested that intergenic hypermethylation may influence lncRNA expression and function. The results of the present study showed that lncRNA n373859 was hypermethylated in the promoter region in paranoid SCZ. As shown in the NONCODE database, lncRNA n373859 is highly expressed in brain tissue. Further study is warranted to investigate whether the aberrant promoter hypermethylation of lncRNA n373859 may lead to downregulation in patients with paranoid SCZ.