Genome-wide DNA methylomic differences between dorsolateral prefrontal and temporal pole cortices of bipolar disorder.

Dorsolateral prefrontal cortex (DLPFC) and temporal pole (TP) are brain regions that display abnormalities in bipolar disorder (BD) patients. DNA methylation - an epigenetic mechanism both heritable and sensitive to the environment - may be involved in the pathophysiology of BD. To study BD-associated DNA methylomic differences in these brain regions, we extracted genomic DNA from the postmortem tissues of Brodmann Area (BA) 9 (DLPFC) and BA38 (TP) gray matter from 20 BD, ten major depression (MDD), and ten control age-and-sex-matched subjects. Genome-wide methylation levels were measured using the 850 K Illumina MethylationEPIC BeadChip. We detected striking differences between cortical regions, with greater numbers of between-brain-region differentially methylated positions (DMPs; i.e., CpG sites) in all groups, most pronounced in the BD group, and with substantial overlap across groups. The genes of DMPs common to both BD and MDD (hypothetically associated with their common features such as depression) and those distinct to BD (hypothetically associated with BD-specific features such as mania) were enriched in pathways involved in neurodevelopment including axon guidance. Pathways enriched only in the BD-MDD shared list pointed to GABAergic dysregulation, while those enriched in the BD-only list suggested glutamatergic dysregulation and greater impact on synaptogenesis and synaptic plasticity. We further detected group-specific between-brain-region gene expression differences in ODC1, CALY, GALNT2, and GABRD, which contained significant between-brain-region DMPs. In each brain region, no significant DMPs or differentially methylated regions (DMRs) were found between diagnostic groups. In summary, the methylation differences between DLPFC and TP may provide molecular targets for further investigations of genetic and environmental vulnerabilities associated with both unique and common features of various mood disorders and suggest directions of future development of individualized treatment strategies.

Introduction

Bipolar disorder (BD) is characterized by recurrent episodes of elevated mood and depression, interspersed with normal mood periods, often accompanied by drastic changes in energy levels that can severely affect the individual’s daily life. While the etiology of BD remains uncertain, its high heritability (Craddock and Sklar, 2013) supports the involvement of genetic and heritable epigenetic factors. DNA methylation is an epigenetic mechanism which can modulate gene expression via the addition of a methyl group at the C5 position of cytosine, mainly in cytosine-phosphate-guanine (CpG) dinucleotides. Since the epigenetic mechanism of DNA methylation is heritable and sensitive to environmental influence, it has received much attention in mood disorder studies (Teroganova et al., 2016).

Previous DNA methylation studies in human brain tissue showed altered expression levels of DNA methylation enzymes (Veldic et al., 2005) and differential methylation in various genes in BD patients (Ludwig and Dwivedi, 2016). Most of these studies focused on the prefrontal cortex (Abdolmaleky et al., 2006; Ghadirivasfi et al., 2011; Kaminsky et al., 2012), with the dorsolateral prefrontal cortex (DLPFC) having important functions in decision making, cognitive and executive functions as well as emotion regulation (Clark and Sahakian, 2008). DLPFC (Brodmann Area [BA] 9 and 46) is functionally connected to the orbital and medial prefrontal networks (Price and Drevets, 2010) and to other cortical regions, particularly various parietal and temporal cortices (Jung et al., 2017). Multiple studies reported structural, functional and neurochemical abnormalities in the DLPFC of BD patients (Birur et al., 2017), including reduced gray matter volume and cortical thickness (Hibar et al., 2017; Price and Drevets, 2010), reduced neuronal and glial density (Rajkowska et al., 2001), reduced resting-state functional connectivity with medial prefrontal cortex (Chai et al., 2011), and reduced activation during executive function tasks (Townsend et al., 2010).

Among the brain regions altered in mood disorders, but remaining unexplored in DNA methylation studies, is the temporal pole cortex (TP; BA38) which generally refers to the rostral end of the temporal lobe. It is highly connected to the amygdala, hippocampus, orbitofrontal cortex, and auditory- and visual-related cortices (Blaizot et al., 2010; Olson et al., 2007). It has been associated with audio-visual information integration (Ohki et al., 2016), semantic memory, fluency and development (Drane et al., 2009; Monzalvo and Dehaene-Lambertz, 2013), recognition of emotions and empathy (Jimura et al., 2009; Pehrs et al., 2017), and social interaction (McGettigan et al., 2017). TP was among the brain regions reported to have reduced cortical thickness in BD and major depressive disorder (MDD) patients (Hibar et al., 2017; van Tol et al., 2014). Gray matter neurite density in TP was significantly reduced in schizophrenia (SCZ) subjects while BD subjects demonstrated a similar trend (Nazeri et al., 2017). BD subjects showed a longer extension of the left arcuate fasciculus towards the TP than controls, suggesting neurodevelopmental abnormalities in this fiber tract could be associated with BD etiology (Sun et al., 2017).

Since BD patients showed structural abnormalities in both DLPFC and TP, we hypothesized that DNA methylomic differences in these brain regions might play a role in BD. To explore whether TP shows DNA methylomic differences (at individual CpG sites and across genomic regions) in BD similar to those of DLPFC, we investigated the genome-wide DNA methylomic differences associated with BD within and between the gray matter of DLPFC (BA9) and TP (BA38). We also included MDD cases to identify BD-specific brain DNA methylomic differences.

Methods

Subjects

