Methylome-wide Association Study of Patients with Recent-onset Psychosis.

Objective: Dysregulation of gene expression through epigenetic mechanisms may have a vital role in the pathogenesis of schizophrenia (SZ). In this study, we investigated the association of altered methylation patterns with SZ symptoms and early trauma in patients and healthy controls.
Methods: The present study was conducted to identify methylation changes in CpG sites in peripheral blood associated with recent-onset (RO) psychosis using methylome-wide analysis. Lifestyle factors, such as smoking, alcohol, exercise, and diet, were controlled.
Results: We identified 2,912 differentially methylated CpG sites in patients with RO psychosis compared to controls. Most of the genes associated with the top 20 differentially methylated sites had not been reported in previous methylation studies and were involved in apoptosis, autophagy, axonal growth, neuroinflammation, protein folding, etc. The top 15 significantly enriched Kyoto Encyclopedia of Genes and Genomes pathways included the oxytocin signaling pathway, long-term depression pathway, axon guidance, endometrial cancer, long-term potentiation, mitogen-activated protein kinase signaling pathway, and glutamatergic pathway, among others. In the patient group, significant associations of novel methylated genes with early trauma and psychopathology were observed.
Conclusion: Our results suggest an association of differential DNA methylation with the pathophysiology of psychosis and early trauma. Blood DNA methylation signatures show promise as biomarkers of future psychosis.

INTRODUCTION

Schizophrenia (SZ), which affects 20 million people worldwide [1], is one of the most common types of serious mental illness. To date, attempts to discover genes that are directly associated with SZ have been largely unsuccessful [2,3]. A recent genome-wide association study (GWAS) of SZ revealed 108 loci associated with the disease [4]. However, the proportion of the variance in liability explained by single nucleotide polymorphisms (SNPs) is small; significant genome-wide loci explain only 3.4% of the variance in liability, and the cumulative effect of common loci expressed as a polygenic risk score (PRS) was estimated to explain only 7%.

Therefore, both genetic and environmental factors must be considered in research on the etiology of SZ. Epidemiological studies have suggested that multiple environmental factors, including prenatal infection/immune activation, paternal age, malnutrition, hypoxia-related obstetric complications, childhood/adolescent trauma and cannabis abuse are associated with an increased risk of SZ [5]. The interaction of risk genes with the environment can result in aberrant epigenetic mechanisms in SZ. After first being proposed by Petronis (2004) [6], the epigenetic component of the pathophysiology of SZ has been widely recognized in recent years. DNA methylation, one of the main epigenetic mechanisms, can be examined by utilizing two complementary approaches. The first involves hypothesis-driven analysis of DNA methylation in candidate genes. Using this approach, abnormal methylation of some genes, such as DRD2 [7], DLGAP2 [8], and COMT [9], has been implicated in the onset and progression of SZ. The second approach involves microarray-based epigenomic profiling. With this microarray-based technique, assessment of DNA methylation in about 450,000 or 853,307 CpG sites can be performed. However, this represents only a small fraction of all CpG sites (∼28 million) [10]. Several methylome-wide association studies (MWAS) have been performed using postmortem brain tissues [11-14] and the peripheral blood cells of patients with SZ [15-20]. These MWAS have revealed differences genome-wide methylation patterns between SZ patients and controls, but the number of CpGs analyzed, and the sample sizes, have been limited. Two large-scale methylome studies linked SZ to hypoxia and infection [20], and some differentially methylated positions were located in a top-ranked SZ region based on GWAS analyses [18]. Few MWAS have been performed in first-episode SZ [19,21-23]. Moreover, with a few exceptions [18,24,25], most previous studies did not control for potential confounders of methylation, such as lifestyle [25] and dietary [26-28] factors.

