Association of DNA Methylation Differences With Schizophrenia in an Epigenome-Wide Association Study

Importance  DNA methylation may play an important role in schizophrenia (SZ), either directly as a mechanism of pathogenesis or as a biomarker of risk.

Objective  To scan genome-wide DNA methylation data to identify differentially methylated CpGs between SZ cases and controls.

Design, Setting, and Participants  Epigenome-wide association study begun in 2008 using DNA methylation levels of 456 513 CpG loci measured on the Infinium HumanMethylation450 array (Illumina) in a consortium of case-control studies for initial discovery and in an independent replication set. Primary analyses used general linear regression, adjusting for age, sex, race/ethnicity, smoking, batch, and cell type heterogeneity. The discovery set contained 689 SZ cases and 645 controls (n = 1334), from 3 multisite consortia: the Consortium on the Genetics of Endophenotypes in Schizophrenia, the Project among African-Americans To Explore Risks for Schizophrenia, and the Multiplex Multigenerational Family Study of Schizophrenia. The replication set contained 247 SZ cases and 250 controls (n = 497) from the Genomic Psychiatry Cohort.

Main Outcomes and Measures  Identification of differentially methylated positions across the genome in SZ cases compared with controls.

Results  Of the 689 case participants in the discovery set, 477 (69%) were men and 258 (37%) were non–African American; of the 645 controls, 273 (42%) were men and 419 (65%) were non–African American. In our replication set, cases/controls were 76% male and 100% non–African American. We identified SZ-associated methylation differences at 923 CpGs in the discovery set (false discovery rate, <0.2). Of these, 625 showed changes in the same direction including 172 with P < .05 in the replication set. Some replicated differentially methylated positions are located in a top-ranked SZ region from genome-wide association study analyses.

Conclusions and Relevance  This analysis identified 172 replicated new associations with SZ after careful correction for cell type heterogeneity and other potential confounders. The overlap with previous genome-wide association study data can provide potential insights into the functional relevance of genetic signals for SZ.

Schizophrenia (SZ) is a neurodevelopmental disorder affecting 1% of the population worldwide, often resulting in lifelong disability.¹ Treatments are not efficacious for the most disabling symptoms such as deficits in working memory, attention, and social cognition.² Understanding the underlying biology of SZ can inform both treatment and prevention strategies. Heritability estimates are high (70%-81%) and some genetic variants have now been identified³ including multiple common variants,^4,5 copy number variants,⁶ and rare mutations.^7,8 The Psychiatric Genomics Consortium has identified 108 loci⁹; despite its size, this study has highlighted modest effect sizes and the importance of regulatory DNA.¹⁰ Mechanistic studies in transgenic animals have significant limitations because animal models have different brain structures and cannot fully recapitulate the phenotypic components of SZ.¹¹

Discordance rates between monozygotic twins are 48% to 50%, reflecting a complex genetic architecture and a role for nongenetic factors in the etiology of SZ.¹² Environmental risk factors include growing up in lower socioeconomic and urban environments,¹³ prenatal infections,¹⁴ obstetric complications,¹⁵ and maternal malnutrition.¹⁶ Some of these environmental factors, such as obstetric complications, can interact with genetic risk factors to alter risk or worsen the course of SZ.¹⁷

Epigenetic processes, which can be modified by environment but also have a genetic component and regulate genetic signals, have been postulated as links between environmental exposures, genetic risk, and SZ. Epigenetics refers to heritable modifications of DNA or associated factors that have information content beyond the primary DNA sequence. DNA methylation (DNAm), a stable modification of cytosines in CpG dinucleotides, regulates neurobiological processes and is implicated in disorders such as Rett syndrome.¹⁸ DNA methylation depends on dietary sources of the essential amino acid methionine, and periods of famine have doubled SZ risk.^19,20 The later onset of SZ, its episodic nature, and the induction of psychotic episodes while undergoing treatment with l-methionine²¹ also suggest a role for epigenetics.