All procedures performed were in accordance with the Declaration of Helsinki and had been approved by the Institutional Review Boards of the University of Mississippi Medical Center, Jackson, MS, and the University Hospitals Cleveland Medical Center, Cleveland, OH. Legally-defined next-of-kin provided informed consent for the collection of tissue, medical records, and interviews. Structured Clinical Interview for DSM-IV Axis I Disorders (First et al., 2002) was administered by a Master-level social worker to knowledgeable informants of the subjects as previously described (Mahajan et al., 2018). Diagnosis of depression based on next-of-kin interview has demonstrated a high degree of agreement with diagnoses made by interviewing patients (Dejong and Overholser, 2009). To determine subjects’ psychopathology, a board-certified clinical psychologist and a board-certified psychiatrist independently reviewed the diagnostic interview scoring notation, the medical examiner’s report, any prior medical records, and a comprehensive narrative that summarized all scores of information about each subject. The social worker, the clinical psychologist and the psychiatrist reached a consensus on the diagnosis. Subjects meeting DSM-IV criteria for BD or MDD, or controls without a psychiatric diagnosis were selected for this study. Subjects were excluded for any neuropathological or neurological disorders. The presence of psychotropic medications and substances of abuse in blood and urine was determined by the medical examiner’s office.

Statistical power statement

Based on unadjusted analyses, two-sided tests with α = 0.05, the minimum detectable effect size in terms of absolute log2 fold change between groups of BD vs MDD/controls was 0.77 while that of MDD vs controls was 0.90. When using α = 0.01, the minimum detectable absolute log2 fold change of BD vs MDD/controls was 0.96, while that of MDD vs controls was 1.15.

Tissue

Postmortem brains were collected at autopsy at the Cuyahoga County Medical Examiner’s Office, Cleveland, OH. Cause of death was determined by the medical examiner. Tissues were dissected and rapidly frozen in 2-methylbutane on dry ice without fixation, and were kept in dry ice during transportation before permanent storage at −80 °C. Gray matter tissues of DLPFC (BA9) and anterior TP (BA38) from 20 BD, ten MDD, and ten control subjects were procured. MDD and control subjects were matched to BD subjects according to the group’s age and sex distributions (i.e., frequency matching), while control subjects were age-and-sex matched to MDD subjects individually. Subject and tissue information are listed in Table 1 (Supplementary Table 1 for details).

Table 1

Subject demographics.

	Controls	MDD	BD	p-value
n	10	10	20
Male (n; %)	8 (80%)	8 (80%)	17 (85%)	0.917
Age (years; mean ± SD)	48.4 ± 13.4	48.7 ± 20.0	47.5 ± 14.1	0.966
Hx Psychosis (n; %)	0	1 (10%)	11 (55%)	0.002
Substance abuse/dependence (n; %)
Nicotine	4 (40%)	5 (50%)	15 (75%)	0.144
Alcohol	0	0	9 (45%)	0.002
Causes of death (n; %)
Disease/accident	10 (100%)	1 (10%)	7 (35%)	< 0.001
Suicide	0	9 (90%)	13 (65%)
Postmortem interval (hr; mean ± SD)	24.0 ± 5.0	24.4 ± 8.3	22.8 ± 7.9	0.821
pH (mean ± SD)	6.5 ± 0.3	6.6 ± 0.2	6.5 ± 0.4	0.783

Note: Controls and MDD subjects were pair-matched by age and sex, while BD subjects were matched with the other two groups by frequency.

Genome-wide DNA methylation profiling

Genomic DNA was extracted from ~40 mg tissue by Gentra Puregene DNA Extraction Kit (QIAGEN Inc., Germantown, MD, USA). Bisulfite conversion and methylomic microarray were performed by the Mayo Clinic Medical Genome Facility Genotyping Core in a single batch. 25 μL of 50 ng/μL DNA underwent bisulfite conversion using EZ DNA Methylation Kit (Zymo Research, Irvine, CA, USA). Bisulfite-treated DNA was then hybridized with the 850 K Infinium MethylationEPIC BeadChip (Illumina Inc., San Diego, CA, USA), which covered over 850 000 methylation sites located in CpG islands, genes, transcription binding sites, open chromatin regions, and enhancers at single-nucleotide resolution. Unmethylated and methylated CpGenome controls were included and duplicates of a pooled DNA sample were included in each array to assess inter-array consistency.

Data processing

Data quality control and statistical analyses were performed in R software v.3.4.1 with Bioconductor package minfi (Aryee et al., 2014). No sample was excluded due to low bisulfite conversion ratios or call rates (detection p < 0.05). Probes that failed in one or more samples, or had a single nucleotide polymorphism at the CpG site, or mapped to non-specific genomic locations were removed. The probes located on X and Y chromosomes were excluded from data quality control procedures but included in subsequent statistical analyses. The initial methylation data contained 866 091 probes. After quality control, 704 741 probes were available for analysis. Data were normalized by functional normalization (preprocessFunnorm command in minfi) which removes between-array variation by regressing out variability in the control probes on each array (Fortin et al., 2014). Probes on X and Y chromosomes were normalized according to the sex of the sample.

Multivariate analyses

Principal component analysis (PCA) was performed on the top 2000 most variable CpG probes considering all samples (CpG probes with largest standard deviations in M-values). Hierarchical clustering was performed on the top 5000 most variable CpG probes using Euclidean distances as a measure of dissimilarity and k = 3 clusters for unsupervised clustering within each brain region. Association between the clusters and clinical factors were tested by chi-square and Fisher’s exact tests.

Identification of differentially methylated positions