The present study aimed to identify methylation changes in CpG sites associated with recent-onset (RO) psychosis using the Infinium Human Methylation BeadChip (Illumina, San Diego, CA, USA), which covers more than 850 K CpG sites. In the methylation analyses, lifestyle factors such as smoking, alcohol, exercise, and diet were controlled for. Furthermore, we evaluated the association of altered methylation patterns with symptoms and early trauma in patients and healthy controls.

METHODS

Participants

As inclusion criteria, it was required that subjects be aged between 19 and 58 years and meet the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition [29] criteria for SZ spectrum disorders (SZ, schizophreniform disorder, psychotic disorder not otherwise specified [NOS]), brief psychotic disorder, or delusional disorder. RO was defined as duration of illness ≤ 2 years. The exclusion criteria were an intelligence quotient < 70, history of head trauma, serious neurological disorder (epilepsy, stroke, Parkinson’s disease, and/or dementia), and significant medical illness. Participants were recruited from the Korean Early Psychosis Cohort Study (KEPS), a nationwide, multicenter, prospective and naturalistic observational study. Age- and sex-matched healthy individuals were recruited for the control group via advertisements. All participants provided written informed consent in accordance with a protocol approved by the Ethics Committee of Jeonbuk National University Hospital (approval no. CUH 2014-11-002).

Clinical Assessment

The severity of symptoms was evaluated within a week of blood sampling using the Positive and Negative Syndrome Scale (PANSS) [30,31]. Childhood trauma was evaluated with the Early Trauma Inventory Self Report-Short Form (ETISR-SF) [32]. Data on factors related to lifestyle were obtained using the Fagerstrom Test for Nicotine Dependence (FTND) [33], Alcohol Use Disorders Identification Test (AUDIT) [34], Dietary Habits Questionnaire (DHQ) [35], and Physical Activity Rating (PA-R) [36]. For simplicity, we used only the fourth item of the FTND (0, non-smoker; 1, ≤ 10 cigarettes/day; 2, 11−20 cigarettes/day; 3, 21−30 cigarettes/day; and 4, ≥ 31 cigarettes/day) and the mean of items 1 and 2 of the AUDIT ([sum of the scores on items 1 and 2] / 2). The DHQ is a 20-item self-administered questionnaire consisting of three subcategories: there are 5 items on diet regularity, 6 on balanced diet, and 9 on unhealthy diet and eating habits. This scale was developed based on dietary guidance published by the Korean Ministry for Health, Welfare and Family Affairs (2010) [37]. Each item is scored on a three-point scale (1, 3, and 5 points) according to the frequency of the dietary habit. Higher scores indicate better dietary habits in the respective categories. The PA-R is a questionnaire for rating physical activity, with scores ranging from 0 (avoids walking or exercise) to 7 (runs more than 10 miles per week or spends more than 3 hours per week in comparable physical activity).

DNA Methylation Assay

DNA methylation refers to a biochemical process in which a methyl group is bound to the cytosine of a CpG dinucleotide. In turn, this alters the levels of gene transcription, without altering the underlying DNA sequence. Genomic DNA was extracted from peripheral blood mononuclear cell samples using a QIAamp DNA Mini Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s instructions. All prepared samples were bisulfite converted according to the EZ DNA Methylation Kit protocols (Zymo Research, Orange, CA, USA). Bisulfite-converted DNA was assessed using the Infinium MethylationEPIC BeadChip (Illumina), to determine the methylation levels of 853,307 CpG sites. After staining, the BeadChips were scanned, and images were extracted using GenomeStudio software (Illumina). The raw intensity data files (idat files) were preprocessed. First, we excluded 486 probes with detection levels of a p value ≥ 0.01 in 25% of the sample, because they were considered unreliable. Then, we removed 19,377 CpG probes located on the sex chromosome. During this process, 59 probes identified as SNPs were removed. The data were background-corrected and adjusted for dye bias using the methylumi and lumi R packages. To reduce bias in the Infinium I and II assays, the corrected signal value was normalized. The beta-mixture quantile method was used to decompose the density profiles of the Infinium I and II assays into three methylation states: unmethylated (close to 0), partially methylated (close to 0.5), and fully methylated (close to 1). Quantile normalization was then performed to fit each β-distribution of the Infinium II profile to the corresponding β-distribution of the Infinium I profile, leaving 846,055 CpGs for subsequent analysis (Supplementary Fig. 1; available online). Raw data were extracted as beta values ranging between 0 and 1 (β value = methylated intensity / [methylated intensity + unmethylated intensity]). All data analyses, and the visualization of differentially methylated CpGs, were conducted using R software (version 3.3.3; www.r-project.org).