Several case-control studies using blood^22-26 and postmortem brain tissues^27-30 have revealed differential genome-wide methylation patterns, but the number of CpGs interrogated and sample sizes have been limited. To our knowledge, only 1 large-scale methylome study of SZ of equivalent size to ours has been published. It used a methylated DNA enrichment method (methyl-CpG binding domain protein-enriched genome sequencing) in a Swedish population³¹ but did not account for cell heterogeneity in blood, which can lead to false-positive results.^32,33 Only 1 of the smaller studies incorporated this critical step in their analysis.²⁶

In this study, we described a comprehensive analysis of genome-scale DNAm in a SZ case-control consortium and its relationship to genetic variation. We profiled genome-wide DNAm in 1334 blood samples using the Illumina HumanMethylation450 BeadChip, accounting for cell type heterogeneity to identify DNAm associated with disease. High-confidence findings were replicated in an independent set of 497 individuals. We established genotypic-methylation relationships and integrated these results with published genome-wide association study data. This study supports the importance of integrating DNAm and quantitative trait analyses into population-based SZ genomic analysis.

Box Section Ref ID

eTable 1 in the Supplement describes the demographic characteristics of the discovery and replication samples. The discovery set consisted of 1334 blood-derived DNA samples (689 SZ cases and 645 controls) from a population of probands and controls derived from the collaborative efforts of 3 ongoing, National Institute of Mental Health–funded, multisite consortia: the Consortium on the Genetics of Endophenotypes in Schizophrenia,³⁴ the Project Among African-Americans to Explore Risks for Schizophrenia,³⁵ and the Multiplex Multigenerational Family Study of Schizophrenia.³⁶ In the replication stage, 497 independent samples (247 SZ cases and 250 controls) from the Genomic Psychiatry Cohort³⁷ were frequency matched with respect to smoking status, age, and sex. Written informed consent was obtained from all individuals, and the study was approved by the Johns Hopkins Medicine Office of Human Subjects Research institutional review board and the Consortium on the Genetics of Endophenotypes in Schizophrenia, the Project Among African-Americans to Explore Risks for Schizophrenia, the Multiplex Multigenerational Family Study of Schizophrenia, and the Genomic Psychiatry Cohort.

Eligible probands were at least 18 years of age and received SZ diagnosis according to the DSM-IV. Participants from the Project Among African Americans To Explore Risks for Schizophrenia self-identified as African American and had family members available for assessment. Diagnostic assessment was performed using the Diagnostic Interview for Genetic Studies along with medical records. The Multiplex Multigenerational Family Study of Schizophrenia probands were of European American descent from multigenerational families with at least 1 affected family member with SZ or schizoaffective disorder. The Consortium on the Genetics of Endophenotypes in Schizophrenia participants had at least 1 unaffected full sibling. Self-identified race/ethnicity was later confirmed with genetic information derived from the Infinium HumanMethylation450 BeadChip (Illumina) (eFigure 1 in the Supplement).

Experimental flow is presented in Figure 1. All analyses were performed in R version 3.1.1 (R Statistics), and raw intensity files were preprocessed and quantile normalized together using the Bioconductor package Minfi, version 1.8.9 (Bioconductor).³⁸ Probes on the sex chromosomes were removed along with probes with annotated single-nucleotide polymorphisms (SNPs) (via the Single Nucleotide Polymorphism database 137 [National Center for Biotechnology Information]) in the CpG site (n = 16 817) or at single base extension sites (n = 7781), which left a total of 456 513 autosomal probes for calculation of beta scale methylation ratios and association analysis. Of the 2891 unique samples in the discovery set (including family members), subsequent samples were dropped according to the following quality control measures: (1) median methylated and median unmethylated signal average less than 10.5 (n = 9); (2) separate clustering in the first principal component of genome-wide DNAm (n = 8); (3) discrepancy between self-reported sex and predicted sex from DNAm and genotypes (n = 18); (4) discrepancy between self-reported race/ethnicity and predicted ancestry from DNAm data (n = 0); and (5) missing age, smoking, or phenotype in the database (n = 84). After removing family members, 1334 samples were left as the discovery arm of the case-control study. For the replication set, the same criteria were applied (1 dropped after principal component analysis and 2 were incongruent for race/ethnicity prediction). The resulting 497 samples constituted the replication set.

