Region-based analysis of altered DNA methylation in schizophrenia

As DNA methylation is often correlated across adjacent CpG sites [20,21], we used the Illumina Methylation Analyser (IMA) [22] package to aggregate and assess DNA methylation across adjacent probes within annotated regions of the genome. Again, we observe widespread evidence for schizophrenia-associated DNA methylation differences, occurring specifically in the PFC (Figure S1 and Table S6 in Additional file 1). Most notably we observe a large (approximately 8 kb) region spanning the gene body of the Neuritin 1 (NRN1) gene across which 29 adjacent CpG sites are consistently hypomethylated in schizophrenia patients compared with controls (Figure 1B). The mean DNA methylation difference across all 29 sites is -1.77% (FDR = 3.1e-2), with a notable region of hypomethylation (-3.7%) spanning three consecutive DMPs (FDR = 2.4e-2). This is of particular interest because NRN1 plays a well-established role in neurodevelopment and synaptic plasticity [23], and genetic variation in the gene has been linked to cognitive phenotypes in schizophrenia [24].

Validation and replication of schizophrenia-associated DNA methylation differences in the prefrontal cortex

To confirm the array data, bisulfite-pyrosequencing was used to validate schizophrenia-associated DMPs in the vicinity of three genes (NRN1 (cg00565348), C8A (cg23844013), and RASA3 (cg24803255)) in the same samples. Bisulfite-PCR amplification was performed in duplicate using the primers and assay conditions in Table S7 in Additional file 1. Fully methylated and fully unmethylated control samples were included in all experiments. For each amplicon we confirmed significant DNA methylation differences in the same direction as reported by the 450K array (Figure 1C; Figure S2 and Table S8 in Additional file 1). We subsequently generated a replication PFC (BA9) 450K dataset using schizophrenia and control brains archived at the Douglas Bell-Canada Brain Bank, Montreal, Canada (n = 33, 18 schizophrenia and 15 controls; Tables S9 and S10 in Additional file 1) using the Illumina 450K array as described above. Strikingly, DNA methylation differences at the top-ranked PFC DMPs identified in the ‘LBBND’ cohort (listed in Table S11 in Additional file 1) were strongly correlated with schizophrenia-associated differences at the same probes in PFC tissue from the Montreal replication samples (r = 0.54, P = 6.8e-09) with a consistent direction of effect observed across both cohorts and significant differences observed for top-ranked DMPs (Figure 1C,D).

Evidence for schizophrenia-associated gene co-methylated modules in the prefrontal cortex

We next employed weighted gene co-methylation network analysis (WGCNA) to undertake a systems-level view of the DNA methylation differences associated with schizophrenia in the PFC [25,26]. Using probe-wise DNA methylation data we identified 110 modules (representing discrete networks of co-methylated sites), and the first principal component of each individual module (termed the ‘eigengene’) was used to assess the relationship with disease status. Twelve PFC modules were significantly associated with schizophrenia, most notably the ‘black’ (n = 6,647 probes, r=-0.52, P = 0.0004) and ‘pink’ modules (n = 4,399 probes, r=-0.45, P = 0.002) (Figure 2A). Module membership in both modules is strongly correlated with probe-level disease significance (black module: r =0.39, P <1e-200; pink module: r = 0.22, P = 2.4e-49) (Figure 2B). Repeating WGCNA on the cerebellum data yielded no schizophrenia-associated modules, providing further evidence of brain region-specific DNA methylation differences in schizophrenia. Furthermore, a comparison of modules across the two brain regions showed that while a subset of 43 PFC modules correlated strongly (r >0.8) with 51 cerebellum modules, none of the PFC schizophrenia-associated modules were correlated with modules in the cerebellum. Finally, despite the low number of samples and lack of power for replicating WGCNA associations in the Montreal replication dataset, the modular structure of DNA methylation in the PFC was found to be strongly preserved across both brain cohorts (Figure S3 in Additional file 1), with both the black and pink modules being highly conserved and negatively correlated with schizophrenia (black module: r=-0.26; pink module: r=-0.10).

Open in a separate window

Figure 2

Schizophrenia-associated gene co-methylation modules in the prefrontal cortex. (A) Heatmap representing correlations between module eigenvalues (ME) and clinical/demographic variables. Each row represents a specific module, indicated by an arbitrary module color labeled on the y-axis. Dark red indicates strong positive correlation, dark green indicates strong negative correlation, and white indicates no correlation, as indicated in color scale bar. P-values are given in parentheses. The black (P = 0.0004) and pink (P = 0.002) WGCNA modules are most significantly associated with schizophrenia in the PFC. (B) Module membership in both modules is strongly correlated with probe significance with disease (black module: r = 0.39, P <1e-200; pink module: r = 0.22, P = 2.4e-49). (C) Ingenuity Pathway Analysis of genes associated with CpG sites in the black module reveals a highly significant enrichment of disease pathways related to schizophrenia and other neuropsychiatric disease. (D) Co-methylation networks between the top 1% of loci in the black and pink modules ranked by their module membership, where thicker lines indicate more connected genes. Black edges represent co-methylation between probes in the black module and pink edges represent co-methylation in the pink module; thicker and deeper colored edges indicate stronger correlations between probes. Interestingly, the two modules are connected by SHANK2, which encodes a molecular scaffold protein that plays a role in synaptogenesis. Furthermore, several of the genes are associated with schizophrenia in recent collaborative genome-wide association studies (GWAS) [27]. GWAS significance for SNPs in the vicinity of each locus is depicted by color.