Differentially methylated positions (DMPs; methylation status at individual site) were identified by the R Bioconductor package limma (Ritchie et al., 2015). For between-group comparisons in each brain area, pairwise differential methylation analyses were performed comparing BD with MDD, BD with controls, and MDD with controls using standard limma workflow with unpaired or paired contrasts (depended on group matching designs) with age and sex as covariates. For between-brain-region comparisons in each subject group, paired differential methylation analyses were performed between BA9 and BA38 using standard limma workflow with paired contrasts. M-values were used for statistical analysis and β-values were reported for interpretability. Statistically significant DMP was considered at adjusted p < 0.05. DMPs were considered as shared between comparisons if the same probe identification number (i.e., genomic location) was found to be significantly different in multiple comparisons. We further performed interaction tests to determine whether the magnitudes of between-region difference were different between groups: the methylation difference between BA9 and BA38 was computed at each CpG probe and then compared between groups with age and sex adjustment. Moderated paired t-test was used for the control-MDD comparison.

Identification of differentially methylated regions

Differentially methylated regions (DMRs; methylation status across a genomic region) were identified by the R Bioconductor package DMRcate (Peters et al., 2015). A DMR was called when a region contained at least two CpG sites with FDR-adjusted p values at default detection level within a lambda of 1000. Statistically significant DMR was considered at Stouffer p < 0.05. The criteria for shared DMRs between comparisons are described in Supplementary Methods.

Functional enrichment analysis

Ingenuity Pathway Analysis (IPA; QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis) was performed on genes annotated to the lists of significant DMPs that were shared by all groups as well as shared by and distinct to BD and MDD. IPA employs Fisher’s exact test to detect significant enrichment. The Illumina Human MethylationEPIC platform was used as the reference population gene set. The content was limited to human brain tissues (“brain”, “cerebral cortex” and “gray matter”). Pathways were considered statistically significant at Benjamini-Hochberg (B–H) corrected p < 0.05.

Gene expression quantification by RT-qPCR

Details of RNA extraction and real-time quantitative polymerase chain reaction (RT-qPCR) are described in Supplementary Methods. Relative gene expression levels were calculated using 2−ΔΔCt method with GAPDH as the housekeeping gene and control BA9 as the reference group, and were compared by repeated measures two-way ANOVA adjusted for age and sex.

Results

Subject and tissue characteristics

The summary of subject and tissue characteristics is presented in Table 1. Sex, age, postmortem interval, tissue pH, and nicotine dependence history were not significantly different among groups. Most of the subjects in the BD and MDD groups died of suicide (65% and 90% respectively) and some had psychosis history (55% and 10% respectively) in contrast to none in the control group (Fisher’s exact tests p < 0.01). The BD group also had more subjects with alcohol abuse/ dependence history (45%; Fisher’s exact test p = 0.002).

Overview of DNA methylomic differences between subject groups and brain regions

The PCA plot showed that the DNA methylomic profiles of samples belonging to the same brain region formed two clusters, but no distinct clustering or ordering was observed between samples belonging to different subject groups (Fig. 1A). K-means clustering analysis showed no associations between resulted clusters and clinical variables including subject groups, smoking history, alcohol abuse history, and postmortem interval, apart from the cause of death (suicide vs not suicide) in BA38 (p = 0.004; Supplementary Fig. 1). Collectively, these unsupervised clustering analysis results suggested that the DNA methylomic profiles of the two brain regions were distinct but not for the profiles of each subject group within a brain region.

Fig. 1.

Epigenome-wide brain regional DNA methylation differences among bipolar disorder (BD), major depression disorder (MDD) and controls. (A) Principal component analysis plot revealed two discrete clusters which corresponded to BA9 (filled circles) and BA38 (open circles). (B) Manhattan plots of methylation differences between BA9 and BA38 at individual probe level in each subject group. Red dotted line denotes p = 1E-6. (C and D) Between-brain-region differentially methylated positions (DMPs; adjusted p < 0.05) and differentially methylated regions (DMRs; Stouffer p < 0.05) are highly overlapped among subject groups. Numbers in the partitions represent the number of DMPs or DMRs that are shared by corresponding groups in overlapped partitions or specific to a certain group in non-overlap partitions. (E) The distribution of genetic elements mapped to between-brain-region DMPs (adjusted p < 0.05) found in all groups was distinct from the distributions of the BD-MDD shared, BD-only, and MDD-only DMPs with fewer DMPs located in or proximal to genes (chi-square test p = 4.04E-42) but more DMPs located in or proximal to CpG islands (chi-square test p = 1.52E-68). *Genetic elements with marked differences in distribution between groups; ExonBnd: exon boundary; TSS: transcription start site; UTR: untranslated region.

Methylation differences in DLPFC (BA9) between subject groups

The Manhattan plots of the DMP comparisons between groups in BA9 are presented in Supplementary Fig. 2. No DMPs or DMRs reached the statistical significance thresholds in any pairwise comparison in BA9 (top between-group DMPs and DMRs are listed in Supplementary Tables 2 and 3). Pathway analyses performed on the list of genes annotated to the top 100 DMPs in any BA9 pairwise comparison detected no significantly enriched pathway (Fig. 2A).

Fig. 2.

Functional enrichment analysis results of the genes mapped to the top 100 differentially methylated probes (DMPs) in the between-group comparisons in (A) BA9 and (B) BA38, and (C) of the genes associated with between-brain-region DMPs shared by all groups, shared by bipolar disorder (BD) and major depressive disorder (MDD) groups, and were distinct to BD and MDD groups. Shared DMPs were defined as significant between-brain-region DMPs (adjusted p < 0.05) having the same probe identification number across various subject groups. All pathways shown are involved in neurotransmitters and other nervous system signaling in Ingenuity Pathway Analysis. Shades of colors represent unadjusted Fisher’s exact test p-values. Dots inside box indicate pathways with Benjamini-Hochberg corrected p < 0.05.

Methylation differences in TP (BA38) between subject groups

