Characterizing the dynamic and functional DNA methylation landscape in the developing human cortex

ABSTRACT

DNA methylation (DNAm) is a key epigenetic regulator of gene expression across development. The developing prenatal brain is a highly dynamic tissue, but our understanding of key drivers of epigenetic variability across development is limited. We, therefore, assessed genomic methylation at over 39 million sites in the prenatal cortex using whole-genome bisulfite sequencing and found loci and regions in which methylation levels are dynamic across development. We saw that DNAm at these loci was associated with nearby gene expression and enriched for enhancer chromatin states in prenatal brain tissue. Additionally, these loci were enriched for genes associated with neuropsychiatric disorders and genes involved with neurogenesis. We also found autosomal differences in DNAm between the sexes during prenatal development, though these have less clear functional consequences. We lastly confirmed that the dynamic methylation at this critical period is specifically CpG methylation, with generally low levels of CpH methylation. Our findings provide detailed insight into prenatal brain development as well as clues to the pathogenesis of psychiatric traits seen later in life.

KEYWORDS:

• DNA methylation
• brain development
• prenatal development
• whole-genome bisulfite sequencing

Introduction

DNA methylation (DNAm) plays an important role in the epigenetic regulation of gene expression, tissue differentiation, and development from prenatal life through adolescence and likely beyond. It has been shown that in the human brain, DNAm is particularly plastic during the first 5 years of postnatal life, both at CpG and CpH (H = A, T, or C) sites [1]. Methylation levels at key sites change over time and these changes lead to adjustments in gene expression and splicing. These key regions have also been linked to neurodevelopmental disorders such as schizophrenia, in which early dysregulation plays a vital role [2,3]. DNAm is an attractive epigenetic mechanism to study in post-mortem human brain tissue, because it represents an interaction between genetic and environmental effects. External factors such as changes in diet [4], exposure to cigarette smoking [5], and exposure to arsenic [6] have been associated with both global and site-specific changes to DNAm levels.

In order to better understand the causes and consequences of deviant DNAm patterns in psychiatric disease development, we must first understand the normal landscape. Illuminating typical methylation changes in prenatal development will both provide insight into gene expression and molecular pathways active in the postnatal developing brain, and provide a baseline for identification of aberrant DNAm in postnatal disease states, as the pathological changes that lead to symptoms of psychiatric disease may precede the onset of illness by several decades [7].

The dorsolateral prefrontal cortex (DLPFC) is a dynamic region of the brain throughout development, essential for motor planning, conceptual organization, and working memory, and is often functionally dysregulated in patients with schizophrenia [8]. Previous studies using microarrays to quantify DNAm have revealed that there are many DNAm changes in the DLPFC around the time of birth [9], but the dynamics of the prenatal brain are far less characterized. Previous work assessing prenatal brain DNAm has largely used microarrays, and sampled from the whole brain rather than a specific region [10]. Furthermore, unlike the adult DLPFC which is a mixture of different cell types [11,12], the majority of cells in the developing brain are neuronal [13]. Research has further previously shown that neuronal DNA methylation is especially dynamic in the earliest stages of postnatal life [9,14]. These data together suggest that profiling the DLPFC in prenatal samples could identify important patterns of neuronal development that could relate to risk for neuropsychiatric disorders [1].

In the present report, we use whole-genome bisulfite sequencing (WGBS) to capture an unbiased map of the DNAm landscape and to characterize both CpG and CpH methylation during prenatal brain development. Dynamic regions of DNAm even in the early developmental context are likely to be associated with psychiatric risk-associated genes, and connection to gene expression data can validate the importance of these genomic regions. These important regions will help point to pathways and mechanisms of normal brain development as well as psychiatric and neurodevelopmental conditions.

Results

Whole-genome bisulfite sequencing in the prenatal human cortex