To identify differentially methylated positions (DMPs) associated with SZ, we fit a linear regression model using the limma Bioconductor package (Bioconductor).³⁹ We tested the association between proportion methylation values (Illumina “Beta” scale) and SZ diagnosis status at each CpG site, adjusting for sex, age, race/ethnicity (African American vs non–African American), current smoking status, estimated cell proportions for the 6 cell types, and the first 2 principal components of the negative control probes. To adjust for multiple testing, a false discovery rate (FDR) cutoff of less than 0.2 was used for genome-wide significance. We used all of the samples and adjusted for race/ethnicity to improve the power of our study, but also performed 3 additional race/ethnicity- and sex-stratified analyses (African American and non–African American only and male only) for clarification of findings.

The eMethods in the Supplement include a detailed description of genome-wide measurement of DNA methylation, assessment of cell heterogeneity and potential confounders, replication analysis, gene ontology and pathway analyses, Illumina genome-wide genotyping and meQTL analysis, and regression modeling formulas.

In epigenetic studies using DNA from heterogeneous tissue, such as blood, it is critical to account for potential confounding by cell type. Each cell type has distinct methylation signatures,⁴⁰ and differences in relative abundance of these cell populations between SZ cases and controls can create associations that are not causal. We implemented a statistical algorithm that estimates cell proportions based on reference DNAm signatures obtained from cell-specific data.^32-37,41 The calculated cell proportions for CD4+ and CD8+ T cells, B cells, monocytes, granulocytes, and natural killer cells for the discovery and replication sets are presented in eFigure 2 in the Supplement. The largest statistically significant differences in cell proportions between cases and controls were observed in CD8+ T cells and monocytes, with SZ showing decreased CD8+ T cells and increased monocytes (CD8+ T cells discovery P value = 3.80E−9, replication P value = 2.95E−5 and monocyte discovery P value = 3.80E−3, replication P value = 8.73E−4). Granulocyte proportions were also higher among SZ cases in the discovery and replication sets (discovery P value = .11 and replication P value = 6.76E−4). CD4+ T cells, natural killer cells, and B cells did not show statistically significant differences in the same direction in both discovery and replication sets.

To our knowledge, only 1 smaller epigenome-wide association study of SZ (N = 105) has accounted for cell heterogeneity.²⁶ Our findings are consistent with the magnitude and direction of their estimated cell proportion differences. We observed no significant association between antipsychotic use and proportions of CD8+ T-cell counts (P value = .62), monocytes (P value = .58), or granulocytes (P value = .54, eFigure 3 in the Supplement).

Our discovery sample included 1334 unrelated individuals (689 SZ cases and 645 controls) from 3 large NIMH consortia: the Consortium on the Genetics of Endophenotypes in Schizophrenia, the Project Among African-Americans to Explore Risks for Schizophrenia, and the Multiplex Multigenerational Family Study of Schizophrenia (see eMethods and eTable 1 in the Supplement for demographic information). Of the 456 513 CpGs examined, 945 CpGs were associated with case-control status at an FDR of less than 0.20.

We found it critical to account for cell type proportions, even after adjusting for sex, age, race/ethnicity, smoking, and batch effects (removing the first 2 principal components of negative control probes). Cell type correction markedly reduced the number of significant results (eFigure 4 and eFigure 5 in the Supplement). To assess residual confounding after adjusting for cell proportion estimates, we plotted adjusted estimates and residual values against cell proportions and showed that the linear models provided an adequate correction for these factors (eFigure 6 and eFigure 7 in the Supplement). In addition to controlling for smoking status, we further removed 294 CpGs directly associated with smoking.^42,43 This eliminated 22 CpGs from the discovery results, leaving 923 SZ-associated DMPs (eTable 2 in the Supplement).