The Manhattan plots of the DMP comparisons between groups in BA38 are presented in Supplementary Fig. 3. No DMPs or DMRs reached the statistical significance thresholds in any pairwise comparison in BA38 (top between-group DMPs and DMRs are listed in Supplementary Tables 2 and 3). Pathway analyses performed on the list of genes annotated to the top 100 DMPs in any BA38 pairwise comparison detected no significantly enriched pathway (Fig. 2B).

DNA methylation differences between BA9 and BA38 are highly shared among subject groups

The Manhattan plots of the DMP comparison between brain regions for each subject group are shown in Fig. 1B. Between BA9 and BA38, 1601 DMPs in the control group, 11 954 DMPs in the MDD group and 39 039 DMPs in the BD group reached statistical significance. Notably, there are a few peaks apparently aligned to the same positions in all subject groups, such as in chromosomes 3, 5 and 15. Many significant DMPs were found in multiple subject groups (Fig. 1C). Full lists of these DMPs included in each distinct or overlap group are presented in Supplementary Table 4. The large number of DMPs detected between BA9 and BA38 and their distribution in the genome across subject groups suggested that some DMPs were densely populated in certain genomic regions. Between-brain-region DMRs also showed a similar degree of overlap among groups (Fig. 1D; full list in Supplementary Table 5). Since some of the shared DMRs did not have identical genomic boundaries and the methylomic differences in the BD group were our primary interest, we focused on the between-brain-region DMPs distinct to BD or observed in both BD and MDD.

Between-brain-region DNA methylation differences distinct to BD subjects only

The significant methylomic differences distinct to the BD group may be associated with clinical manifestation in mood states, mania, and higher frequency of substance abuse/dependence in these subjects. The top 20 between-brain-region DMPs distinct to BD were listed in Table 2 (full list in Supplementary Table 4). The genes associated with these DMPs were significantly enriched in 15 nervous system-associated pathways (Fig. 2C; details in Supplementary Table 6), including the top-ranked “axonal guidance signaling” pathway (223/319 hits). The between-brain-region methylation differences of some of these DMPs also showed trended differences between BD and controls and between BD and MDD but not between MDD and controls (interaction tests at nominal p < 0.05). We selected a few genes (ODC1, TF, CALY) which were annotated to these DMPs since we expected their BD-specific between-brain-region methylation levels might mirror their gene expression levels. The gene expression levels of ODC1 and CALY showed significant region × group interaction, indicating that the between-brain-region differences were dependent on subject groups (Fig. 3E and F). For TF, No significant difference was detected between brain regions or between groups (Supplementary Fig. 4A).

Table 2

Top differentially methylated positions significantly different (adjusted p < 0.05) between BA9 and BA38 found in the BD groups only. For each subject group, the methylation differences between regions (BA38 vs BA9) are presented in the Log2 fold change (FC) columns, and beside them are the corresponding adjusted p values. Interaction tests were performed to test whether the between-brain-region methylation differences were different between groups; their unadjusted p-values are listed on the right-most columns since no adjusted p < 0.05. The top 30 DMPs are listed in descending adjusted p values of the BD group. DMPs in bottom rows were selected for gene expression analysis based on interaction tests nominal p < 0.05 in controls vs BD and BD vs MDD, i.e., the BD between-brain-region difference is most likely to be disparate from controls and MDD).

Probe ID	Chr	Position	Gene	Elements	Log2 FCBA38/BA9			BA38vsBA9 within group Adjusted p-value			BA38-BA9 between groups p-value[a]
					Ctl	MDD	BD	Ctl	MDD	BD	CtlvsBD	BDvsMDD	CtlvsMDD