We performed whole-genome bisulfite sequencing (WGBS) to better characterize the shifting DNAm landscape in the developing human dorsolateral prefrontal cortex (DLPFC) in 20 prenatal samples during the second trimester in utero (Table S1). After data processing and quality control (see Methods), we analysed 28,612,402 CpG and 13,166,520 CpH (H = A, T, or C) sites across the epigenome. We first focused on CpG DNAm and performed a series of quality control analyses. We quantified read coverage at every CpG, which after processing and alignment resulted in an average of 17.7 reads per CpG (Figure 1(a)). Most of these sites were methylated (>80% DNAm) and a minority were unmethylated (<20% DNAm) (Table S2). We then performed principal component analysis (PCA) on the DNAm levels across these high-coverage CpGs and found that the top principal components were most associated with sex and estimated embryonic stem cell (ESC) fractions (which likely represent radial glia) via statistical deconvolution (Figure 1(b); Figure S1; see Methods), in line with our previous work [9]. We lastly compared the DNAm levels from WGBS to levels from Illumina 450k microarray measured on 15 of the 20 samples for the subset of CpGs in common to both technologies (N = 456,305), and found that overall, CpG methylation levels were highly concordant regardless of CpG read coverage (Figure 1(c, d)). These analyses together suggested that our DNAm data were of high quality for subsequent methylation analyses.

Figure 1. WGBS Data. (a) Coverage levels across all CpGs. The mean coverage level was 17.7 reads per CpG. (b) PCA on all CpG DNAm levels shows that the top principal component – accounting for 10.3% of variance, is sex. However we can also see that there is another factor separating the two clusters of females, shown here to be embryonic stem (ES) cell fraction. (c) Root-mean-square deviation (RMSD) of WGBS measurements from 450k measurements for each sample. Mean RMSD is shown. (d) WGBS measured DNAm for each site on 450k, organized by 450k methylation measurement, of sample LIBD399, shown in the red circle in (C)

DNAm associations and their characterization

To understand the changes in the prenatal epigenetic landscape, we performed linear modelling across all CpG sites. We found that DNAm changes were abundant even during this relatively restricted period in prenatal development, with 36,546 CpG sites associating with age across 14–20 post-conception weeks (PCW, at FDR < 0.05, Table S3). On average, each week of development was associated with a 1.8% change in DNAm (interquartile range (IQR): 1.18–2.33%, Figure 2(a)) among these significant age-associated CpGs, with some significant sites changing as much as 8.7% per week. These age-associated CpGs were not evenly distributed across the genome, and there was a noted depletion of signal on the sex chromosomes (0 on chrY and 339 on chrX, see Table S4). There was also some non-random distribution across the autosomes, with some chromosomes having up to 2.5x higher frequency of age-associated CpGs than others (Table S4). With regard to genomic features, 2,743 of these CpGs lie in previously annotated CpG islands, with 5,933 in shores, defined as 2 kb from island ends. Unlike microarray platforms, which are biased towards gene promoters, the vast majority (73%) of our age-associated CpGs lie within genes, consistent with previous findings on age-associated methylation using WGBS [1], and mostly in introns. Additionally, 2,932 age-associated CpGs lie within 1 kb of a transcription start site and thus potentially in promoters. A small percentage of these age-associated CpGs (6.5%) were annotated to intergenic sequences. Sites were fairly evenly split in the direction of methylation change by age, with 46% increasing in methylation as the cortex develops. These sites were fairly robust to confounding by ethnicity, as sensitivity analyses further adjusting for self-reported race showed little impact on age-associated CpG sites (Table S5). Additionally, less than half of these age-associated CpG sites are significantly associated with the estimated ESC proportion in the sample, suggesting that many of these CpG sites have true prenatal age changes in methylation rather than reflecting effects of maturing cell type composition.

Figure 2. Methylation associations across ages and between sexes. (a) Percent change in methylation per post-conception week (PCW) across all significantly age-associated CpGs (FDR < 0.05). A positive % change indicates hypermethylation with increasing age. (b) An example of an age-associated CpG on chromosome 10. (c) Percent change in methylation between sexes at sex-associated CpGs, showing that the X chromosome tends to have greater percentage differences than autosomal CpGs. A positive % change indicates hypermethylation in males. (d) An example of a sex-associated CpG on chromosome 2