Identification of Differentially Methylated CpG Sites

To identify 2,912 differentially methylated CpG sites associated with RO psychosis, we performed linear regression adjusted for confounding factors: β = α0 + α1 age + α2 sex + α3 education + α4 alcohol + α5 smoking + α6 exercise + α7 diet + ε. To adjust for multiple testing, a false discovery rate (FDR) cutoff p value of < 0.001 was used. Correlations between the methylation values of 2,912 selected CpG sites and total and subscale scores on the PANSS and ETISR-SF were calculated. p values < 0.05 were considered significant.

Functional Enrichment Analysis and Power Calculation

To identify functions associated with the methylation-altered genes (2,535 genes corresponding to 2,912 CpG sites), we performed Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/). We performed KEGG enrichment analysis to identify the biological processes, pathways, and networks shared by these genes. GO terms and KEGG pathways that contained 10 genes and had an Expression Analysis Systematic Explorer score < 0.05 were considered significant. The top 15 GO terms and KEGG pathways were selected.

The study was powered to detect mean differences in DNA methylation between cases and controls at a p value < 3.44E-06, following a previously described method [36]. The power thresholds (ranging between 0 and 1) of the probes, based on the numbers of cases and controls, are shown in Supplementary Fig. 2 (available online). Our study was adequately powered to detect a 0.05% mean difference at ∼40.6% of all sites.

RESULTS

Distribution of Differentially Methylated CpG Sites

Demographic and clinical characteristics of the participants (61 RO psychosis patients and 47 controls) are described in Table 1. Diagnoses included SZ (29.51%), schizophreniform disorder (31.15%), other specified SZ spectrum and psychotic disorders (22.95%), brief psychotic disorder (9.84%), and delusional disorder (6.56%). We identified 2,912 differentially methylated CpG sites associated with RO psychosis (FDR < 0.001) after controlling for age, sex, education level, alcohol use, smoking status, exercise, and diet. Of the 2,912 CpGs, 1,509 sites (51.82%) were hypomethylated, and 1,403 (48.18%) were hypermethylated. Of the 1,509 hypomethylated CpGs, 951 (63.02%), 308 (20.41%), and 11 (0.73%) sites were located in the promoter regions, gene bodies, and 3'-UTRs, respectively (Fig. 1A). When classified according to the CpG contents in the genes, 970 (64.28%), 193 (12.79%), and 49 sites (3.25%) were located in the CpG islands (CGIs), CGI shores, and CGI shelves, respectively (Fig. 1B). Of the 1,403 hypermethylated CpGs, 505 (35.99%), 644 (45.90%), and 51 (3.64%) sites were located in the promoter regions, gene bodies, and 3'-UTRs, respectively (Fig. 1C). In addition, 526 (37.49%), 375 (26.73%), and 148 sites (10.55%) were located in the CGIs, CGI shores, and CGI shelves, respectively (Fig. 1D). The top 20 differentially methylated sites with the lowest FDR values are presented in Table 2. Corresponding gene information of top 20 sites are provided in the online Supplementary Table 1. Ten sites from the ACIN1, PDCD10, SERPINI1, GNAS, HIST1H4J, ZNF668, SSRP1, ANKRD13B, HIST1H4K, CDKN2A, CDKN2B-AS1, and FKBP2 genes, and ten from the ETFA, TYRO3P, RNF31, PSME2, NUMBL, MARCH1, ANP32C, Autophagy-related 4B cysteine peptidase (ATG4B), IGF2BP2, and METTL16 genes, showed that patients were more likely to have a hypo- or hypermethylated status than controls. The full names of these genes, and the top 100 differentially methylated sites, are provided in the online Supplementary Table 2.