cg20803618	2	100210271	AFF3	Body	−0.671	−0.530	−0.929	1.79E-01	6.66E-02	1.12E-08	0.105	0.013	0.970
cg11864566	X	99667214			−0.400	−0.510	−0.600	4.19E-01	7.10E-02	3.19E-07	0.117	0.562	0.313
cg00610310	16	1128740	SSTR5	TSS200	0.619	0.427	0.551	7.02E-02	5.25E-02	4.09E-07	0.516	0.284	0.327
cg09043518	17	80290106	SECTM1	5′UTR	0.934	0.600	0.959	1.50E-01	8.41E-02	4.23E-07	0.900	0.053	0.617
cg18989491	10	6185502	PFKFB3	TSS1500	−0.489	−0.448	−0.739	2.92E-01	9.40E-02	5.34E-07	0.075	0.048	0.749
cg09560062	18	41143802			−0.627	−0.469	−0.742	1.63E-01	1.17E-01	7.30E-07	0.451	0.094	0.441
cg21242356	17	46667683	LOC404266	TSS200	−0.424	−0.332	−0.485	2.50E-01	6.96E-02	8.03E-07	0.580	0.116	0.827
cg22915785	8	144120111	C8orf31	TSS1500	0.444	0.511	0.734	5.96E-01	1.24E-01	9.07E-07	0.258	0.178	0.136
cg11714334	22	39638092	PDGFB	Body	−1.109	−0.742	−1.352	1.76E-01	2.17E-01	9.70E-07	0.376	0.046	0.408
cg16233888	8	61842553			0.708	0.561	1.011	1.62E-01	5.69E-02	1.14E-06	0.128	0.026	0.814
cg19318393	1	223936508	CAPN2	Body	0.465	0.547	0.791	4.89E-01	6.15E-02	1.16E-06	0.133	0.087	0.786
cg25344672	9	135285207	C9orf171	TSS1500	−0.562	−0.359	−0.492	1.17E-01	8.35E-02	1.62E-06	0.531	0.225	0.338
cg17696563	2	242756362	NEU4	Body	1.398	1.126	1.712	1.71E-01	1.41E-01	1.68E-06	0.411	0.168	0.641
cg10303842	5	22853730	CDH12	5′UTR	−0.344	−0.386	−0.606	3.24E-01	6.34E-02	1.87E-06	0.047	0.077	0.575
cg05760107	16	81450558			−0.352	−0.277	−0.536	2.92E-01	1.16E-01	1.88E-06	0.094	0.011	0.941
cg27613737	21	47038433			0.925	0.520	0.947	6.49E-02	7.89E-02	2.03E-06	0.913	0.021	0.124
cg00303183	2	162283559			−0.438	−0.308	−0.509	9.95E-02	1.38E-01	2.19E-06	0.480	0.073	0.823
cg02326285	14	79111719	NRXN3	5′UTR	0.569	0.396	0.700	1.51E-01	1.43E-01	2.41E-06	0.392	0.047	0.272
cg21300267	1	26547059			0.549	0.359	0.568	9.37E-02	7.27E-02	2.42E-06	0.900	0.065	0.134
cg23097006	20	25063817	VSX1	TSS1500	−0.875	−0.624	−1.088	1.08E-01	6.08E-02	2.59E-06	0.371	0.044	0.892
cg21811204	5	132158625	SHROOM1	Body	−0.564	−0.475	−0.601	9.32E-02	6.44E-02	2.62E-06	0.758	0.375	0.948
cg19696441	1	20512463	UBXN10	TSS200	−0.668	−0.697	−1.014	3.12E-01	9.36E-02	2.91E-06	0.142	0.198	0.179
cg12509733	10	131647800	EBF3	Body	0.457	0.331	0.493	2.61E-01	2.19E-01	3.14E-06	0.767	0.175	0.715
cg15042460	11	63994295	TRPT1	TSS1500	0.246	0.497	0.894	5.97E-01	1.68E-01	3.15E-06	0.002	0.031	0.187
cg13812927	9	89517863			−0.427	−0.557	−0.833	3.38E-01	1.37E-01	3.24E-06	0.021	0.188	0.799
cg07755196	10	102995440	FLJ41350	Body	−0.368	−0.276	−0.431	1.39E-01	6.74E-02	3.38E-06	0.468	0.062	0.350
cg24931954	11	3164851	OSBPL5	5′UTR	0.205	0.681	0.846	6.25E-01	8.89E-02	3.76E-06	0.002	0.454	0.011
cg20282967	17	61916136	SMARCD2	Body	0.710	0.537	0.785	3.91E-01	1.40E-01	3.88E-06	0.746	0.189	0.541
cg06809252	1	110612044	ALX3	Body	−0.411	−0.489	−0.848	3.98E-01	2.46E-01	3.93E-06	0.030	0.102	0.169
cg07952815	14	104073479			0.554	0.231	0.616	1.27E-01	3.01E-01	4.03E-06	0.631	0.004	0.036
Selected for Gene Expression Analysis
cg17231641	2	10587759	ODC1	5′UTR	−0.257	−0.346	−0.648	5.62E-01	1.13E-01	4.82E-06	0.017	0.015	0.392
cg24662360	3	133465136	TF	5′UTR	0.172	0.237	0.811	7.95E-01	7.78E-01	1.32E-04	0.016	0.016	0.950
cg26583481	10	135139100	CALY	3′UTR	0.406	0.383	0.868	4.59E-01	1.08E-01	1.46E-04	0.043	0.005	0.920

All adjusted p-values in these comparisons were > 0.985; bold values: unadjusted p < 0.05.

Fig. 3.

Regional plots of between-brain-region methylation differences of (A) ODC1, (B) CALY, (C) GALNT2, and (D) GABRD. Each point represents the Log2 methylation fold change BA38 vs BA9 at a CpG site (gray for controls, blue for MDD, and red for BD). Dotted lines indicated the differentially methylated positions (DMPs) on which gene selection was based. The detailed plots of between-brain-region methylation differences in relation to CpG islands and DNaseI hypersensitivity regions are shown in Supplementary Figs. 5–8. Relative gene expression levels (estimated mean ± standard deviation with age and sex adjustments) of these genes are shown in corresponding bar charts. (E) ODC1 gene expression level of BD was significantly higher in BA38 than BA9 (two-way repeated measure ANOVA F2,32 = 3.351, pgroup×region = 0.048; posthoc Bonferroni p = 0.032), and such difference opposed the trends shown in the controls (BA38 < BA9) and MDD (BA38 ≈ BA9).(F) CALY gene expression level of controls was significant lower in BA38 compared to BA9 (two-way repeated measure ANOVA F2,32 = 8.513, pgroup×region = 0.001; posthoc Bonferroni p = 0.001), in contrast to the between-brain-region differences of BD (BA38 > BA9) and MDD (BA38 < BA9 with smaller difference). In BA38, BD showed a marginally higher CALY gene expression level than controls (posthoc Bonferroni p = 0.080). (G) GALNT2 gene expression level of BD was significantly higher in BA38 than BA9 (two-way repeated measure ANOVA F2,32 = 8.815, pgroup×region = 0.001; posthoc Bonferroni p = 0.048). In contrast, GALNT2 gene expression level of controls was significantly lower in BA38 than BA9 (posthoc Bonferroni p = 0.024), while that of MDD was not significantly different between regions. GALNT2 gene expression level was significantly higher in BD compared to controls in BA38 (posthoc Bonferroni p = 0.040) while the reverse was observed in BA9 with marginal significance (posthoc Bonferroni p = 0.058). (H) GABRD gene expression levels of controls and MDD were significantly lower in BA38 than BA9 (two-way repeated measure ANOVA F2,32 = 5.138, pgroup×region = 0.001; posthoc Bonferroni pCtl < 0.001 and pMDD = 0.045) but no significant difference was observed between brain regions in BD. The BD group showed a relatively lower GABRD gene expression level in BA9, but a higher level in BA38, compared to controls. *Bonferroni posthoc test p < 0.05; #Bonferroni posthoc test p < 0.10.