We further explored whether these sites could be organized into differentially methylated regions (DMRs), which have been shown to be more functionally relevant than individual CpGs. We, therefore, used a ‘bumphunting’ technique adapted to WGBS data [15] to identify regions of methylation with 1.5% DNAm changes per week (corresponding to 10.5% changes across our developmental window) change in DNAm levels across adjacent CpGs, and calculated statistical significance using 1000 bootstrap-based permutations [16]. Using a cut-off of family-wise error rate (FWER) <0.2, there were 34 DMRs across prenatal development, though a less conservative cut-off (p < 0.05) identified 3,446 DMRs (Table S6, Figure S2A for DMR plot). The DMRs are similarly unevenly distributed throughout the genome as the CpGs, being far less frequent on chrX and variable among the autosomes. The DMRs on average had a width of 1696 bp (IQR: 921–2181), 75% overlap with annotated genes, and 12% overlap with annotated CpG islands. Like the CpGs, the DMRs were split between hypo- and hypermethylation, wherein 49% of the DMRs, methylation increased with age.

We also found differentially methylated sites and regions between the sexes in these prenatal samples. There were 618,978 significantly associated CpGs by sex and as expected, the vast majority (95%) of these were on the sex chromosomes. There were still 30,823 significant autosomal CpGs (Figure 2(c,d), Table S7), and while 993 of them were in regions homologous to the sex chromosomes, the majority (97%) had no homology to chrX or chrY. Among these, we found conservatively 28 DMRs with a 10% methylation change between sexes (Table S8, Figure S2B for DMR plot). Again we saw that these were not global, but regional changes in methylation, with equal numbers of DMRs being hypo- and hypermethylated in males and females. These data show that there are many differences in the DNAm landscape throughout second-trimester development, and between sexes prenatally.

Previously, many studies of brain DNA methylation have used the Illumina Infinium® HumanMethylation450 BeadChip (‘450k’) and more recent Infinium MethylationEPIC (‘850k’) microarray technologies. While these platforms can sensitively measure DNAm levels without the high coverage/sequencing requirements of WGBS, they assay a limited number of sites. To better contrast the tradeoffs between breadth (assaying more sites) and depth (assaying a given site more accurately), we compared our WGBS findings to analogous statistics calculated using these arrays on the same DNA extractions. First, using the probe coordinates alone, the 450k array does not measure DNAm levels at 35,401 (97%) of the significant age-correlated CpGs or 610,578 (98.7%) of the significant sex-correlated CpGs we identified. The 450k DNAm data validated 566 (49%) out of the 1,145 age-correlated sites that were covered. Performing the same age model on the 450k data, our WGBS data validated 279/423 (66%) of the significant results. The effects of age on methylation were generally directionally concordant between the two datasets (Figure S3). Though DMRs have much wider spans, the 450k does not cover 10 (30%) and 12 (43%) of DMRs for age and sex, respectively. The newer 850k microarray had almost twice the number of probes, yet we found that it does not measure 95% of our significant sites for age and 98% of sites for sex. Additionally, this microarray design failed to capture 24% of age DMRs and 36% of sex DMRs. Prevalent microarray platforms are therefore potentially missing a great deal of dynamic DNAm sites across the genome, and further highlight the tradeoffs between breadth, depth, cost, and complexity of generating and analysing DNAm data for epigenomic analyses.

Functional characterization of associated sites

While methylation analyses provided detailed information about changes in the methylation landscape across brain prenatal development and between the sexes, we wanted to better understand the molecular consequences of these changes. We performed gene set enrichment analyses using Gene Ontology (GO) sets to better understand in which processes the genes containing significant CpGs were involved. The top biological processes associated with genes containing age-associated CpGs were related to axon development and guidance, and regulation of neuron projection (Figure 3, Table S9). These sites were also shown to be enriched for enhancers in prenatal brain tissue as well as enriched for transcription start site regions in brain tissue from the Epigenome Roadmap project [17], implying that these areas are likely functional in prenatal brain [17] (Figure 4(a)). To better assess the putative functionality of the methylation at these CpGs and DMRs, we correlated these DNAm associations with nearby expression levels using RNA-seq data from the same cortical dissection. Among the 36,546 age-related CpGs, 59% were within an RNA-seq-measured gene (corresponding to 1,648 unique genes), and 7,983 (37%) of these CpGs had methylation levels significantly correlated to expression of the gene they lie directly within (at p < 0.05, corresponding to 653 unique genes, 40%, see Methods). Correlations were seen in both positive and negative directions, depending on the CpG, with a slight skew towards negative correlation, though we see that in CpGs that lie in promoters, correlation is mostly negative, as one would expect (Figure 4(c)). These expression-correlated CpGs did not have any differential gene set enrichment from the overall group of proximal genes. These sites were mostly in weak transcription, enhancer, and quiescent chromatin states in prenatal brain tissue (Figure 4(b)) [17]. A total of 24,679 (66%) CpGs were within an expanded window of 5 kb around genes, and 8,759 (35%) correlated to the nearest gene’s expression levels. It is important to note that this procedure only tested sites within actual genes, and only the effect on the same gene, so this did not account for potential enhancer, upstream, and trans-acting effects.