We sought to confirm our results in a replication set of 497 independent samples (247 SZ cases and 250 controls). Of the 923 DMPs from the discovery set, 625 (68%) showed association in the same direction and 172 (19%) of these had a replication P value of less than .05. Quadrant correspondence (proportion of replication CpG effect estimates in the same direction as the discovery estimate) was even stronger when considering more stringent FDR thresholds among the discovery sample, with 80% correspondence at an FDR of less than 0.05 (Figure 2A).

The location, effect size (magnitude of covariate-adjusted methylation change comparing SZ cases with controls), and statistical support for each of the 172 replicated SZ-DMPs are reported in eTable 3 in the Supplement. The first 15 DMPs with the lowest FDR values are presented in the Table. Among the genes near DMPs, those previously associated with SZ include RPS6KA1,⁴⁴MAD1L1,⁵KLF13,⁴⁵SULT4A1,⁴⁶NPDC1,⁴⁷AUTS2,⁴⁸PIK3R1,⁴⁹ and PNPO.⁵⁰ We found genes previously reported to play a role in neuronal development and function (S100A2, SULT4A1, HES1, NTM, MARK4, KIFC3, NPDC1, SCAP, DHCR24, HS6ST1, GAS7, CARD10, DYRK2, RPS6KA4, SLC38A10, CUTA, FOXO1, LRFN3, and IRS1) and genes involved in both T-cell development and neuronal-cell fate (ZC3H12D, TCF3, MAP3K1, and IKZF4). Genes with multiple DMPs nearby included 2 genes previously associated with SZ (DDR1⁵¹ and IFITM1⁵²), a component of the mechanistic target of rapamycin signaling pathway that facilitates synaptic plasticity (RPTOR⁵³), and a gene critical in neurogenesis (NCOR2,⁵⁴Figure 2B and 2C). Genes with SZ-DMPs near them have the same average length as genes on the whole array, but do have more probes per gene.

In male-only analyses, 23 replicated SZ-DMPs were identified (eTable 4 in the Supplement), 16 of which were among the replicated DMPs in the sex-combined analyses. Stratified analysis in non–African American samples revealed similar results, although with less statistical significance given the reduction in sample size (eTable 5 in the Supplement). The African American group showed a greater number of statistically significant DMPs, despite the lower sample size, potentially owing to uncorrected population stratification among samples of African descent because ethnicity adjustment in these data was categorized as African American vs non–African American (eTable 6 in the Supplement). Preliminary post hoc analysis of the influence of antipsychotic medications on DNAm revealed that atypical antipsychotics exerted the biggest effect. We then grouped the SZ cases into 2 subgroups: 1 that received atypical antipsychotics, and 1 that received typical antipsychotics only or no treatment. Four of the 172 SZ-DMPs were affected by atypical antipsychotic treatment (P value ≤ .05), with 1 probe inside of the gene MARCH11 and 3 near the promoters of the genes SATB1 and IFITM1 (eFigure 8 in the Supplement).

Statistically significant (FDR < 0.2) functional pathway categories among SZ-DMPs included those related to regulation of signaling and signal transduction (eTable 7 in the Supplement). While gene ontology is centered on genes, network analysis focuses on their physical interactions. Consistent with the gene ontology results, the top 3 most relevant regulatory networks from our pathway analysis included (1) amino acid metabolism, energy production, and post-translational modification; (2) hereditary disorder, neurological disease, and psychological disorders (eFigure 9 in the Supplement); and (3) neurological disease, organismal injury and abnormalities, and embryonic development.