Between-brain-region DNA methylation differences shared between BD and MDD subjects

The significant methylomic differences shared by both BD and MDD groups may be associated with the characteristics shared by these illnesses, such as depression and suicide. The top 20 between-brain-region DMPs shared by BD and MDD were listed in Table 3 (full list in Supplementary Table 4). The genes annotated to these DMPs were significantly enriched in seven nervous system-associated pathways (Fig. 2C; details in Supplementary Table 6), including the top-ranked “axon guidance signaling” pathway (89/319 hits). The between-brain-region methylation differences of some of these DMPs also showed trended differences between BD and controls and between MDD and controls but not between BD and MDD (interaction tests at nominal p < 0.05). We again selected a few genes (GALNT2, RUNX1T1, GABRD) annotated to these DMPs to compare whether between-group differences in their between-brain-region methylation levels might mirror their between-brain-region gene expression levels. The gene expression levels of GALNT2 and GABRD showed significant region × group interaction, indicating that the between-brain-region differences were dependent on subject groups (Fig. 3G and H). In RUNX1T1, no significant difference was detected between brain regions or between groups (Supplementary Fig. 4B).

Table 3

Top differentially methylated positions significantly different (adjusted p < 0.05) between BA9 and BA38 found in both BD and MDD groups. For each subject group, the methylation differences between regions (BA38 vs BA9) are presented in the Log fold change (FC) columns, and beside them are the corresponding adjusted p values. Interaction tests were performed to test whether the between-brain-region methylation differences were different between groups; their unadjusted p-values are listed on the right-most columns since no adjusted p < 0.05. The top 30 DMPs are listed in descending adjusted p values of the BD group. DMPs in bottom rows were selected for gene expression analysis based on interaction tests nominal p < 0.05 in controls vs BD and controls vs MDD, i.e., the between-brain-region differences in BD and MDD are most likely to be disparate from controls).

Probe ID	Chr	Position	Gene	Elements	Log2 FCBA38/BA9			BA38vsBA9 within group Adjusted p-value			BA38-BA9 between groups p-value[a]
					Ctl	MDD	BD	Ctl	MDD	BD	CtlvsBD	BDvsMDD	CtlvsMDD
cg26536894	22	43659884	SCUBE1	Body	−1.008	−1.112	−1.114	2.16E-01	7.49E-05	3.98E-11	0.546	0.960	0.899
cg22018084	2	69038738	ARHGAP25	Body	−1.284	−1.383	−1.296	5.18E-02	6.09E-05	4.35E-10	0.932	0.621	0.830
cg10828175	6	159528564			−0.845	−1.144	−1.033	8.59E-02	2.97E-07	6.44E-10	0.214	0.295	0.134
cg06717565	12	114836487	TBX5	Body	−1.098	−1.062	−1.130	1.48E-01	5.76E-05	6.77E-10	0.877	0.745	0.315
cg20099830	12	114842031	TBX5	5′UTR	−1.247	−1.285	−1.332	6.17E-02	1.20E-05	4.24E-09	0.709	0.856	0.684
cg02554050	17	72462691	CD300A	5′UTR	−0.539	−0.588	−0.616	1.16E-01	1.42E-03	5.16E-09	0.454	0.764	0.854
cg20626840	2	38201011	FAM82A1	Body	−0.451	−0.723	−0.687	3.81E-01	1.07E-03	5.26E-09	0.084	0.638	0.272
cg26219540	5	92937739			−1.529	−1.492	−1.691	5.79E-02	5.64E-07	6.62E-09	0.569	0.273	0.941
cg15114001	2	136899351			−0.354	−0.537	−0.679	3.64E-01	1.24E-02	8.49E-09	0.005	0.162	0.291
cg13386858	8	127236702			1.435	1.461	1.356	9.50E-02	1.08E-04	1.03E-08	0.770	0.549	0.829
cg16709475	12	113904823	LHX5	Body	1.095	1.088	1.343	1.48E-01	1.01E-02	1.11E-08	0.272	0.198	0.934
cg03589296	2	66666292	MEIS1	Body	−1.021	−0.840	−0.849	7.44E-02	3.40E-03	1.23E-08	0.320	1.000	0.573
cg24984195	16	1361783	UBE2I		0.865	0.944	1.036	1.08E-01	7.19E-05	1.43E-08	0.348	0.411	0.653
cg06048910	11	2925473	SLC22A18	Body	0.994	0.942	1.379	1.13E-01	2.93E-02	1.50E-08	0.109	0.071	0.729
cg24298860	11	92600238	FAT3	Body	0.582	0.965	0.908	7.67E-02	6.26E-03	1.52E-08	0.019	0.738	0.090
cg12449515	1	230333800	GALNT2	Body	0.603	1.099	1.121	1.48E-01	2.54E-04	1.56E-08	0.004	0.973	0.016
cg02362103	1	119532773	TBX15	TSS1500	−0.804	−1.023	−1.095	1.27E-01	1.96E-05	1.63E-08	0.128	0.544	0.014
cg04658543	12	113904848	LHX5	Body	1.032	1.104	1.236	6.00E-02	2.88E-04	1.72E-08	0.345	0.410	0.516
cg25112291	18	56936699	RAX	Body	−0.705	−0.791	−0.720	8.13E-02	2.95E-04	1.86E-08	0.899	0.364	0.177
cg18150909	13	50706830	DLEU1	Body	−0.585	−0.692	−0.745	1.92E-01	2.50E-03	2.02E-08	0.256	0.665	0.185
cg14489049	9	135317080	C9orf171	Body	1.028	1.077	1.012	1.41E-01	2.75E-04	2.07E-08	0.958	0.694	0.795
cg13628942	12	114845747	TBX5	5′UTR	−0.975	−1.110	−1.269	8.82E-02	1.21E-06	2.11E-08	0.190	0.328	0.194
cg15117372	6	35992219	SLC26A8	5′UTR	−0.349	−0.638	−0.621	3.37E-01	1.38E-03	2.47E-08	0.014	0.737	0.009
cg21006325	21	44480624	CBS	Body	1.254	1.151	1.324	9.60E-02	1.11E-02	2.65E-08	0.815	0.471	0.914
cg04701618	12	131443231	GPR133	Body	1.543	1.529	1.367	9.95E-02	4.26E-04	2.66E-08	0.542	0.415	0.749
cg17869960	6	12718530	PHACTR1	5′UTR	−0.834	−0.926	−0.992	2.05E-01	2.52E-04	3.33E-08	0.385	0.661	0.998
cg22821289	11	115581830			0.848	1.435	1.412	3.74E-01	1.35E-04	3.34E-08	0.066	0.840	0.069
cg08732352	7	35298391			−1.488	−1.725	−1.712	5.34E-02	1.47E-03	3.89E-08	0.432	0.957	0.127
cg03706358	15	93043867	C15orf32	3′UTR	0.722	0.870	1.044	2.30E-01	1.53E-03	4.20E-08	0.041	0.208	0.325
cg26108046	3	183535040	MAP6D1	3′UTR	1.080	1.115	1.186	6.70E-02	1.62E-03	4.34E-08	0.585	0.777	0.150
Selected for Gene Expression Analysis
cg12449515	1	230333800	GALNT2	Body	0.603	1.099	1.121	1.48E-01	2.54E-04	1.56E-08	0.004	0.973	0.016
cg26959358	8	93107604	RUNX1T1	TSS200	−0.296	−0.766	−0.668	4.41E-01	5.66E-03	2.21E-06	0.011	0.498	0.024
cg14522944	1	1961743	GABRD	3′UTR	0.044	0.501	0.598	8.79E-01	4.82E-02	1.54E-04	0.002	0.503	0.040