Figure 3. Gene Ontology (GO) and gene set enrichment for age-associated CpGs. The top 20 GO hits for all categories are plotted in order of significance (indicated by bar colour), showing the gene count for each GO category

Figure 4. Functional consequences of differential methylation (a) Chromatin states of age-associated CpGs in various tissues. We show the distribution of chromatin states from the 15-state model of our CpGs in various tissues from the Roadmap Epigenomics Consortium. In Roadmap prenatal brain tissues, our sites are enriched for enhancer regions, and in all brain tissues, our sites are enriched for TSS sites. (b) Chromatin state of CpGs correlated to expression in prenatal brain tissue. We see that in prenatal brain tissue, most of our age-associated and expression-correlated CpGs are in weak transcription or enhancer states. Chromatin states are defined as follows. 1: Active TSS 2: Flanking Active TSS 3: Transcription at gene 5ʹ and 3ʹ 4: Strong transcription 5: Weak transcription 6: Genic enhancers 7: Enhancers 8: ZNF genes & repeats 9: Heterochromatin 10: Bivalent/Poised TSS 11: Flanking Bivalent TSS/Enhancer 12: Bivalent Enhancer 13: Repressed PolyComb 14: Weak Repressed PolyComb 15: Quiescent/Low (c) An example of correlation between methylation level at a promoter and gene expression. Gene expression from RNA-seq is plotted against DNAm level at a correlated site, showing a decrease in gene expression with increasing methylation in example gene FEZF2

The correlation to expression was much lower for sex-related CpGs, at 5% and 9% significant correlation depending on inclusion or exclusion of sex chromosomes. The lack of correlation to expression here is unsurprising because there are almost no transcriptomic differences between the sexes in these samples. Gene set enrichment on autosomal sex-related CpGs found genes involved with synaptic transmission and signalling, regulation of GTPase activity, and the glutamate receptor signalling pathway (Table S10).

Age-related DMRs were most enriched for genes related to stem cell proliferation and various cell fate specifications (Table S11). Among the DMRs across age, 21% were significantly correlated with gene expression. The genes represented by this are OTX1, AC246817.1, CYP2E1, PLEKHH3, and DUX4L32. We also find that 33% of the DMRs overlap with one or more enhancers defined by psychENCODE [18]. There was no enrichment for any gene category among the sex-related DMRs, but it is worth noting that many of these genes coded for lncRNAs. Twenty-five per cent of sex-related DMRs were correlated to gene expression within their own gene and these genes were LINC01606, LINC01166, SPATC1L, and RFPL2. To compare this rate of correlation to the rest of the genome, we looked at non-DMRs: regions that fit bumphunter’s clustering criteria but were not differentially methylated. Correlation to gene expression overall was not different between differentially methylated regions and non-differentially methylated regions; the DMRs we found were not enriched for highly correlated regions. These results provide potential starting points for further understanding the DNA methylation landscape of the developing brain, allowing us to better understand which processes are active during normal and disordered development.

CpH methylation changes

While DNA methylation occurs exclusively in the CpG context in almost all somatic cell types, neurons in the human brain uniquely have high levels of DNA methylation in other cytosine contexts (predominantly CpA) [19]. We, therefore, investigated the potential role of CpH methylation across brain development. Only 2.4% (IQR: 2.33–2.49) of sites in a CHG context and 1.71% (IQR: 1.61–1.79) of sites in a CHH context were methylated. For comparison, around 10% of CpH sites are methylated in postnatal neurons, so it is likely as previously suggested [1] that CpH methylation accumulates beginning around the time of birth. This contrasted with CpG sites, which were predominantly methylated across the genome (Table S2). Because there are many more CpH than CpG sites in the genome, there were actually comparable numbers of methylated CpG and CpH sites (N = 13,166,520), despite the vastly different proportions. However, the DNAm levels at these sites were much lower in the CpH context (mean = 10%), likely resulting from the low proportion of mature post-mitotic neurons that contribute to this signal.