Table 1

Demographic and clinical characteristics of participantsa

Variable	Patients (n = 61)	Control (n = 47)	p value
Age	31.820 ± 12.075	30.872 ± 8.882	0.598a
Sex			1.205b
Female	38 (62.2)	34 (72.3)
Male	23 (37.7)	13 (27.7)
Education	13.934 ± 2.658	15.745 ± 0.988	<0.001c
DI (mo)	18.62 ± 41.637
AUDIT	0.984 ± 0.152	0.181 ± 0.152	0.369a
CPZ dose (mg/d) (n = 23)	324.68 ± 220.2
DHQ (n = 54)	60.222 ± 1.715	56.638 ± 1.904	0.164a
FTND (n = 59)	0.373 ± 0.087	0.149 ± 0.061	0.037a
PAR (n = 60)	1.775 ± 0.218	1.660 ± 0.241	0.724a
PANSS
Positive	19.328 ± 5.732
Negative	13.164 ± 6.336
General	34.295 ± 8.405
Total	66.787 ± 15.641
ETISR-SF (n = 58)
Emotional	1.810 ± 1.791	0.489 ± 1.081	0.000a
General	1.810 ± 1.791	0.745 ± 1.073	0.000a
Physical	2.069 ± 1.664	1.149 ± 1.367	0.003a
Sexual	0.397 ± 1.025	0.277 ± 0.615	0.834a
Total	6.138 ± 4.803	2.660 ± 2.656	0.000a

Values are presented as mean ± standard deviation or number (%).

FTND, Fagerstrom Test for Nicotine Dependence; AUDIT, Alcohol Use Disorders Identification Test; DHQ, Diet History Questionnaire; PAR, Physical Activity Ratio; CPZ, chlorpromazine; DI, duration of illness; ETISR-SF, Early Trauma Inventory Self Report-Short Form; PANSS, Positive and Negative Syndrome Scale.

aMann–Whitney U test; bchi-square test; cValues are presented as mean ± standard deviation.

Fig. 1