All adjusted p-values in these comparisons were > 0.985; bold values: unadjusted.

Effect of suicide on methylation differences in BD group

Since the cause of death was significantly associated with methylation clusters (see section 3.2), we checked the potential effect of suicide in the BD group in BA38 and between-brain-region differences. Control and MDD groups were not included since all controls died of illness in contrast to nearly all MDD subjects died of suicide in addition to inter-group variations. No significant DMPs or DMRs were detected between those died of suicide and of illnesses in BA9, BA38, or between-brain-region differences.

In an attempt to screen for a methylomic profile associated with suicide completion, we also compared the between-brain-region methylomic differences in BD + MDD cases that died of suicide to controls, and between-brain-region methylomic differences in BD + MDD cases that died of illness to controls. No significant DMPs and DMRs were detected in either of these comparisons, which might be partly attributed to unmatched critical group characteristics (e.g., age, sex, methods of suicide).

Discussion

In identifying DNA methylation differences associated with BD in BA9 and BA38, we discovered that the most striking between-group disparities lay in the between-brain-region methylation differences. Some of these disparities could potentially contribute to group-specific between-brain-region gene expression differences as demonstrated by ODC1, CALY, GALNT2, and GABRD. The lists of genes mapped to the DMPs distinct to BD and shared by MDD and BD were enriched primarily in pathways related to neurodevelopment and synaptic plasticity.

Structural abnormalities have been observed in both BA9 and BA38 of mood disorder subjects including reduced gray matter volume and cortical thickness (Hibar et al., 2017; Price and Drevets, 2010), supporting the association between these two brain regions and mood regulation. Despite these regions not being directly connected (Catani and Thiebaut de Schotten, 2008), they are connected via lateral orbitofrontal cortex (BA11 and BA47), frontal pole cortex (BA10), and parahippocampal gyrus (Catani et al., 2002; Catani and Thiebaut de Schotten, 2008; Price and Drevets, 2010; Yeterian et al., 2012). These indirect connections via prefrontal and limbic areas may allow communication between BA9 and BA38 for regulating emotion processing, verbal expression, reward response, impulse control, decision making, and socio-affective processing, dysfunctions of which are closely related to BD symptoms. Thus, comparing these regions may elucidate the underlying neural networks contributing to BD symptomology and neuropathology.

Among the top genes of the BD-only list are genes that participate in neurodevelopment and had demonstrated genetic associations with BD, e.g., PDGFB, CDH12, and NRXR3 (Kataoka et al., 2016; Kuo et al., 2014; Redies et al., 2012), while among the top genes of BD-MDD shared list are genes previously associated with BD, circadian rhythmicity, or MDD, e.g., LHX5, MEIS1, and GALNT2 (Davidson et al., 2006; Ferguson et al., 2018; Gao et al., 2015; Green et al., 2003). In selected genes of these lists (BD-only: ODC1 and CALY; BD-MDD shared: GALNT2 and GABRD), we detected group-specific between-brain-region gene expression differences of which BD showed an opposite gene expression balance between BA9 and BA38 compared to controls while MDD lay in between, suggesting the between-brain-region imbalance of methylation levels of these genes may contribute to between-brain-region imbalance of their gene expression levels that are associated with BD and MDD. These results support BD and MDD etiology and/or neuropathology may involve between-brain-region dysregulations of polyamine biosynthesis (ODC1 function), receptor endocytosis and recycling (CALY function), lipid metabolism and glucose homeostasis (GALNT2 function), and tonic inhibition via GABAA receptor (GABRD function). Nonetheless, it should be noted that not all between-brain-region methylation differences among groups detected would demonstrate differential gene expression (as in TF and RUNX1T1) and the between-brain-region gene expression differences of MDD might not be similar to BD even though the corresponding DMPs belonged to the BD-MDD shared list.