Using the same linear model as on CpGs above, no CpHs were significantly associated with age within this time period when accounting for the large multiple testing burden. However, 51 autosomal CpHs were significantly associated with sex (at FDR < 0.05, Table S12). Additionally, there were 40 associated CpHs on the sex chromosomes, though most of these were on chrY, overall very proportionally different from CpG methylation. There was no enrichment for any more specific trinucleotide context in these significant sites over the whole genome. These CpHs were less likely to be in or near genes than the CpGs, but of the 17 autosomal CpHs that were in or near genes, only 1 was significantly correlated (p < 0.05) with expression of the nearest gene within 10 kb. It is notable that the majority of these CpHs are relatively far from a gene. Additionally, the effect size of the CpH methylation at these significant sites was independent of the effect of CpGm in the surrounding area in most cases. Five of the CpHs’ effects were reduced when accounting for the methylation of the nearest CpG (Figure S4), showing that these few were dependent on the local CpG landscape, but most function independently. These loci were found both in and outside of genes and their surrounding regions. Overall, CpH methylation did not seem to be very dynamic or functional in prenatal second-trimester development, in contrast to its dynamic nature after birth [1].

Enrichments for neuropsychiatric conditions and their genetic risk

Lastly, we assessed the clinical relevance of age- and sex-associated DNA methylation changes across prenatal brain development. We therefore performed enrichment testing of the 1,648 genes represented by prenatal age-associated CpGs, and found enrichment for bipolar disorder- (p = 0.014, odds ratio (OR) = 2.15), neurodevelopmental disorder- (p = 1.9e-4, OR = 7.66), and schizophrenia (p = 9.7e-6, OR = 2.96)-associated genes [20]. The set of age-associated CpG genes was particularly enriched for autism spectrum disorder (ASD)-associated genes (p = 1.6e-40, OR = 4.14) [21]. To refine this finding, we looked to a recent exome study by Satterstrom et al. [22] which stratifies its gene set by clinical symptoms. While our genes are enriched for genes associated with ASD-dominant traits (p = 2.5e-4, OR = 4.59), they are far more strongly enriched for genes associated with neurodevelopmental delay (p = 7.7e-7, OR = 10.97). Sex-associated CpGs were also enriched for ASD-associated genes (p = 6e-29, OR = 2.9), but not for genes associated with other psychiatric disorders. We also tested enrichment for differentially expressed genes in ASD, bipolar disorder, and schizophrenia [23] and found that our age-associated genes are enriched for genes that are both up- and down-regulated in ASD and schizophrenia (Table S13). We additionally assessed whether the subset of genes that were correlated with age-associated CpG DNAm levels (N = 653) showed differential enrichment compared to the full set of 1,648 proximal genes among these brain disorders. We found no differential enrichment, suggesting that many of DNAm changes might be functional without directly influencing proximal gene expression, in line with Roadmap Epigenomics enrichment analyses above.

Our age-associated DMRs represented a greater portion of the genome than CpGs, and may further have different functional effects. To further understand what clinical phenotypes may be linked to our DMRs, we performed stratified LD score regression [24]. Our DMR sets represented a fairly small portion of the genome, but at 5.8 Mb, the less stringent age-based DMRs were wide enough to detect enrichment of BMI- and Subjective Well-Being (SWB)-associated markers (Table S14). These DMRs were also shown to be enriched for brain-linked traits overall, compared to non-brain-linked traits. Additionally, the DMRs overlap with 15 of 108 schizophrenia loci from Psychiatric Genetics Consortium 1 [25], which we are particularly interested in due to the neurodevelopmental hypothesis. These results suggest that dynamic methylation states in prenatal development may play a role in many disease and non-disease phenotypes that manifest later in life.

Discussion

Here we demonstrate that CpG methylation changes are associated with prenatal cortical development, are likely functionally relevant, and may play a role in developmental psychiatric disorders such as autism spectrum disorder (ASD).