Percentage of the 2,912 CpG sites according to their location and contents in the genes. Of the 2,912 CpGs, 1,509 sites (51.82%) were hypomethylated and 1,403 sites (48.18%) were hypermethylated. Of the hypomethylated CpGs, (A) 951 (63.02%) 308 (20.41%) and 11 sites (0.73%) were located in the promoter regions, gene bodies and 3'-UTR; (B) 970 (64.23%), 193（12.79%）and 49 sites (3.25%) were located in the CpG islands (CGIs)，CGI shores and CGI shelves. Of the hypermethylated CpGs, (C) 505 (35.99%), 644 (45.90%) and 51 sites (3.64%) were located in the promoter regions, gene bodies and in 3'-UTR; and (D) 526 (37.49%), 375 (26.73%）and 148 sites (10.55%) were located in the CGIs, CGI shores and CGI shelves.

Table 2

Top 20 differentially methylated CpG sites associated with recent onset psychosis

Probe	Chromosome	Position	Illumina gene annotation	Genes within 50 kb of associated CpG	∆b	FC	p value	FDR
cg23712970	14	23540735	ACIN1	PSMB5; PSMB11; CDH24; C14orf119; LMLN2; CEBPE	−0.0387	−1.1941	9.98E-20	8.45E-14
cg18128437	3	167453397	PDCD10; SERPINI1	PDCD10; SERPINI1	−0.0324	−1.2067	1.03E-16	2.91E-11
cg19259030	15	76553926	ETFA;TYRO3P	C15orf27	0.0202	1.0436	7.07E-17	2.91E-11
cg05874478	14	24615768	RNF31;PSME2	PCK2; NRL; DCAF11; FITM1; EMC9; AK307150; PSME1; IRF9; REC8; IPO4; TM9SF1	0.0357	1.0862	6.85E-16	1.45E-10
cg14094693	13	49345282	NA	PSME2P2	0.0310	1.0540	3.22E-15	5.44E-10
cg23496597	20	57463725	GNAS	GNAS-AS1	−0.0301	−1.0615	1.61E-14	2.27E-09
cg26313233	2	114415511	NA	RPL23AP7; RABL2A; SLC35F5	0.0202	1.0418	4.44E-14	5.37E-09
cg26933384	19	41188655	NUMBL	ADCK4; ITPKC	0.0487	1.0627	6.28E-14	6.64E-09
cg25335190	6	27791899	HIST1H4J	HIST1H2BL; HIST1H3H; HIST1H2AJ; HIST1H2BM; HIST1H4K; HIST1H2AK; HIST1H2BN; HIST1H2AL; HIST1H1B; HIST1H3I; HIST1H4L; HIST1H4PS1; HIST1H2AI; HIST1H2BPS2;	−0.0456	−1.4136	9.95E-14	8.42E-09
cg26850290	16	31072493	ZNF668	STX4; BC039500; ZNF646; PRSS53; VKORC1; BCKDK	−0.0214	−1.0293	9.36E-14	8.42E-09
cg02688348	11	57103429	SSRP1	TNKS1BP1; P2RX3; PRG3	−0.0222	−1.2388	1.29E-13	9.17E-09
cg24132580	4	165118938	MARCH1; ANP32C	MARCH1;ANP32C	0.0298	1.0815	1.30E-13	9.17E-09
cg25447461	17	27920715	ANKRD13B	TAOK1; ABHD15; GIT1; CORO6; ABHD15-AS1; SSH2; TP53I13; ANKRD13B	−0.0180	−1.1321	1.56E-13	1.02E-08
cg15028047	2	242607778	ATG4B	THAP4; ING5; DTYMK	0.0184	1.0268	2.16E-13	1.15E-08
cg24376286	2	198245629	NA	SF3B1a	0.0253	1.0552	1.92E-13	1.15E-08
cg25863503	6	27799295	HIST1H4K	HIST1H2BL; HIST1H3H; HIST1H2AJ; HIST1H2BM; HIST1H4J; HIST1H2AK; HIST1H2BN; HIST1H2AL; HIST1H1B; HIST1H3I; HIST1H4L	−0.0344	−1.2981	2.17E-13	1.15E-08
cg05303981	9	21995305	CDKN2A; CDKN2B-AS1	CDKN2A-DT; MTAP; CDKN2B-AS1; CDKN2B	−0.0132	−1.1116	2.34E-13	1.16E-08
cg12394426	11	64009745	FKBP2	STIP1;FERMT3; TRPT1; VEGFB; DNAJC4; PLCB3; BAD; GPR137; KCNK4; NUDT22; DNAJC4; KCNK4-TEX40; PPP1R14B; PPP1R14B-AS1	−0.0264	−1.2046	3.34E-13	1.57E-08
cg24450631	3	185542709	IGF2BP2	IGF2BP2	0.0298	1.0617	8.16E-13	3.14E-08
cg02030802	17	2409437	METTL16	METTL16	0.0171	1.0205	9.68E-13	3.56E-08

Positive fold change indicates hypermethylation in patients compared to controls.

FC, fold change; FDR, false discovery rate; NA, not available; ∆β, group mean β of patients–group mean β of controls.

aGWAS candidate genes.

Functional Enrichment Analysis

The top 15 significantly enriched KEGG pathways included the oxytocin signaling pathway, long-term depression pathway, axon guidance, endometrial cancer, long-term potentiation, MAPK signaling pathway, and glutamatergic pathway, among others (Table 3). In total, 341 genes were annotated with 15 KEGG pathways. The top 15 significantly enriched GO terms are provided in the online Supplementary Table 3.

Table 3

Top 15 KEGG pathways

Pathway ID	Description	Gene count	p value	FDR
hsa04921	Oxytocin signaling pathway	31	8.84E-04	1.16E+00
hsa04730	Long-term depression	16	1.82E-03	2.39E+00
hsa04360	Axon guidance	26	2.89E-03	3.76E+00
hsa05213	Endometrial cancer	13	9.96E-03	1.24E+01
hsa04720	Long-term potentiation	15	1.20E-02	1.48E+01
hsa04010	MAPK signaling pathway	41	1.25E-02	1.53E+01
hsa04724	Glutamatergic synapse	22	1.29E-02	1.58E+01
hsa04520	Adherens junction	15	2.23E-02	2.58E+01
hsa04144	Endocytosis	38	2.43E-02	2.78E+01
hsa04120	Ubiquitin mediated proteolysis	24	2.72E-02	3.06E+01
hsa04510	Focal adhesion	33	2.98E-02	3.30E+01
hsa04270	Vascular smooth muscle contraction	21	3.18E-02	3.48E+01
hsa04370	VEGF signaling pathway	13	3.34E-02	3.62E+01
hsa01100	Metabolic pathways	155	3.34E-02	3.62E+01
hsa04725	Cholinergic synapse	20	3.53E-02	3.78E+01

KEGG, Kyoto Encyclopedia of Genes and Genomes; FDR, false discovery rate.

Relationship between the Beta Values of the Selected CpG Sites and Clinical Indicators

Significant correlations are shown in Figure 2. In the patient group, positive correlations were found between the beta value of cg13562874 and the ETISR-SF total score, while negative correlations were identified between the beta values of cg13810931 and cg18128437 and the PANSS total score. In the control group, the ETISR-SF total score exhibited negative correlations with the beta values of cg23207361 and cg23408615, and a positive correlation with the beta value of cg25548986.

Fig. 2

Correlation between the beta-values of selected CpG sites and clinical indicators. In patient group, (A) positive correlation between the beta-value of cg13562874 near ZSWIM8 gene (p = 0.00097) and total score of the ETISR-SF; and negative correlations (B) between the beta-values of cg13810931 near PITPNM1 gene (p = 0.00053) and (C) cg18128437 near PDCD10 and SERPINI1 genes (p = 0.00074) and total score of the PANSS. In control group, negative correlations (D) between the beta-value of cg23207361 near BAT2 gene (p = 0.0053) and (E) cg23408615 near intergenic region (p = 0.0044), and total score of the ETISR-SF and (F) positive correlation between the beta-value of cg25548986 near ONECUT1 gene (p = 0.00019) and ETISR-SF total score.

DISCUSSION

Epigenetics provides insight into the molecular link between the genetic and environmental factors that contribute to the complex etiology of diseases. The goal of the present MWAS was to identify methylation biomarkers in the peripheral blood of patients with RO psychosis. We found associations of several novel methylated genes with psychosis and early trauma. The results provide new insight into the epigenetic causes and treatment of psychotic disorders.

Analysis of the distribution of differentially methylated CpG sites within genes demonstrated that hypomethylation was more common in the promoter region, and hypermethylation in the gene body. Although the role of DNA methylation in gene bodies remains unclear, gene body DNA methylation in highly expressed genes is often seen in human cells [39]. CGIs in gene bodies reportedly act as alternative promoters [40,41], and tissue- or cell type-specific CGI methylation is prevalent in gene bodies [41,42]. Similarly, in this study, hypomethylation was more common in CGIs, while hypermethylation was more abundant outside CGIs. This indicates that psychosis-associated CpG sites often do not reside in CGIs. Accumulating evidence suggests that phenotypically relevant variation in DNA methylation often occurs outside CGIs [43].

Several of the 22 genes associated with the top 20 differentially methylated CpG sites have been implicated in the pathogenesis of SZ. Two genes, apoptotic chromatin condensation inducer 1 (ACIN1) and programmed cell death 10 (PDCD10), are involved in apoptosis. Our finding that hypomethylation of cg23712970 in ACIN1 and PDCD10 may result in increased expression of these genes is partially consistent with the recent finding that expression levels of the BCL2 and CASP3 apoptosis genes were significantly higher in the peripheral blood lymphocytes of schizophrenic patients [44]. Analysis of mRNA expression in the post-mortem brain tissue of patients with SZ also revealed altered mRNA expression in the IFITM2, TIAL1, PDCD6, and MCL genes, which are involved in programmed cell death [45]. ATG4B is another gene involved in autophagy. Furthermore, the hypomethylation of cg15028047 seen in the present study is in line with the decreased expression of autophagy-related genes in Brodmann area 22 in previously reported in SZ patients [46,47]. Importantly, genome-wide profiling of DNA methylation in induced pluripotent stem cells demonstrated that 23 differentially methylated genes between SZ patients and controls were functionally annotated as ‘Autophagy’ [48]. It should be noted that, as antipsychotic drugs are known to induce autophagy [49], our finding could be related either to disease processes or uncontrolled environmental factors. Hypomethylation of the serpin family I member 1 (SERPINI1), ring finger protein 31 (RNF31), and FKBP prolyl isomerase 2 (FKBP2) genes, which appear to enhance axonal growth, prevent neuroinflammation, and accelerate protein folding and trafficking, respectively, seems counterintuitive given that increased expression of these genes could be beneficial. A possible explanation could be that these methylation changes may reflect either compensatory mechanisms in response to disease processes or the effects of antipsychotic drugs. Hypermethylation of methyltransferase-like 16 (METTL16) is notable, in that a methionine-induced animal model of SZ [50], and the methylation hypothesis of SZ, provide solid support for a role of epigenetics in understanding the pathophysiology of SZ [51]. Histones play a central role in transcription regulation, DNA repair, DNA replication, and chromosomal stability. Our findings regarding the hypomethylation of H4 clustered histone 11 (HIST1H4J) and H4 clustered histone 12 (HIST1H4K) are novel; most previous studies of SZ have considered histone methylation at H3K9 [52,53] and global histone H4 acetylation [54]. Hypermethylation of the membrane-associated ring-CH-type finger 1 (MARCH1) gene, which is involved in regulating vesicular transport and endocytosis, merits further investigation given the novel hypothesis that altered clathrin-mediated endocytosis (CME) is responsible for the pathophysiology of SZ [55]. When we compared our top 100 CpG sites (with 120 genes) to the 172 differentially methylated positions reported by a large-scale methylome study of SZ [18], only four genes, caspase recruitment domain family member 14 (CARD14), nuclear factor of activated T cell 1 (NFATC1), phosphofurin acidic cluster sorting protein 2 (PACS2), and SKI proto-oncogene (SKI), overlapped. Moreover, comparison with the 100 top-ranked SZ-associated genes in the study by Dempster et al. [15] revealed no overlap. This lack of overlap may be attributable to differences in disease stage and confounding factors, such as the study subjects, total number of CpG sites examined, and sample size. Overall, it should be noted that fold- change and △β values for the top 20 CpG sites were low. These low values were also reported in previous studies [18,24]. This suggests that marked differences in methylation status between patients and controls may be difficult to detect in peripheral blood samples, or by microarrays.

Some of the top 15 KEGG pathways have important implications for SZ. Oxytocin plays a central role in birthing, mother-infant bonding, and a wide range of related social behaviors, and is considered a candidate susceptibility gene for SZ. Altered oxytocin receptor (OXTR) methylation and OXTR polymorphisms have been reported in SZ, and in individuals at ultra-high risk for psychosis [22,56]. As for the axon guidance pathway, a recent analysis of European and American GWAS of SZ revealed that 15 axon guidance pathway related genes were associated with the disease [57]. The mitogen-activated protein kinase (MAPK)-associated pathway activates transcription factors related to learning, memory, cell proliferation, and apoptosis. Hypermethylated genes related to MAPK signaling pathways were reported in SZ-affected monozygotic twins [58]. In the endocytosis pathway, 38 genes were enriched, including MARCH1. Interestingly, DNA methylation-mediated modulation of endocytosis has been proposed as a mechanism underlying synaptic function [59].

Higher trauma scores in our patient group were associated with hypermethylation of the zinc finger SWIM- type containing 8 (ZSWIM8) gene. The C2H2-type zinc finger proteins (C2H2-ZNFs) are known to play significant roles in the regulation of neural stem cell activity and subsequent brain development [60,61]. Accumulating evidence suggests a role for C2H2-ZNFs in the development of SZ [62-64]. Therefore, this finding suggests that earlier trauma may decrease the expression of ZNFs, in turn leading to aberrant brain development. However, in the control group, higher trauma scores were associated with hypomethylation of the branched-chain amino acid transaminase (BAT2) gene and the intergenic region, as well as hypermethylation of the one cut homeobox 1 (ONECUT1) gene. The BAT2 and ONECUT1 genes are involved in glucose metabolism and diabetes mellitus. It is difficult to determine the implications of these findings given that the controls were all healthy. Another interesting finding was that, in the patient group, higher PANSS scores were correlated with hypomethylation of the phosphatidylinositol transfer protein membrane associated 1 (PITPNM1), PDCD10, and SERPINI1 genes. PITPNM1 plays a role in maintaining normal diacylglycerol levels in the Golgi apparatus, while the PDCD10 and SERPINI1 genes are involved in apoptosis and axonal growth. Hence, these results suggest that more severe psychopathology may be associated with more active changes in gene expression at both the cellular and axonal levels.

Some limitations should be considered when interpreting the findings of this study. First, the small sample size and diagnostic heterogeneity among the participants may limit the generalizability of the findings. However, it should be noted that the present study was adequately powered to detect a 0.05% mean difference at ∼40.6% of all sites. Second, as some of the patients were on antipsychotics (19 were antipsychotic-naïve, 19 were antipsychotic-free, and 23 were medicated), potential effects of medication on the epigenome cannot be ruled out. Considering several reports on the effects of antipsychotics on methylation profiling in animals [65] and patients with SZ [66], this factor should be controlled for in future studies. Third, significant CpG sites were not validated with independent samples. Fourth, as we did not have RNA samples, the effects of methylation changes on gene expression were not explored. Fifth, cellular heterogeneity of blood samples was not controlled for. Given that cellular composition explains much of the observed variability in DNA methylation [67], the cellular heterogeneity of our samples may have affected the results. Sixth, the SNPs were not all removed. As sequence polymorphisms may affect the readout of methylation signals in Infinium arrays [68], this should be considered in future studies.

As strengths of the present study, the patients were all in a relatively early disease stage, and various lifestyle factors were controlled for. Most importantly, this was the first study to investigate the association between early trauma and DNA methylation changes in RO psychosis. In summary, using a methylome-wide array, we identified differentially methylated CpG sites in patients with RO psychosis compared to controls. The genes associated with the top 20 differentially methylated sites were largely novel, had not been reported in previous methylation studies, and were involved in apoptosis, autophagy, axonal growth, neuroinflammation, protein folding, etc. In the patient group, significant associations of these novel methylated genes with early trauma and psychopathology were observed. These findings support the notion that altered DNA methylation could be involved in early trauma and the pathophysiology of psychosis.

The data set analyzed in this study is available in the Gene Expression Omnibus (GEO) repository (GSE 157252, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE157252).