Hypermethylation in the promoter region is generally associated with gene down-regulation while hypermethylation in the gene body is generally associated with gene up-regulation (Suzuki and Bird, 2008; Yang et al., 2014); hypermethylation in 5′untranslated region (UTR) and 3′UTR was reported to be negatively and positively correlated with gene expression respectively (Mishra and Guda, 2017). Located in the 5′UTR of ODC1, cg17231641 had a negative log2 fold change (i.e., hypomethylation in BA38 compared to BA9) in all groups and was particularly low in BD, hence a higher gene expression level of ODC1 in BA38 compared to BA9 was expected. Such trend was only observed in BD, while MDD showed comparable between-brain-region gene expression levels and controls displayed the opposite trend of lower gene expression in BA38 than BA9. cg26583481 and cg14522944 were located in the 3′UTRs of CALY and GABRD respectively, and cg12449515 was located in an intron of GALNT2. These DMPs had positive log2 fold changes (i.e., hypermethylation in BA38 compared to BA9) in all groups with distinct differences either in BD only or in both BD and MDD, hence these genes were expected to show increased gene expression levels in BA38 compared to BA9. However, such increase in gene expression was only observed in GALNT2 in BD, while controls and MDD demonstrated the opposite trend in all three genes. These results imply that these between-brain-region DMPs might regulate gene transcription accompanying other mechanisms, such as miRNA action and interaction with other upstream regulators and transcription factors. A wider coverage of CpG sites may reveal other sites with stronger correlations between methylation and gene expression levels.

Most of the canonical pathways significantly enriched by the genes associated with DMPs shared by BD and MDD and those distinct to BD are critical in neurodevelopment, including axon guidance, CREB signaling, and synaptic long-term depression. Axon guidance signaling is the top enriched pathway for both BD-MDD and BD-only gene lists but not for the list shared by all groups. Many of the molecules involved in axon guidance continue to be expressed in the brain throughout life to maintain the integrity and plasticity of neuronal connections (Lin et al., 2009; Shen and Cowan, 2010). The enrichment of axon guidance signaling supports the notion that the etiologies of BD and MDD may share some neurodevelopmental deficits in BA9 and/or BA38 which may persist into adulthood (Parellada et al., 2017). A significant enrichment in methylation of the GABAergic receptor signaling pathway was observed only in the BD-MDD shared list, which is important given the association between GABAergic systems and depression and between cortical inhibition and suicidality (Lewis et al., 2018). Meanwhile, a significant enrichment in glutamate receptor signaling was observed solely in the BD-only list. These differences in how GABAergic and glutamatergic pathways are affected could be a consequence of some commonality between BD and MDD (i.e., depression), of the unique features of BD (i.e., mania), or of the higher prevalence of certain clinical characteristics in our BD samples (i.e., psychosis). Other pathways enriched in the BD group that are involved in synaptogenesis and synaptic plasticity, such as netrin signaling and synaptic long-term potentiation/ depression, emphasize the deficits of these mechanisms in BD pathology in BA9 and/or BA38 or in networks which they participate.

A few limitations need to be noted. Firstly, our results lacked cell type specificity. Future methylomic studies could apply technologies, such as laser capture microdissection and fluorescence-activated cell sorting, to separate the cell type of interest with microscopic regional specificity (McCullumsmith and Meador-Woodruff, 2011; Ruzicka et al., 2018). Secondly, our method could not differentiate DNA methylation from DNA hydroxymethylation, which upregulates gene expression in contrast to DNA methylation (Feng et al., 2017; Song et al., 2011), thus reducing the predictability of gene expression levels inferred by our results. Thirdly, our sample size was limited and unbalanced among groups. Our current sample size only allowed the detection of large variations. With more subjects, the BD group had more statistical power to detect smaller significant differences than the MDD and control groups. Future studies should include larger brain sample cohorts and balanced sample sizes between groups. Fourthly, some clinical factors known to modulate DNA methylation, such as last recorded mood state, nature of death, medications, and substance use history, could not be perfectly matched between subject groups (Dell’Osso et al., 2014; Guidotti and Grayson, 2014; Huzayyin et al., 2014; Jia et al., 2017; Lohoff et al., 2017; McCullumsmith and Meador-Woodruff, 2011), thereby preventing statistical adjustment. Our findings might be partly attributed to these factors which were particularly complicated in BD and MDD cases. However, it should be noted that our unsupervised hierarchical clustering analysis did not detect significant associations between clusters and the history of smoking and alcohol abuse (Supplementary Fig. 1). Lastly, despite the MethylationEPIC BeadChip covering over 850 000 CpG sites located in coding regions and enhancers, this only represents ~3% of the estimated total 28.3 million CpG dinucleotides in the human genome (Luo et al., 2014). Such coverage is limited compared to other genome-wide methods (e.g., whole genome bisulfite sequencing) and potentially precludes the discovery of many DMPs and DMRs.

In summary, we identified DNA methylomic differences between DLPFC and TP that are associated with BD and MDD and associated with BD only but absent in controls. These methylomic differences were related to genes and pathways involved in axon guidance, glutamatergic and GABAergic neurotransmissions, synaptic plasticity and various neurodevelopmental mechanisms. These findings imply that a methylation imbalance between DLPFC and TP may contribute to BD symptomology. If so, methylomics could prove to be important targets for further investigations of genetic and environmental vulnerabilities associated with BD.