DNA methylation patterns of protein-coding genes and long non-coding RNAs in males with schizophrenia.

Schizophrenia (SCZ) is one of the most complex mental illnesses affecting ~1% of the population worldwide. SCZ pathogenesis is considered to be a result of genetic as well as epigenetic alterations. Previous studies have aimed to identify the causative genes of SCZ. However, DNA methylation of long non-coding RNAs (lncRNAs) involved in SCZ has not been fully elucidated. In the present study, a comprehensive genome-wide analysis of DNA methylation was conducted using samples from two male patients with paranoid and undifferentiated SCZ, respectively. Methyl-CpG binding domain protein-enriched genome sequencing was used. In the two patients with paranoid and undifferentiated SCZ, 1,397 and 1,437 peaks were identified, respectively. Bioinformatic analysis demonstrated that peaks were enriched in protein-coding genes, which exhibited nervous system and brain functions. A number of these peaks in gene promoter regions may affect gene expression and, therefore, influence SCZ-associated pathways. Furthermore, 7 and 20 lncRNAs, respectively, in the Refseq database were hypermethylated. According to the lncRNA dataset in the NONCODE database, ~30% of intergenic peaks overlapped with novel lncRNA loci. The results of the present study demonstrated that aberrant hypermethylation of lncRNA genes may be an important epigenetic factor associated with SCZ. However, further studies using larger sample sizes are required.

Introduction

Schizophrenia (SCZ) is considered a heritable disorder that affects brain responses to external stimuli (1,2). This abnormality may cause certain psychological disorders and cause a patient to withdraw from society. The precise etiological factors associated with SCZ remain unknown. In clinical practice, SCZ diagnosis predominantly depends on symptoms that are numerous and debilitating (3). A number of SCZ sub-types are distinguishable according to the symptoms that a patient exhibits. Paranoid and undifferentiated SCZ are the most common SCZ sub-types (4). However, the causes of these two SCZ sub-types have remain elusive.

Studies have demonstrated that a number of susceptibility genes may be associated with specific illness sub-types (4). For example, the dystrobrevin-binding protein 1 gene has been shown to be associated with negative and cognitive symptoms in patients with SCZ (4). SCZ is a complex disorder that may be caused by genetic and epigenetic factors (5). The symptoms associated with SCZ are variable. Therefore, there is a requirement to identify biomarkers that may indicate a patients' genetic pre-disposition to SCZ and the possible environmental exposure risks associated with SCZ.

DNA methylation, one of the most common and important epigenetic factors, is influenced by environmental factors and regulates gene expression (6). Numerous studies have demonstrated that DNA methylation is involved in SCZ (7–18). For example, Reelin is a protein which is responsible for neuronal connectivity and synaptic plasticity (19). Abnormal expression of Reelin may result in cognitive deficits in SCZ (19). In brain samples from patients with SCZ, mRNA and protein levels of Reelin were shown to be reduced by ~50% in response to Reelin promoter hypermethylatation compared to those in control samples (8). However, to the best of our knowledge, studies on DNA methylation in SCZ have predominantly focussed on protein-coding genes. Few studies have investigated the involvement of long non-coding RNAs (lncRNAs) in SCZ (20,21).

lncRNA is a type of non-coding RNA that is >200 nucleotides long. Sequence characteristics are similar to those of mRNAs; however, they lack encoding capacity. In mammals, a number lncRNAs have been shown to be associated with SCZ (22,23). The functions of the majority of lncRNAs are unknown. However, studies have demonstrated that lncRNAs are regulators of expression, sub-cellular location and activity of certain protein-coding genes. lncRNAs are associated with a number of biological processes via protein-coding gene regulation, including genomic imprinting control, cell differentiation, immune responses, human diseases and tumorigenesis (24–26). Therefore, alterations in lncRNA gene methylation may result in differential expression patterns and, therefore, affect the expression of target protein-coding genes. Usually, these lncRNAs are anti-sense transcripts and their promoters are located in the intronic region of the corresponding protein coding genes. DNA methylation in lncRNA gene promotor regions may affect the expression of their corresponding target genes by influencing histone modification, which is mediated by binding to histone modification proteins, including polycomb repressive complex 2 and protein regulator of cytokinesis 1 (27). DNA methylation patterns of lncRNA genes have been investigated in developmental biology and in diseases, such as cancer (28,29). However, to the best of our knowledge, no studies on DNA methylation patterns of lncRNA genes in SCZ have been preformed. In the present study, it was hypothesized that epigenetic alterations of lncRNAs may be associated with SCZ.