Schizophrenia-associated co-methylated modules are enriched for neurodevelopmental pathways and loci previously implicated in the disorder

To test whether the identified schizophrenia-associated modules were biologically meaningful, pathway and gene ontology analyses were performed using Ingenuity Pathway Analysis (IPA) and EASE (DAVID) software [28]. Enrichment for specific pathways and biological functions was determined relative to the relevant microarray gene list using a right-tailed Fisher’s test. IPA of genes associated with CpG sites in the top-ranked PFC modules reveals a highly significant enrichment of disease pathways related to ‘schizophrenia and other neuropsychiatric disorders’ in the black module (P = 3.67e-15 to 1.27e-03; Figure 2C), and ‘nervous system development and function’ in both the black (P = 3.80e-18 to 1.37e-03; Figure S4 in Additional file 1) and pink (P = 4.04e-15 to 3.16e-03; Figure S5 in Additional file 1) modules. Consistent with the IPA results, gene ontology analysis using EASE (DAVID) [28] shows that both modules are significantly enriched for functions related to nervous system and neuron development, synaptic transmission and calcium ion binding (Table S12 in Additional file 1). An analysis of co-methylation networks between the top-ranked ‘hub’ probes in the black and pink modules shows that they are connected by DNA methylation in the vicinity of SHANK2, a molecular scaffold protein that plays a critical role in neurodevelopment and synaptogenesis, and has been widely implicated in several neurodevelopmental disorders, including schizophrenia [29] (Figure 2D). Furthermore, several of the probes are located in genes that have been strongly implicated in schizophrenia by a recent large collaborative genome-wide association study (GWAS; for example, CACNA1C, TSNARE1, PITPNM2, and HLA-L) [27].

Genome-wide quantification of DNA methylation

Bisulfite-pyrosequencing

Independent verification analysis was performed on three PFC schizophrenia-associated DMPs (in the vicinity of NRN1, C8A, and RASA3), based on results from both the probe-wise and region-level analysis. Pyrosequencing assays were designed using the PyroMark Assay design software (Qiagen, Hilden, Germany). Bisulfite-PCR amplification was performed in duplicate using the primers and assay conditions in Table S7 in Additional file 1. Fully methylated and fully unmethylated control samples were included in all experiments.

Replication analysis

DNA samples from the Montreal cohort were analyzed using the Illumina Infinium HumanMethylation450 BeadChip using the same processing and quality control steps as described above. For the top 100 schizophrenia-associated DMPs in the LBBND cohort (Table S4 in Additional file 1), we extracted β-values for matched probes in the Montreal dataset. Pearson’s correlation was used to assess the relationship between disease-associated DNA methylation differences in the two datasets at these probes.

Weighted gene co-methylation network analysis

Network analysis was performed on normalized probewise DNA methylation data (PFC: n = 445,617 probes; cerebellum: n = 440,836 probes) using WGCNA [26]. For each brain region pair-wise correlations were used to identify modules of highly co-methylated probes, independent of disease status. Specifically, an unsigned network was created with a soft threshold parameter of 6 using the blockwiseModules function based on a block size of 20,000. Each module was then labeled with a unique color name. The first principle component of the methylation matrix of each module was calculated to give a ‘module eigengene’ (ME) for each module (a weighted average methylation profile). To identify modules associated with schizophrenia, each of the ME values was regressed on disease status, in addition to other sample-specific variables: pH, age, gender, brain weight and cerebellum weight. To test whether the identified schizophrenia-associated modules were biologically meaningful, pathway and gene ontology analyses were performed using IPA [44] and EASE (DAVID) [28,45]. Correlations were performed between the ME and DNA methylation values to calculate the module membership of each probe (the extent to which a given probe contributes to the ME). This module membership value was used to identify the ‘hub’ genes in each module. The function modulePreservation in the WGCNA package was used to calculate module preservation statistics for the Montreal PFC data based on the PFC modules built from the discovery dataset.

This work was supported by grants from the UK Medical Research Council (MRC) (grant number MR/K013807/1) and US National Institutes of Health (grant number AG036039) to JM. RP and HS were funded by MRC PhD studentships. The human embryonic and fetal material was provided by the Joint MRC (grant number G0700089)/Wellcome Trust (grant number GR082557) Human Developmental Biology Resource.