The shifting DNAm landscape during this restricted prenatal time period – the second trimester – likely plays a role in neurogenesis of the developing human brain. Following current dogma, changes in DNAm at these critical sites presumably lead to altered gene expression which promotes the formation of new neural connections in the cortex. At this point in time, we see that it is CpGs that are the important sites in the DNAm landscape, as opposed to CpHs. This is in line with previous research which finds that most CpHm is established around birth and postnatally [1,14], later in neuron differentiation and development [26]. We also see that much – but not all – of this regulation is occurring in cis – within the site’s own gene. We do not currently understand the trans effects of our developmentally dynamic sites, but they may explain the effects of the methylation that does not correlate to its nearest gene’s expression. The fact that the age-associated CpG sites and DMRs we identified are often in enhancer regions in prenatal brain is further evidence that these sites may be regulating from a distance, reflecting the growing importance of higher-order chromatin structures in the developing brain [27]

Our findings are in line with previous findings based on microarray-derived data [10], but using WGBS allows us to assess methylation at far more loci than the more commonly used Illumina microarrays. We have shown that these microarrays do not even measure the majority of sites that are dynamic in neural development. Despite potentially less precision, we also show that WGBS has enough coverage to assess DNAm dynamics and that its findings were consistent with deeper-coverage assays. Additionally, we present novel loci due to our narrower developmental age range, which means that we can detect more subtle changes in the second trimester but also means that we do not replicate all the findings of studies done across wider timespans [10]. Confining our samples to just cortex tissue rather than whole brain tissue also allows us to detect more regionally specific changes. Substantial methylation variation in our homogenate data likely comes from profiling mixtures of neuronal cells at different maturation states in the developing cortex, which are becoming identified and profiled with increasing resolution [13]. As this cellular heterogeneity was incorporated into our methylation analyses of age, we were able to identify that most age-associated sites were not associated with cellular composition, suggesting these age affects may be shared across many neuronal cell states and not reflective of the shifting cellular landscape of the brain. Future single nuclei WGBS profiling can ideally improve our broad deconvolution-estimated cell classes, which was based on a very small number of samples in culture using microarrays [28].

We also find that there are autosomal differences in DNAm between the sexes, though by our assessment they seem less functionally significant. They may be acting more in trans which makes it more difficult to assess significance, but there are also just generally far fewer transcriptomic differences between sexes than between ages. Sex-differential DNAm later in life has been linked to psychiatric genes, which has been proposed as a mechanism for sex differences in psychiatric disorders [29], so perhaps the early changing DNAm is laying the groundwork for future effects. Like Xia et al. [29], we found that ASD-associated genes are differentially methylated between the sexes, but this does not readily translate into transcriptomic differences.

Given the consistent enrichment in our developmentally associated sites for psychiatric-linked genes, these data support the notion of early neurodevelopmental components for disorders such as schizophrenia [7,30,31]. Genes implicated in bipolar disorder, schizophrenia, and ASD are clearly dynamic at this time point in life, so dysregulation at this stage could lead to vulnerabilities to these disorders later in life. Our data also imply that BMI and subjective well-being – non-disease traits – could be linked to neural development at this stage of life.

As more data are generated, particularly through genome-scale methods like WGBS, we will be able to establish normal ranges of DNAm at all ages, which will undoubtedly provide insight into molecular dynamics in this hard-to-study period and organ, as well as give clues to where deviation from the norm is important. Unfortunately, WGBS still does not distinguish between methylation (5mC) and hydroxymethylation (5hmC), another chemical modification thought to be an important state in active demethylation. However, previous research has suggested that levels of 5hmC are very low in this developmental time period [32], suggesting that 5mC is likely driving the vast majority of our age-associated methylation differences.

Our findings are only the beginning of what may be found given the limited sample size of our study, but even here we reveal important processes. There is room for much more characterization of these epigenetic marks, but it is clear that they are worth understanding. DNA methylation serves as an exciting potential avenue to understand neural development and psychiatric disorders. There are clear and functional changes in the neuronal DNAm landscape over this important window of brain development and between sexes, and further investigation will help elucidate unknown mechanisms in the brain.

Materials and methods

Study samples