In the present study, methylated DNA-binding domain-sequencing (MBD-Seq) was used in order to identify DNA hyper-methylated regions in male patients with paranoid or undifferentiated SCZ. DNA methylation levels of protein-coding genes and lncRNAs were measured in peripheral blood samples in order to explore epigenetic biomarkers for diagnosis and provide a foundation for the development of novel anti-psychotic drugs for SCZ.

Materials and methods

Sample collection

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).

Table I

Demographic and clinical data.

Characteristic	Control	Case 1	Case 2
Gender	Male	Male	Male
Age (years)	25	34	22
Diagnostic subtype	NA	Paranoid	Undifferentiated
Age at onset	NA	19	11
Family history	NA	Negative	Positive
Antipsychotic drug	NA	Quetiapine	Clozapine
Psychotic trauma	NA	Negative	Negative

MBD enrichment and next-generation sequencing

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.

Data processing of MBD-seq

MBD-seq reads were filtered using FASTX-Toolkit software with a cutoff of 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.

Annotation of DNA methylation peaks

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).

CpG island analysis

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.

Gene ontology (GO) and pathway enrichment analysis

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.

Results

Statistical information from the MBD-seq dataset

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%).

Table II

Methylated DNA-binding domain-sequencing datasets.

Data	Control	Paranoid SCZ	Undifferentiated SCZ
Total reads (n)	30,901,427	8,720,205	20,384,242
Reads post-Q20 filter, n (% of total)	13,800,627 (45)	4,017,085 (47)	9,003,647 (45)
Read/match ratio (%)	75.69	81.43	76.34
Unique read/match ratio (%)	66.48	72.83	67.61
Peaks (n)		1,236	1,552

SCZ, schizophrenia.

Figure 1

Chromosomal distribution of DNA hyper-methylation peaks. SCZ, schizophrenia; M, mitochondria.

Figure 2

Distribution of gene function classes for the DNA hypermethylation peaks. UTR, untranslated region; TTS, transcript termination site; SCZ, schizophrenia.

Hypermethylation peaks in SCZ are predominantly enriched in genes characterized by neuron function

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).

Figure 3

Most enriched gene ontology terms of protein-coding genes exhibiting DNA hypermethylation. SCZ, schizophrenia.

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.

Table III

Hypermethylation genes shared between paranoid and undifferentiated SCZ samples.

Gene ID	Gene name	Gene definition
2139	EYA2	Eyes absent homolog 2 (Drosophila)
2272	FHIT	Fragile histidine triad
11240	PADI2	Peptidyl arginine deiminase, type II
100507421	TMEM178B	Transmembrane protein 178B
5332	PLCB4	Phospholipase C, β 4
5184	PEPD	Peptidase D
84691	FAM71F1	Family with sequence similarity 71, member F1
400941	LINC00487	Long intergenic non-protein coding RNA 487
23261	CAMTA1	Calmodulin binding transcription activator 1
84131	CEP78	Centrosomal protein 78kDa
29994	BAZ2B	Bromodomain adjacent to zinc finger domain, 2B
121256	TMEM132D	Transmembrane protein 132D
57221	KIAA1244	KIAA1244
3680	ITGA9	Integrin, α 9
643650	LOC643650	Uncharacterized LOC643650
3899	AFF3	AF4/FMR2 family, member 3
10395	DLC1	Deleted in liver cancer 1
3786	KCNQ3	Potassium voltage-gated channel, KQT-like subfamily, member 3
29119	CTNNA3	Catenin (cadherin-associated protein), α 3
1002	CDH4	Cadherin 4, type 1, R-cadherin (retinal)
54715	RBFOX1	RNA binding protein, fox-1 homolog (Caenorhabditis elegans) 1
27086	FOXP1	Forkhead box P1
169044	COL22A1	Collagen, type XXII, α 1
256380	SCML4	Sex comb on midleg-like 4 (Drosophila)
25771	TBC1D22A	TBC1 domain family, member 22A
10426	TUBGCP3	Tubulin, gamma complex associated protein 3
124540	MSI2	Musashi homolog 2 (Drosophila)
51151	SLC45A2	Solute carrier family 45, member 2
10404	CPQ	Carboxypeptidase Q
125336	LOXHD1	Lipoxygenase homology domains 1
80216	ALPK1	α-kinase 1
92293	TMEM132C	Transmembrane protein 132C

All genes are protein-coding with the exception of LINC00487 and LOC643650. SCZ, schizophrenia.