Our replicated SZ-DMPs overlapped 3 genetic regions identified in a 2014 Psychiatric Genomics Consortium mega-analysis⁹: regions 7, 11, and 37 from that publication list. Enrichment analysis using a Fisher exact test was not statistically significant (P = .47); however, the 450K array covers only 87 of the 105 autosomal PGC regions. In region 11 (chromosome 15), 2 SZ-SNPs (rs7177338 and rs2071382) are associated with methylation at the SZ-DMP (cg25647583), located in an intergenic region and equidistantly located (approximately 500 base pairs) between FURIN and FES genes (Figure 3A-C and eTable 8 in the Supplement). Notably, we observed the same genetic risk association to these SNPs in our data (Figure 3D). We further assessed the potential for DNAm mediation of this genetic risk. We estimated the effect of the SNP genotype G on SZ (SZ affected by G), and the effect of the SNP on SZ, adjusting for methylation (M) at the CpG locus (SZ affected by G + M). Conditioning on M did not significantly reduce the effect of G on SZ (β coefficient estimate for G alone, −0.301; P = .04; β coefficient estimate for G conditioning on M, −0.280; P = .06; Figure 3E), suggesting that the genotype is controlling M and SZ independently.

To our knowledge, our study is the largest epigenome-wide association study of SZ to date, using methods robust to cell-type confounding and with supporting genetic data. We identified and replicated 172 SZ-associated DMPs in 689 cases and 645 controls. The small but significant changes in DNAm are of similar effect size to other methylome studies of complex disorders^33,55,56 and highlight the potential usefulness of using blood as a cost-effective surrogate tissue.

The SZ-DMPs discovered in this study include genes important in neuronal function (eg, NCOR2, SULT4A1, and HES1), genes previously associated with SZ (eg, RPS6KA1, MAD1L1, and DDR1), and genes also involved in T-cell development (ZC3H12D, TCF3, and IKZF4). Several of the neuronal genes have been implicated in Alzheimer disease pathology (SORL1,⁵⁷HES1,⁵⁸MARK4,⁵⁹SCAP,⁶⁰DHCR24,⁶¹GAS7,⁶²CARD10,⁶³CDC37,⁶⁴CUTA,⁶⁵IRS1,⁶⁶ and FOXO1⁶⁷), supporting a common network of dysregulation in SZ and Alzheimer disease.⁶⁸

Schizophrenia cases exhibited statistically significant lower estimated CD8+ T cells and higher monocyte and granulocyte counts than the controls, and these cell proportion differences constitute an important confounder of blood epigenetic studies. The direction of bias can be positive (as is typically presumed), inflating association signals, or negative, reducing ability to see true associations. Cell-specific studies using sorted purified samples are ideal but not always practical in epidemiologic settings, and this adjustment method can alleviate some bias associated with this limitation.

To understand the genetic architecture underlying our results, we integrated methylome-wide data with genotyping information from our samples. The Psychiatric Genomics Consortium association results revealed SZ-SNPs controlling the methylation level of an SZ-DMP (cg25647583), where the risk allele was associated with lower methylation levels in SZ cases. Although our genotyping data set was relatively small (n = 374), we confirmed findings from the large psychiatric consortium meta-analysis.

We attempted to replicate previously published methylome-wide SZ association results. One study reported 450K platform findings among 2 small Japanese study samples (n = 105 discovery and n = 48 replication).²⁶ Although this study adjusted for cell proportions, the samples were not adequately randomized by diagnosis on the arrays, creating batch effects. Attempting to remove this bias using aggressive batch effect correction methods would render the studies noncomparable. We also attempted to replicate the 2014 large (n = 1497) case-control study of SZ samples from Sweden,³¹ but the samples were assayed using methyl-CpG binding domain protein-enriched genome sequencing, a different technology from the one we used in this study, containing largely nonoverlapping CpGs. Only 1 of their 15 replicated methylation blocks overlapped probes available on the 450K array. For their top gene, FAM63B, the closest 450K CpG (cg21149266) lies 126 bp away from their closest validated position (chr15: 59146756). We observed no statistically significant difference between SZ cases and controls at this position in our discovery set (P value = .12) but detected a difference in the same reported direction in our replication set (P value = .002) (eFigure 10 in the Supplement).