Brain tissue from these second-trimester prenatal samples was obtained via a Material Transfer Agreement with the National Institute of Child Health and Human Development Brain and Tissue Bank. All specimens were flash-frozen, then brain pH was measured and post-mortem interval (PMI, in hours) was calculated for every sample (mean = 2.85 hours). Post-mortem tissue homogenates of the dorsolateral prefrontal cortex (DLPFC) were obtained from all subjects. Samples were obtained from the developing prefrontal cortex from the dorsolateral convexity of the frontal lobe half-way between the frontal pole and temporal pole, 4 mm lateral to the central sulcus. Specimens extended from the surface of the brain to the ventricular zone. There were 10 of each male and female subjects, and 17 subjects were of African American ancestry and 3 were of European ancestry (Table S1), and the samples had no observable neuropathological insults.

Data generation

Genomic DNA was extracted from 100 mg of pulverized DLPFC tissue with the phenol-chloroform method. DNA was subjected to bisulfite conversion followed by sequencing library preparation using the TruSeq DNA methylation kit from Illumina. Lambda DNA was spiked in prior to bisulfite conversion to assess its rate, and we used 20% PhiX to better calibrate Illumina base calling on these lower complexity libraries. Resulting libraries were pooled and sequenced on an Illumina HiSeq X Ten sequencer with paired end 150bp reads (2x150bp), targeting 90Gb per sample. This corresponds to 30x coverage of the human genome as extra reads were generated to account for the addition of PhiX. These WGBS libraries were constructed and sequenced independently from the same DNAs described in Price et al., 2019 [1]

Data processing

The raw WGBS data were processed using FastQC to control for quality of reads [33], Trim Galore to trim reads and remove adapter content [34], Arioc for alignment to the GRCh38.p12 genome (obtained from ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.27_GRCh38.p12/GCA_000001405.27_GRCh38.p12_assembly_structure/Primary_Assembly/assembled_chromosomes/) [35], duplicate alignments were removed with SAMBLASTER [36], and the Bismark methylation extractor to extract methylation data from the sequencing data [37]. We then used the bsseq R/Bioconductor package (v1.18) to process and combine the DNA methylation proportions across the samples for all further manipulation and analysis [15]. After initial data metrics were calculated, the methylation data were smoothed using BSmooth for modelling. CpGs were filtered to those that had ≥3 coverage in all samples, and CpHs were filtered to those that had ≥3 coverage and non-zero methylation in at least half (≥10) of the samples.

Comparison to 450k

In order to compare WGBS methylation levels to 450k methylation levels, we used data from the same samples using the two methods (from GEO: GSE74193), and compared their methylation levels at the same sites both graphically (Figure 1d) and assessing the mean differences and root-mean-square deviation (Figure 1c). We then validated our model’s findings by applying the same model to the microarray data to overlapping sites and considered loci significant at FDR < 0.05 and validated if its association had p < 0.05 in the validating set (Figure S3).

Gene set enrichment

We annotated our data using Gencode v. 29 on hg38. We performed gene set enrichment using clusterProfiler (v3.12) [38] for Gene Ontology [50] sets with a p-value cut-off of 0.01 and q-value cut-off of 0.05. We used SFARI 2.0, accessed on 5/16/2019 [21] for an autism spectrum disorder gene set, and a set of clinical gene sets defined by Birnbaum et al. for other neuropsychiatric and neurodevelopmental disorders: four schizophrenia sets based on copy number variants, GWAS, and meta-analysis, two ASD sets based on GWAS and copy number variants, and one set each for neurodegenerative disorders, neurodevelopmental disorders, intellectual disability, and bipolar disorder [20]. We also tested for enrichment using an ASD-exome gene set [22], differentially expressed gene sets [23,39], and transcriptome-wide association study signals [23,40]. Enrichment was calculated on a background of genes expressed in our samples, to avoid brain bias. We performed LD score regression as described by Finucane et al. [24] and with data from recent GWAS [41,42].

Data analysis