Protein-coding genes with hypermethylated promoters

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.

Figure 4

Hypermethylated promoter of (top) the PED8B gene in the undifferentiated SCZ samples and the negative control and (bottom) RBFOX1 gene in the paranoid SCZ samples and the negative control. Methylation levels surrounding the peaks from −5,000 to +50,000 bp; the methylation peaks are located in the center. PDE8B, phosphodiesterase 8B; RBFOX, RNA-binding protein, fox-1 homolog; SCZ, schizophrenia.

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).

lncRNA genes and DNA hypermethylation

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.

Discussion

DNA methylation is associated with epigenetic modifications and may be an important factor for SCZ. Aberrant hypermethylation may contribute to abnormal gene expression and influence the development of SCZ. Although a number of studies have reported the effects of numerous important coding genes on DNA methylation in SCZ, including reelin (8,9), catechol-O-methyltransferase (10,11), sex determining region Y (SRY)-box 10 (12), brain-derived neurotrophic factor (16), dopamine receptor D4 (17) and cytotoxic T-lymphocyte-associated protein 4 (17), there are limited studies on lncRNAs. In the present study, a number of hypermethylated loci were identified in protein-coding genes and lncRNAs in male patients with SCZ. The results have provided novel insights for the understanding of SCZ etiology.

In the present study, male patients were selected in order to investigate candidate hypermethylated lncRNAs associated with SCZ. Gender differences are common for a number of aspects of SCZ, such as incidence and prevalence, the presentation and course of the illness, and the various effects on brain development (52). The onset of SCZ is delayed in females and there is typically a second peak in female patients at >45 years of age (53). Estrogen level reduction is considered to be associated with this process (54). Males typically suffer from SCZ more severely and are less responsive to medication than females (52). A number of reasons to explain the early onset in men compared with women have been proposed, including slower maturation and differences in lateralization of the brain, greater incidence of head trauma and the lack of neuroprotective effects from female hormones (52). Therefore, males are more susceptible to SCZ than females (52,55–57). Studies of hormones in female mammals have revealed that psychotic symptoms may vary throughout the menstrual cycle, and that estrogen may enhance neurocognitive performance and improves symptoms of SCZ (58,59). Due to the gender differences associated with SCZ, gender-specific characteristics and endocrine gland side effects require further research in order to develop anti-psychotic treatments. Advances in epigenetic research in SCZ will provide novel insights into SCZ pathogenesis and the development of medication for patients with SCZ.

SCZ pathophysiology research is often based on individual genes and the environment (60). An epidemiological study suggested that interactions between genes and the environment are involved in the molecular mechanisms underlying the onset of SCZ (60). Epigenetic factors, including DNA methylation, genomic imprinting, histone modifications and non-coding RNA expression are involved in SCZ pathogenesis (60–63). These factors may contribute to the mechanisms associated with gender-specific genes and gender dimorphism in SCZ (61).

DNA extracted from peripheral blood samples has been shown to be informative for SCZ research. Dysregulation levels of a number of SCZ candidate genes in peripheral blood samples are associated with a number of clinical symptoms of SCZ (63,64). In the present study, which used genome-wide methylation analyses, DNA methylation levels of protein-coding genes and lncRNAs in peripheral blood samples were investigated. Epigenetic biomarkers were identified that may be used for the diagnosis and development of anti-psychotic drugs for patients with SCZ.

The number of peaks that exhibited DNA hypermethylation were similar between the different SCZ subtypes (paranoid and undifferentiated). Hypermethylated genes were shown to exhibit functions associated with the nervous system in paranoid and undifferentiated SCZ. This suggested that DNA methylation may affect the activity or expression of neuron-associated genes, and thereby influence the development of SCZ. Hypermethylation of promoters may affect the transcription of corresponding genes. In the present study, 7 and 12 peaks were identified in the promoter regions of Refseq genes in paranoid and undifferentiated SCZ samples, respectively. The expression of these genes may be affected by changes to DNA methylation and, therefore, influence biological processes associated with SCZ. With the exception of the peaks corresponding to coding genes, a number of peaks were shown to correspond with known lncRNAs. Furthermore, ~30% of the intergenic peaks that were identified corresponded with novel lncRNAs, according to the RefSeq database. These results suggested that lncRNAs may represent an important regulator for SCZ pathogenesis.

In conclusion, hypermethylated coding and non-coding genes were identified in samples from males with paranoid and undifferentiated SCZ. The results of the present study may clarify the mechanisms underlying SCZ pathogenesis.