In this emerging area of SZ epigenetics, there is bound to be an initial period of valid and interesting differences of opinion and gaps in our knowledge. Moreover, in the case of SZ, the brain is inaccessible, making blood studies valuable. While the effect sizes are small, that is also the case with genetic variation.⁶⁹ Other issues will require clarification including the relative causal and risk weight of factors associated with SZ (diet, medications, and stress) and how these factors increase risk or whether they are epiphenomena. The more causal interpretation is supported by the well-known, biologically based brain developmental effects of environmental factors such as the Dutch Hunger Winter and influenza on the developing brains of embryos. A single study cannot answer all of these challenging questions. We will be looking at the relationship of the blood results reported in this study and related brain and neurocognitive data in subsequent studies. Future a priori planned Schizophrenia Psychiatric Genome-Wide Association Study and sequencing of these Consortia cohorts and integration with epigenetic analysis in these same participants will add to the value of both approaches. This study will hopefully be a foundational element raising as many questions as it resolves in an exciting but novel area of SZ research.

We identified novel differential methylation loci between SZ and control individuals from a large, diverse population, using a widely used and cost-effective, high-throughput technology. Using a robust analysis pipeline, with cell proportion correction and sequential replication in independent data sets, we provided candidate loci that contribute to the knowledge of SZ and possibly other psychoses. These findings lay a firm groundwork for the implementation of genome-wide methylation approaches to the understanding of neuropsychiatric disorders.

Corresponding Author: Andrew P. Feinberg, MD, MPH, Johns Hopkins University Schools of Medicine, Engineering, and Public Health, 855 N Wolfe St, Rangos 570, Baltimore, MD 21205 (afeinberg@jhu.edu).

Submitted for Publication: August 9, 2015; final revision received January 20, 2016; accepted January 22, 2016.

Published Online: April 13, 2016. doi:10.1001/jamapsychiatry.2016.0144.

Author Contributions: Ms Montano and Dr Taub had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. Dr A. P. Feinberg and Dr Fallin cosupervised this work and contributed equally as principal investigators.

Study concept and design: Fallin, A. P. Feinberg.

Acquisition, analysis, or interpretation of data: Montano, Taub, Jaffe, Briem, J. I. Feinberg, Trygvadottir, Idrizi, Runarsson, Berndsen, R. C. Gur, Perry, Sabunciyan, Yolken, Hyde, Kleinman, Sobell, C. N. Pato, M. T. Pato, Go, Nimgaonkar, Weinberger, Braff, R. E. Gur, Fallin, A. P. Feinberg.

Drafting of the manuscript: Montano, Fallin, A. P. Feinberg.

Critical revision of the manuscript for important intellectual content: Montano, Taub, Jaffe, Moore, Perry, Sabunciyan, Hyde, Kleinman, Sobell, C. N. Pato, M. T. Pato, Go, Nimgaonkar, Weinberger, Braff, R. E. Gur, Fallin, A. P. Feinberg.

Statistical analysis: Montano, Taub, Jaffe.

Obtained funding: R. C. Gur, Nimgaonkar, Perry, Braff, A. P. Feinberg.

Administrative, technical, or material support: Briem, J. I. Feinberg, Trygvadottir, Idrizi, Runarsson, Berndsen, Fugman, Hyde, Kleinman, C. N. Pato, M. T. Pato, Sobell.

Study supervision: Fallin, A. P. Feinberg.

Conflict of Interest Disclosures: None reported.

Funding/Support: This work was supported by the National Institutes of Health grant “A Genome-Wide Methylation Scan for Epigenetic Contributions to Schizophrenia” under the U01 cooperative agreement mechanism: grants 5U01MH085270 (Dr A. P. Feinberg), 5U01MH085264 (Dr R. E. Gur), 5U01MH085269 (Dr Nimgaonkar), 5U01MH085281 (Dr Perry), 5U01MH085265 (Dr Braff), and 5U24MH068457-13 (Jay Tischfield, Rutgers University).

Role of the Funder/Sponsor: The sponsors had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

Additional Contributions: We thank Roxann Ashworth, MHS (Genetic Resources Core Facility, Johns Hopkins Institute of Genetic Medicine), for her expertise on genotyping. She was compensated for her work.