We performed cellular deconvolution on our WGBS data using microarray-based reference datasets, which was adapted from Jaffe et al., 2015 [43] and Jaffe et al., 2016 [9]. Briefly, this approach uses 5 different cell classes from Illumina 450k microarray data: (1) ESCs, (2) NPCs derived from ESCs, and (3) dopamine neurons derived from ESCs, from Kim et al. [28] and (4) neurons and (5) non-neurons obtained by NeuN+ sorting from adult DLPFC from Kaminsky et al., 2013 [12]. Importantly, these five cell types were estimated in a single deconvolution model, and we previously showed in the above manuscripts that DNAm and gene expression in prenatal samples was strongly associated with the estimated ESC and NPC cell fractions. Here, we first filtered those sites on the Illumina 450k that were well covered in our WGBS data (N = 324,225), and then constructed the deconvolution model (CpGs of interest and corresponding model parameters) following the strategy of Jaffe & Irizarry 2014 [44]. We then applied this model to our WGBS data at these sites to estimate the fraction of each of the five cell types in each sample.

Differential methylation analyses for age were performed using linear regression modelling, accounting for sex, age, and estimated ESC fraction. The association between sex and DNAm was extracted from this same model, those significant sex-associated CpGs (and their surrounding sequence, 150 bp window) were re-aligned to the genome with Bowtie2 [45] to check for homology to ensure that detected autosomal regions were not potentially misaligned homologous regions from the sex chromosomes. The regression analyses above were formed using limma (v3.30) [46,47] which employed empirical Bayes and returned moderated T-statistics, which were used to calculate p-values and estimate the false discovery rate (FDR, via Benjamini-Hochberg approach [48]).

RNA-seq comparison analysis

The homogenate RNA-seq samples were also part of a larger study of RNA-seq data from homogenate DLPFC tissue (BrainSeq Consortium Phase I) [49]. We compared a site’s methylation or mean methylation (for DMRs) to the expression level from RNA-seq using a linear model adjusting for the strongest covariates, including sex, age, and estimated embryonic stem cell (ESC) fraction.

Authors’ contributions

KAPM, AJP, JEK, and AEJ conceptualized the project and methodology. KAPM and AEJ investigated and analyzed the data and wrote the paper. LCT, ASS, RW, and NJE provided software support. AJP, RT, TMH, and RW curated data. JEK, RT, and TMH provided resources. AJP, LCT, AEJ, and DRW reviewed and edited the paper. DRW provided supervision, and AEJ was the principal investigator and oversaw experimental design and analysis. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Postmortem brain tissue was obtained via a Material Transfer Agreement with the National Institute of Child Health and Human Development Brain and Tissue Bank which collected the tissue under an IRB-approved protocol and consent of the next-of-kin (https://www.medschool.umaryland.edu/btbank/).

Acknowledgments

The authors thank the UMB Brain Bank at the Department of Pediatrics in the University of Maryland School of Medicine for the tissue provided. This project was supported by the Lieber Institute for Brain Development and by NIH grants R21MH102791, R01MH112751, and T32GM781437. Finally, we are indebted to the generosity of the families of the decedents, who donated the brain tissue used in these studies.

Availability of data and materials

Raw and processed nucleic acid sequencing data generated to support the findings of this study are part of the PsychENCODE Consortium and the Brainseq Consortium data releases. Specifically, WGBS data have been deposited at www.Synapse.org along with the other PsychENCODE data, under the accession code syn5842535. The homogenate RNA-seq samples were also part of a larger study of RNA-seq data from homogenate DLPFC tissue (BrainSeq Consortium Phase I) [50], which was also deposited at www.synapse.org and summarized in http://eqtl.brainseq.org/phase1. The processed, homogenate RNA-seq data for this study have additionally been deposited via Globus under the jhpce#brainepi-cellsorted collection at the following location: http://research.libd.org/globus/jhpce_brainepi-cellsorted/. NeuN-sorted RNA-seq data were originally published as part of phase II of the Brainseq Consortium (http://eqtl.brainseq.org/phase2/) and have also been deposited via Globus under the jhpce#brainepi-polyA collection at the following location: http://research.libd.org/globus/jhpce_brainepi-polyA/. Publicly available data reprocessed in support of the conclusions in this work were downloaded from the Gene Expression Omnibus under GEO accession GSE47966 and GSE74193.

Disclosure statement

The authors declare no competing interests.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This project was supported by The Lieber Institute for Brain Development and by NIH grants R21MH102791, R01MH112751, and T32GM781437; National Institutes of Health [R21MH102791]; National Institutes of Health [T32GM781437]; National Institutes of Health [R01MH112751].