Adaptation of the targeted capture Methyl-Seq platform for the mouse genome identifies novel tissue-specific DNA methylation patterns of genes involved in neurodevelopment

Abstract

Methyl-Seq was recently developed as a targeted approach to assess DNA methylation (DNAm) at a genome-wide level in human. We adapted it for mouse and sought to examine DNAm differences across liver and 2 brain regions: cortex and hippocampus. A custom hybridization array was designed to isolate 99 Mb of CpG islands, shores, shelves, and regulatory elements in the mouse genome. This was followed by bisulfite conversion and sequencing on the Illumina HiSeq2000. The majority of differentially methylated cytosines (DMCs) were present at greater than expected frequency in introns, intergenic regions, near CpG islands, and transcriptional enhancers. Liver-specific enhancers were observed to be methylated in cortex, while cortex specific enhancers were methylated in the liver. Interestingly, commonly shared enhancers were differentially methylated between the liver and cortex. Gene ontology and pathway analysis showed that genes that were hypomethylated in the cortex and hippocampus were enriched for neuronal components and neuronal function. In contrast, genes that were hypomethylated in the liver were enriched for cellular components important for liver function. Bisulfite-pyrosequencing validation of 75 DMCs from 19 different loci showed a correlation of r = 0.87 with Methyl-Seq data. We also identified genes involved in neurodevelopment that were not previously reported to be differentially methylated across brain regions. This platform constitutes a valuable tool for future genome-wide studies involving mouse models of disease.

Keywords:

• brain
• DNA methylation
• epigenetics
• genome-wide
• methylation array

Abbreviations

Apcdd1	=	Adenomatous Polyposis Coli Down-Regulated 1
ChIP	=	Chromatin immunoprecipitation
DMCs	=	Differentially methylated cytosines (DMCs)
DMRs	=	Differentially methylated regions
DNAm	=	DNA methylation
FDR	=	False discovery rate
GFAP	=	Glial fibrillary acidic protein
Gb	=	Gigabases
GO	=	Gene ontology
H3K4me1	=	Histone marks histone 3 lysine 4 monomethylation
H3K27ac	=	Histone 3 lysine 27 acetylation
KEGG	=	Kyoto Encyclopedia of Genes and Genomes
MAP	=	Mitogen activated protein
Msx1	=	msh homeobox1
PAVIS	=	Peak Annotation and Visualization
RV	=	Range of variation
TFBS	=	Transcription factor binding sites
UTR	=	Untranslated regions

Introduction

DNA methylation (DNAm) is an epigenetic mechanism that has been observed to contribute to a variety of phenotypes. Studies have noted its role in development¹ and its contribution to cell^2,3 and tissue-specific^4,5 differentiation. DNAm patterns have also been observed to vary between individuals.^6,7 Importantly, DNAm is not static, but can be altered by environmental influences, and thus provides a mechanism through which external stimuli or exposures can affect phenotypes. For example, studies on early life stressors have observed that changes to DNAm patterns in genes related to the hypothalamic-pituitary-adrenal axis, which is involved in the physiological stress response, are associated with anxiety and mood disorders.^8,9 Further, abnormalities in the DNAm pattern can also contribute to a variety of other diseases, including cancer and imprinting disorders, such as Prader-Willi syndrome and Angelman syndrome.¹⁰ Due to its involvement in both development and disease, a variety of methods have been developed to study DNAm at a genome-wide level.

Currently, 3 different approaches are commonly used in conjunction with next-generation sequencing platforms to investigate genome-wide DNAm patterns. These approaches rely on either protein-affinity enrichment of methylated regions, bisulfite conversion, or methylated cytosine-sensitive restriction enzymes. Although widely used, these methods suffer from several limitations. For example, protein affinity-based enrichment of methylated cytosine approaches have been reported to differentially enrich for regions by varying CpG density, and they do not provide base-pair resolution of methylated cytosines.^11,12 In contrast, methods using methylated cytosine-sensitive restriction enzymes can provide base-pair resolution, but enrich for CpG dense regions.¹¹ Although these 2 approaches can be used complementarily, limited starting material and increased cost of sequencing may be prohibitive. An alternative unbiased method to detecting methylated CpGs at base-pair resolution is whole genome bisulfite sequencing. However, for this approach, obtaining adequate sequencing depth requires a substantial elevation in sequencing cost. In addition, 70–80% of the sequenced reads do not provide useful information on DNA methylation.¹³ Given these limitations, a method that provides cost-effective detection of CpG methylation at base-pair resolution without bias of CpG dense and poor regions would have advantages for the study of genome-wide DNAm. Genome-wide targeted capture technology offers these advantages.

The genome-wide targeted capture approach uses DNA or RNA baits that contain complementary sequences of targeted regions to enrich these candidate regions for next-generation sequencing. This approach has been widely used in exome sequencing studies^14,15 to identify rare variations. It has more recently been adapted to investigate DNAm at base-pair resolution through bisulfite conversion of the human genome.^16,17 As baits can be designed against regions of interest, this eliminates any bias for CpG dense or poor regions while enabling sequencing reactions to be focused on regions of interest, thus improving yield per sequencing cost. Despite these benefits, only a small number of studies to date have used this approach in humans,^16,17 and no study has yet used it in mouse. This study describes the adaptation of Methyl-Seq for the mouse methylome and demonstrates its ability to reliably detect DNAm by investigating differences between tissue types. It has also revealed novel tissue-specific differences in DNAm patterns of genes involved in neurodevelopment. This is currently the only bait-based assay for the mouse methylome. This study provides the foundation for future use of this assay in studying genome-wide DNAm signatures in mouse models of disease.

Results

Single-base resolution maps of 5-methylcytosines from targeted bisulfite sequencing of liver and cortex and hippocampus

Single base-pair resolution maps of 5-methylcytosine was generated by bisulfite-sequencing of targeted captured DNA from 3 different tissues (liver, cortex, and hippocampal) each performed in triplicate. Genomic regions captured by this assay are summarized in Table 1. This generated an average of ∼81 million, 77 million, and 77 million reads for cortex, hippocampus, and liver, respectively, which uniquely aligned to the bisulfite-converted mouse reference genome (mm9/NCBI37). Eighty-five percent (cortex), 67% (hippocampus) and 70% (liver) of aligned reads remained after deduplication, which covered ∼75% of targeted regions. Approximately 25% of targeted regions had zero coverage, suggesting that baits were unable to isolate these regions. From these filtered reads, 6.4 gigabases (Gb) (cortex), 4.9 Gb (hippocampus) and 5.0 Gb (liver) were present. The mean coverage of targeted regions was 54× (cortex), 42× (hippocampus) and 44× (liver). About 73% of targeted bases were observed to have coverage of at least 10×. As a measure of efficiency of the targeted capture assay, targeted regions were observed to be ∼20-fold more abundant than genomic background for all tissues. An overview of these results is described in Table 2.

Table 1. Mouse methylome target capture design. There was approximately 10 Mb of target overlap between data repositories resulting in a final hybridization capture array targeting 99 Mb across the murine genome. TFBS: transcription factor binding sites.

Site Classification	Number of Targets	Total Bases Covered (bps)
CpG Islands	16,027	10,512,276
Tissue Specific DMRs	33,456	10,452,692
Ensembl Regulatory Features	171,796
-CpG shores and shelves (±4kb )
-Dnase I Hypersensitive sites		91,799,015
-Histone Modifications
-TFBS
-Polymerase
Open Regulatory Annotation (ORegAnno)	14,951
-Promoters
-Enhancers		9,983,957
-TFBS
-Regulatory Polymorphisms
Total Merged Intervals	201,788	109,147,819

Table 2. Overview of Methyl-Seq data for different tissues. Values shown are the average of 3 independent samples from each tissue. Values shown were obtained from Picard HsMetrics and FastQC analysis. Filtering mentioned in this table involves deduplication of aligned reads.

	Cortex	Hippocampus	Liver
Total aligned reads	81 million	77 million	77 million
Remaining aligned reads after filtering	69 million	51 million	54 million
% of total aligned reads after filtering	85	67	70
% of duplicated sequence from filtered reads	5	5	5
Number of aligned bases after filtering	6.4 billion	4.9 billion	5.0 billion
% of aligned bases in target region after filtering	79	80	82
Mean target coverage from filtered reads	54	43	44
% of target bases with ≥ 10×	75	71	72
Fold enrichment of targeted region	20	20	20
% of targets with zero coverage	23	24	24

We report here on our evaluation of tissue-specific differentially methylated cytosines (DMCs) in the CpG context. A total of ∼1.3 million CpGs were commonly captured in all tissues with a minimum coverage of 10× and good quality base calls (< Q30). Correlation analysis showed a tight reproducibility of methylation levels between replicates of cortex (Pearson's r = 0.97, P < 2.2 × 10⁻¹⁶) (Fig. 1A), hippocampus (Pearson's r = 0.97, P < 2.2 × 10⁻¹⁶) (Fig. 1B), and liver (Pearson's r = 0.96, P < 2.2 × 10⁻¹⁶) (Fig. 1C). A high correlation in methylation was observed between cortex and hippocampus (Pearson's r = 0.98, P < 2.2 × 10⁻¹⁶) (Fig. 2A). In contrast, reduced correlation was observed between brain tissues and the liver (Pearson's r = 0.84, P < 2.2 × 10⁻¹⁶) (Fig. 2B-C). About 85% of CpGs showed some degree of methylation in both brain and liver tissues. Between tissues, CpG methylation was observed to be significantly higher in the cortex (42% ± 35%, mean ± SD) than hippocampus (40% ± 35%, mean ± SD) (P = 1.15 × 10⁻¹¹). CpG methylation in both brain tissues were higher than in the liver (38% ± 34%, mean ± SD) (P < 2.2 × 10⁻¹⁶). We also analyzed tissue differences using hierarchical clustering and principle component analysis approaches. These, seen in Figure S1 (A–B), further demonstrated distinct tissue-specific clustering of methylation patterns. A histogram reflecting the variance within each tissue type is provided in Figure S2 (A–C).

Figure 1. Reproducibility of cortex, hippocampus, and liver replicates. (A) Cortex. (B) Hippocampus. (C) Liver. Figures shown are scatter plots of comparisons between DNA methylation levels of tissue replicates. Data shown are a representative of one of the 3 comparisons between replicates of each tissue.

Figure 2. Comparison of DNA methylation level between different tissues. (A) Hippocampus vs. Cortex. (B) Cortex vs. Liver (C) Hippocampus vs. Liver. Figures shown are scatter plots between different tissues using average DNA methylation levels from triplicates of each tissue tissue.

To determine the technical variation for this assay, targeted capture followed by bisulfite-sequencing was performed on DNA derived from cortex and liver in technical triplicates. Sequencing data was aligned to the bisulfite-converted mouse genome and deduplicated as mentioned above. Methylated cytosines in targeted regions with a minimum coverage of 10× and of good base quality (<Q30) across triplicates were used to determine the range of variation (RV) between technical triplicates. This was calculated by determining the difference between the maximum and minimum methylation level for every methylated CpG. The RV was similar across the different tissues (10% ± 10%, median ± IQR). An RV in the 95th percentile, reflecting a range of 27%, was used as a threshold for detecting tissue-specific methylation differences as this removes 95% of technical variation from the assay. Using these screening criteria, tissue-specific DMCs in the CpG context were identified between hippocampus and liver (N = 186,213), cortex and liver (N = 197,519), and cortex and hippocampus (N = 14,519) respectively. The absolute methylation difference for DMCs in the CpG context was not significantly different between the cortex vs. liver (44% ± 13%, mean ± SD) as compared to the hippocampus vs. liver (44% ± 13%, mean ± SD) after Bonferroni correction for multiple testing. However, the absolute methylation difference was significantly higher (P < 2.2 × 10⁻¹⁶) for these 2 comparisons as compared to the cortex vs. hippocampus (33% ± 6%, mean ± SD).

Distribution of differentially methylated cytosines in genomic and CpG contexts

Previous studies have shown that the distribution of DMCs between cell and tissue types can vary between genomic contexts.^4,18 To determine the genomic distribution of DMCs among these regions, DMCs were assigned to their targeted capture baits that had been mapped to these functional genomic contexts. For all tissue comparisons, DMCs in the CpG context were present in significantly greater numbers than expected in the introns and intergenic regions downstream of genes (≥1 bp downstream of transcriptional termination site as previously described¹⁹) (Fig. 3A-C). For 3′ untranslated regions (UTR), DMCs in the CpG context were present in significantly greater numbers than expected in the brain vs. liver comparisons, but not in the cortex vs. hippocampus comparison. In contrast, DMCs in promoters, exons, and 5′ UTR were significantly fewer than expected in all tissue comparisons (Fig. 3A-C). Likewise, DMCs in upstream intergenic (>2,000 bps upstream of transcriptional start site) regions were present in significantly smaller numbers than expected in the brain vs. liver comparisons, but not in the cortex vs. hippocampus comparison which showed no significant difference between observed and expected numbers (Fig. 3A-C). Similarly, the distribution of DMCs was assessed according to previously defined CpG contexts, namely CpG islands, CpG shores, CpG shelves, and CpG open sea.^18,20 The fewest DMCs were observed in CpG islands for all tissue comparisons, and these numbers were significantly smaller than expected (Fig. 4A-C). In contrast, more DMCs were present in CpG contexts outside of CpG islands such as CpG shores, CpG shelves and CpG open sea, and these numbers were significantly greater than expected (Fig. 4A-C). However, between cortex and hippocampus, the number of DMCs in CpG shores was significantly less than expected (Fig. 4C).

Figure 3. Distribution of DMCs in the CpG context across genomic annotations. Chi-Square was performed to determine if the number of observed DMCs deviated from the expected number of DMCs for each genomic category. (A) Hippocampus vs. Liver, (B) Cortex vs. Liver, and (C) Hippocampus vs. Cortex. ***: P < 0.001, N.S.: not significant. P-values remained significant after Sidak's post-hoc multiple comparison test.

Figure 4. Distribution of DMCs across CpG contexts. Chi-Square was performed to determine if the number of observed DMCs deviated from the expected number of DMCs for each CpG context. (A) Hippocampus vs. Liver, (B) Cortex vs. Liver and (C) Hippocampus vs. Cortex. ***: P < 0.0001. P-values remained significant after Sidak's post-hoc multiple comparison test.

Functional annotation of differentially methylated cytosines between liver and brain tissues

To determine how tissue-specific DMCs might contribute to cellular differences between brain and liver tissues, DMCs in the CpG contexts were assigned to the nearest gene and analyzed by gene ontology (GO) term enrichment of cellular components. After correcting the enrichment P-values for multiple testing, genes that were more methylated in the CpG context in the liver compared to cortex and hippocampus had significant enrichment of GO terms for neuronal components, as would be expected given the correlation of increased methylation with reduced expression and thus functional deactivation. These included (only cortex vs. liver result is shown for simplicity): neuronal projection (P = 5.9 × 10⁻¹⁷), microtubule cytoskeleton (P = 9.5 × 10⁻⁹), synapse part (P = 1.1 × 10⁻¹⁰), dendrite (P = 5.8 × 10⁻¹⁰), postsynaptic density (P = 3.3 × 10⁻⁷), and axon (P = 1.3 × 10⁻⁵). In line with these cellular differences, molecular signaling pathways were also observed to be different between the liver and brain tissues. Pathway analysis showed that genes which were more methylated in the CpG context in the liver were significantly enriched in signaling pathways known to be involved in neuronal function including the mitogen activated protein (MAP) kinase signaling pathway (P = 3.6 × 10⁻¹⁰), the Wnt signaling pathway (P = 1.7 × 10⁻⁹), axon guidance (P = 1.7 × 10⁻⁸), the Erb signaling pathway (P = 1.7 × 10⁻⁶), and the neurotrophin pathway (P = 1.9 × 10⁻⁵). These differences remained significant after multiple corrections.

Distribution of differentially methylated cytosines between enhancer regions

Many studies have shown that enhancers play an important role in regulating tissue-specific gene expression. To determine if differential methylation patterns at enhancer regions contribute to phenotypic differences, the distribution of DMCs were mapped to enhancers shared between the cortex and liver, cortex-specific enhancers, liver-specific enhancers, and regions outside of these categories. Enhancers were identified using the histone marks histone 3 lysine 4 mono-methylation (H3K4me1) and histone 3 lysine 27 acetylation (H3K27ac) > 2 kb from the transcriptional start site.^21,22 Histone marks based on cortex and liver of adult C57/Bl6 mice were retrieved from the UCSC Genome Browser. DMCs were present in significantly greater numbers than expected in enhancers shared between the cortex and liver (P = 2.6 × 10⁻¹⁵) (Fig. 5A), cortex-specific enhancers (P < 1.0 × 10⁻¹⁵) (Fig. 5B) and liver-specific enhancers (P = 1.0 × 10⁻¹⁵) (Fig. 5C). In contrast, DMCs were present in significantly lower numbers than expected in regions outside of these enhancer categories (Fig. 5D). Of those DMCs mapped to cortex-specific enhancers, 86% were more highly methylated in the liver, while 14% were more highly methylated in the cortex. Similarly, of those DMCs mapped to liver-specific enhancers, 96% were more highly methylated in the cortex, while 4% were more highly methylated in the liver. Of those present in enhancers shared between the cortex and liver, 68% were more highly methylated in the cortex while 32% were more highly methylated in the liver. We did not include the hippocampus for this analysis, as the data for H3K4me1 and H3K27ac marks for adult C57/Bl6 hippocampus is unavailable.

Figure 5. Distribution of DMCs across enhancer regions. Chi-Square was performed to determine if the number of observed DMCs deviated from the expected number of DMCs for each enhancer context. (A) Enhancers shared between cortex and liver. (B) Enhancers present in the cortex. (C) Enhancers present in the liver. (D) Enhancers in other regions. ***: P < 0.0001. P-values remained significant after Sidak's post-hoc multiple comparison test.

Validation of differentially methylated cytosines by bisulfite-pyrosequencing

To evaluate whether Methyl-Seq reliably detects DNAm differences between tissues, a selection of methylated CpGs that were the most highly statistically significant were interrogated by bisulfite-pyrosequencing using Methyl-Seq samples. CpGs investigated in this study were selected from regions containing clusters of DMCs, i.e., differentially methylated regions (DMRs), ranked by number of DMCs (Table 3). As expected, comparison of these DMRs to histone modifications (H3K4me1, H3K4me3 and H3K27ac) and DNase I hypersensitive sites annotated in the UCSC Genome Browser showed that every one of the 19 regions contained hallmarks of at least on of these regulatory genomic features.²³ A good correlation in DNAm differences between tissues was observed between Methyl-Seq and bisulfite-pyrosequencing data (Pearson's r = 0.84, P < 2.2 × 10⁻¹⁶). In a separate experiment, another cohort of mouse tissues was used to validate tissue-specific DMCs detected by Methyl-Seq. From a total of 75 tissue-specific DMCs, ∼86% showed significant tissue-specific DNAm differences; all in the same direction as observed by Methyl-Seq (Table 4A-C), with an overall correlation of 0.87 (P < 2.2 × 10⁻¹⁶) when comparing tissue-specific DNAm differences between the 2 platforms. There was, however, variation in the magnitude of tissue-specific differences between the 2 platforms. For example, tissue-specific methylation differences observed in Gm13345 and Zfp2865 varied by ∼8% between Methyl-Seq and bisulfite-pyrosequencing detection. Tissue-specific methylation differences observed in Sox9 showed variation of 17%–24% between Methyl-Seq and bisulfite-pyrosequencing. The average variation detected between the 2 approaches was 20% ± 12% (mean ± SD). To determine the relationship between methylation levels and gene expression levels for validated differentially methylated genes, expression data from the mouse cortex, hippocampus and liver from GNF Expression Atlas 2 was obtained using the UCSC Genome Browser. Expression data was available for 15 of the 17 validated genes. Of the 15 genes, 11 genes showed reduced gene expression in tissues that had increased methylation. In contrast, 4 genes (Neat1, Nr1d1, Espn, Potla) did not follow this trend. Neat1 did not show any difference in gene expression between the cortex and liver even though the cortex was more highly methylated than the liver. Similarly, Pot1a did not show any difference in gene expression between the cortex and hippocampus even though the hippocampus was more highly methylated than the cortex. Nr1d1 and Espn showed higher expression in the liver than cortex although their methylation level was higher in the liver.

Table 3. Differentially methylated regions (DMRs) identified from Methyl-Seq for bisulfite-pyrosequencing validation. Differentially methylated cytosines (DMCs) were clustered as DMRs and ranked by number of DMCs; DMC density is calculated by the number of DMCs divided by the length of the DMR.

	Chr	Start	End	Dist. to TSS	Gene	Number of DMCs.	DMC density	Absolute Meth. Diff.	P-value	FDR
Cortex vs. Liver	chr19	5842873	5846952	1954	Neat1	41	1.00E-02	32	5.62E-09	5.44E-08
	chr2	172980623	172981822	2042	Pck1	22	1.83E-02	49	4.76E-11	6.09E-10
	chr6	117122068	117123267	3475	Cxcl12	17	1.42E-02	59	7.05E-32	1.12E-29
	chr11	98642984	98643343	−6338	Nr1d1	7	1.94E-02	42	5.90E-10	6.56E-09
	chr15	59482973	59485852	2767	Trib1	38	1.32E-02	41	4.23E-17	1.12E-15
Hippocampus vs. Liver	chr10	79594075	79595077	2350	Dos	42	4.19E-02	36	6.11E-09	5.04E-08
	chr4	151499194	151501113	1922	Espn	63	3.28E-02	36	2.91E-08	2.18E-07
	chr11	112645778	112646857	0	Sox9	31	2.87E-02	35	1.66E-08	1.29E-07
	chr3	55586231	55588510	0	Mab21l1	52	2.28E-02	38	7.02E-10	6.57E-09
	chr16	18585137	18586816	247	Tbx1	19	1.13E-02	33	4.33E-08	3.17E-07
Cortex vs. Hippocampus	chr11	48639149	48639868	0	Trim7	13	1.81E-02	25	7.68E-08	1.52E-05
	chr2	52427056	52427535	−12970	Cacnb4	10	2.08E-02	41	2.99E-29	5.77E-25
	chr6	25758574	25760493	0	Pot1a	8	4.17E-03	10	2.77E-03	5.38E-02
	chr9	72704496	72704977	49417	Prtg	15	3.11E-02	31	3.80E-09	1.27E-06
	chr18	63083911	63085175	1930	Apcdd1	12	9.49E-03	30	3.64E-11	2.39E-08
	chr5	38217039	38217694	−1218	Msx1	14	2.13E-02	29	3.22E-10	1.60E-07
	chr2	35346662	35347111	0	Gm13445	14	3.11E-02	32	6.43E-07	8.87E-05
	chr13	74671394	74671993	−51004	Zfp825	7	1.17E-02	25	9.09E-06	7.20E-04
	chr8	106832961	106833918	0	Cmtm1	10	1.04E-02	15	2.82E-04	1.02E-02

Table 4. Comparison of Methyl-Seq with bisulfite-pyrosequencing validation. Bisulfite pyrosequencing was performed on a separate cohort of liver and brain tissues at CpG dinucleotide locations described above. (A) Cortex vs Liver, (B) Hippocampus vs Liver and (C) Cortex vs. Hippocampus data shown of average percentage methylation are mean ± SEM of methylation level detected at Cs in the CpG context at the corresponding location in cortical and hippocampal tissues. Methyl-Seq, N = 3. Bisulfite-pyrosequencing, N = 8. P-values for Methyl-Seq and bisulfite-pyrosequencing data were corrected for multiple testing by the sliding linear model method and Holm-Šídák test, respectively.

(A)
		Methyl-Seq				Bisulfite-pyrosequencing
Gene	Location	Average % methylation (Cortex)	Average % methylation (Liver)	Average % methylation difference	P-value	Average % methylation (Cortex)	Average % methylation (Liver)	Average % methylation difference	P-value
Neat1	chr19: 5843973	66.0 ± 3.8	7.3 ± 2.2	58.7	<1.0 × 10^–15	97.6 ± 0.4	60.6 ± 13.3	37.0	0.014
	chr19: 5843978	80.7 ± 0.9	7.3 ± 2.7	73.4	<1.0 × 10^–15	98.2 ± 0.2	75.8 ± 7.0	22.5	N.S.
	chr19: 5843984	51.0 ± 2.6	3.7 ± 1.2	47.3	<1.0 × 10^–15	98.8 ± 0.4	79.5 ± 6.7	19.3	N.S.
Pck1	chr2: 172980725	87.3 ± 2.2	21.7 ± 10.9	65.6	1.9 × 10^–13	97.6 ± 1.1	44.9 ± 0.8	52.8	<0.0001
	chr2: 172980736	78.0 ± 6.0	24.0 ± 11.5	54.0	3.9 × 10^–10	86.8 ± 1.0	41.1 ± 0.8	45.6	<0.0001
	chr2: 172980741	84.7 ± 6.3	24.0 ± 9.2	60.7	4.6 × 10^–14	89.4 ± 0.7	40.0 ± 0.5	49.4	<0.0001
	chr2: 172980808	92.7 ± 0.7	44.3 ± 2.2	48.4	<1 × 10^–15	96.9 ± 0.1	64.4 ± 1.3	32.5	<0.0001
	chr2: 172980818	97.3 ± 1.5	47.3 ± 8.4	50.0	2.4 × 10^–15	97.0 ± 0.8	70.7 ± 1.2	26.4	<0.0001
	chr2: 172980825	89.7 ± 1.7	39.3 ± 3.5	50.4	1.8 × 10^–14	95.6 ± 0.7	59.4 ± 1.5	36.3	<0.0001
Cxcl12	chr6: 117122940	90.3 ± 2.0	20.0 ± 4.7	70.3	<1 × 10^–15	95.3 ± 0.25	37.6 ± 0.8	57.6	<0.0001
	chr6: 117122942	91.0 ± 2.1	24.7 ± 3.9	66.3	<1 × 10^–15	93.0 ± 0.3	33.3 ± 0.6	59.8	<0.0001
Nr1d1	chr11: 98643223	21.3 ± 0.3	58.7 ± 4.9	−37.4	5.4 × 10^–8	41.5 ± 3.1	85.9 ± 0.3	−44.4	<0.0001
	chr11: 98643236	7.0 ± 2.1	62.3 ± 1.3	−55.3	<1 × 10^–15	51.9 ± 0.8	86.1 ± 0.1	−34.3	<0.0001
Trib1	chr15: 59485570	87.7 ± 2.0	26.7 ± 11.2	61.0	<1 × 10^–15	70.9 ± 0.5	33.5 ± 0.6	37.4	<0.0001
	chr15: 59485576	83.7 ± 2.0	18.7 ± 4.8	65.0	<1 × 10^–15	83.1 ± 0.5	32.4 ± 0.7	50.8	<0.0001
	chr15: 59485587	67.0 ± 2.1	23.7 ± 3.7	43.3	<1 × 10^–15	86.4 ± 2.0	35.1 ± 0.7	51.2	<0.0001
(B)
		Methyl-Seq				Bisulfite-pyrosequencing
Gene	Location	Average % methylation (Hippocampus)	Average % methylation (Liver)	Average % methylation difference	P-value	Average % methylation (Hippocampus)	Average % methylation (Liver)	Average % methylation difference	P-value
Dos	chr10: 79594501	40.7 ± 3.0	84.3 ± 4.1	−43.6	<1.0 × 10^–15	87.4 ± 0.4	96.0 ± 0.2	−8.6	<0.0001
	chr10: 79594505	45.7 ± 4.4	89.0 ± 2.1	−43.3	<1.0 × 10^–15	98.9 ± 0.8	99.0 ± 0.8	−0.1	N.S.
	chr10: 79594513	27.0 ± 1.2	87.3 ± 1.3	−60.3	<1.0 × 10^–15	79.4 ± 1.3	89.1 ± 1.0	−9.8	<0.0001
Espn	chr4: 151499268	6.7 ± 1.5	44.0 ± 6.0	−37.3	<1.0 × 10^–15	19.6 ± 1.2	84.1 ± 1.0	−64.6	<0.0001
	chr4: 151499277	1.7 ± 1.7	37.0 ± 3.6	−35.3	<1.0 × 10^–15	24.9 ± 1.9	90.1 ± 0.6	−65.3	<0.0001
	chr4: 151499288	8.7 ± 2.0	52.0 ± 2.6	−43.3	<1.0 × 10^–15	8.9 ± 1.1	62.3 ± 0.9	−53.4	<0.0001
	chr4: 151499290	6.0 ± 1.5	51.3 ± 2.2	−45.3	<1.0 × 10^–15	16.5 ± 1.8	62.9 ± 1.1	−46.4	<0.0001
Sox9	chr11: 112645850	13.0 ± 3.5	53.0 ± 14.2	−40.0	2.3 × 10^–08	65.3 ± 2.6	86.1 ± 0.7	−20.9	<0.0001
	chr11: 112645852	17.3 ± 0.3	60.0 ± 6.5	−42.7	5.1 × 10^–9	69.8 ± 1.7	88.6 ± 1.0	−18.9	<0.0001
	chr11: 112645861	8.0 ± 2.5	43.7 ± 4.3	−35.7	2.7 × 10^–09	56.0 ± 1.7	72.8 ± 1.1	−16.8	<0.0001
	chr11: 112645863	10.7 ± 2.3	52.0 ± 3.5	−41.3	10.0 × 10^–11	72.1 ± 1.6	88.9 ± 0.8	−16.8	<0.0001
	chr11: 112645869	11.3 ± 2.2	55.3 ± 3.7	−44.0	6.3 × 10^–12	63.6 ± 3.2	82.9 ± 1.2	−19.3	<0.0001
Mab21l1	chr3: 55587792	70.7 ± 3.3	12.7 ± 5.2	58.0	<1.0 × 10^–15	63.8 ± 0.9	20.0 ± 0.6	43.8	<0.0001
(B)
		Methyl-Seq				Bisulfite-pyrosequencing
Gene	Location	Average % methylation (Hippocampus)	Average % methylation (Liver)	Average % methylation difference	P-value	Average % methylation (Hippocampus)	Average % methylation (Liver)	Average % methylation difference	P-value
	chr3: 55587807	70.0 ± 4.2	14.3 ± 11.8	55.7	<1.0 × 10^–15	61.6 ± 1.5	23.1 ± 0.6	38.5	<0.0001
	chr3: 55587811	80.3 ± 0.7	21.0 ± 5.0	59.3	<1.0 × 10^–15	68.9 ± 1.0	23.8 ± 0.8	45.1	<0.0001
	chr3: 55587819	84.3 ± 2.7	21.0 ± 8.0	63.3	<1.0 × 10^–15	66.1 ± 1.4	27.4 ± 0.7	38.8	<0.0001
Tbx1	chr16: 18585838	81.0 ± 1.2	32.0 ± 10.5	49.0	<1.0 × 10^–15	93.9 ± 0.2	62.1 ± 1.3	31.8	<0.0001
	chr16: 18585845	76.7 ± 1.9	34.3 ± 4.1	42.4	<1.0 × 10^–15	92.0 ± 0.3	62.8 ± 1.4	29.3	<0.0001
	chr16: 18585857	85.0 ± 3.5	35.3 ± 8.1	49.7	<1.0 × 10^–15	94.4 ± 0.3	58.9 ± 2.1	35.5	<0.0001
(C)
		Methyl-Seq				Bisulfite-pyrosequencing
Gene	Location	Average % methylation (Cortex)	Average % methylation (Hippocampus)	Average % methylation difference	P-value	Average % methylation (Cortex)	Average % methylation (Hippocampus)	Average % methylation difference	P-value
Trim7	chr11: 48639597	3.3 ± 1.5	32.0 ± 1.2	−28.7	<1.0 × 10^–15	17.4 ± 1.6	19.1 ± 1.6	−1.8	N.S.
Cacbn4	chr2: 52427198	24.1 ± 2.4	54.7 ± 1.9	−30.4	<1.0 × 10^–15	51.3 ± 4.9	62.0 ± 2.3	−10.8	<0.0001
	chr2: 52427219	45.0 ± 2.6	88.3 ± 2.6	−43.3	<1.0 × 10^–15	84.8 ± 3.3	87.1 ± 2.4	−2.4	N.S.
	chr2: 52427253	42.3 ± 1.2	89.7 ± 2.2	−47.4	<1.0 × 10^–15	89 ± 1.6	92.6 ± 1.1	−3.6	N.S.
	chr2: 52427269	45.0 ± 1.0	94.3 ± 0.7	−49.3	<1.0 × 10^–15	86.4 ± 1.7	90.5 ± 0.8	−4.1	N.S.
	chr2: 52427357	42.0 ± 2.1	84.7 ± 4.3	−42.7	<1.0 × 10^–15	79.0 ± 2.4	85.5 ± 1.2	−6.5	0.0003
	chr2: 52427382	41.0 ± 1.7	83.0 ± 2.9	−42.0	<1.0 × 10^–15	78.4 ± 7.1	85.5 ± 5.8	−7.1	<0.0001
	chr2: 52427404	37.7 ± 2.8	76.3 ± 1.2	−38.6	<1.0 × 10^–15	79.0 ± 6.6	85.6 ± 2.3	−6.6	0.0002
Pot1a	chr6: 25760244	28.7 ± 1.5	62.0 ± 1.5	−33.3	<1.0 × 10^–15	43.6 ± 10.8	68.1 ± 11.4	−24.5	0.0006
	chr6: 25760252	25.7 ± 2.2	63.3 ± 2.3	−37.6	<1.0 × 10^–15	45.0 ± 9.9	68.4 ± 14.2	−23.4	0.0012
	chr6: 25760294	21.7 ± 0.9	51.0 ± 0.6	−29.3	<1.0 × 10^–15	45.0 ± 10.6	65.4 ± 13.2	−20.4	0.0056
	chr6: 25760303	23.3 ± 1.9	61.0 ± 2.6	−37.7	<1.0 × 10^–15	42.5 ± 11.3	64.4 ± 12.4	−21.9	0.0026
Prtg	chr9: 72704618	92.0 ± 1.5	63.0 ± 5.0	29.0	1.2 × 10^–10	97.4 ± 2.6	88.1 ± 3.6	9.3	0.0005
	chr9: 72704629	92.7 ± 1.2	59.7 ± 2.4	33.0	5.6 × 10^–14	95.0 ± 2.9	79.0 ± 4.6	16.0	<0.0001
	chr9: 72704638	94.0 ± 1.0	63.3 ± 1.9	30.7	1.1 × 10^–13	94.6 ± 1.6	83.0 ± 2.8	11.6	<0.0001
	chr9: 72704658	93.7 ± 2.4	56.0 ± 2.1	37.7	<1.0 × 10^–15	99.5 ± 0.7	82.0 ± 7.1	17.5	<0.0001
	chr9: 72704669	95.7 ± 1.9	52.0 ± 4.2	43.7	<1.0 × 10^–15	74.3 ± 1.3	66.1 ± 8.6	8.1	0.003
	chr9: 72704672	90.7 ± 0.9	60.3 ± 2.3	30.4	4.7 × 10^–11	74.8 ± 3.4	68.6 ± 7.3	6.1	0.049
	chr9: 72704722	90.3 ± 0.7	60.3 ± 3.0	30.0	6.5 × 10^–9	92.4 ± 1.5	77.9 ± 4.5	14.5	<0.0001
Apcdd1	chr18: 63083916	67.3 ± 3.8	22.0 ± 1.7	45.3	<1.0 × 10^–15	76.3 ± 4.3	41.6 ± 5.4	34.6	<0.0001
	chr18: 63083934	60.7 ± 2.3	23.0 ± 6.1	37.7	7.2 × 10^–11	67.4 ± 2.5	39.6 ± 5.4	27.8	<0.0001
	chr18: 63084105	76.0 ± 2.5	39.7 ± 1.9	36.3	<1.0 × 10^–15	77.6 ± 3.6	50.1 ± 3.8	27.5	<0.0001
Msx1	chr5: 38217514	4.0 ± 1.2	33.0 ± 2.3	−29.0	<1.0 × 10^–15	68.8 ± 25.3	86 ± 3.2	−17.3	N.S.
	chr5: 38217519	5.7 ± 0.9	38.3 ± 3.5	−32.6	<1.0 × 10^–15	65.3 ± 30.4	88.5 ± 2.5	−23.3	N.S.
	chr5: 38217538	7.0 ± 1.2	47.7 ± 2.3	−40.7	<1.0 × 10^–15	89.3 ± 18.4	96.8 ± 1.2	−7.5	N.S.
Gm13445	chr2: 35346878	76.3 ± 3.3	38.3 ± 2.0	38.0	1.6 × 10^–6	85.1 ± 1.7	38.0 ± 5.1	47.1	<0.0001
	chr2: 35346892	96.7 ± 1.8	54.0 ± 1.0	42.7	1.8 × 10^–12	94.5 ± 0.7	42.5 ± 6.0	52.0	<0.0001
	chr2: 35346894	92.7 ± 1.3	48.7 ± 0.9	44.0	7.5 × 10^–13	94.0 ± 0.9	53.9 ± 6.4	40.1	<0.0001
	chr2: 35346924	85.3 ± 5.6	36.7 ± 2.9	48.6	2.2 × 10^–12	88.1 ± 1.2.	54.0 ± 4.4	34.1	<0.0001
	chr2: 35346926	88.0 ± 0.6	53.0 ± 3.8	35.0	1.9 × 10^–8	92.0 ± 1.2	63.1 ± 3.1	28.9	<0.0001
Zfp825	chr13: 74671656	80.3 ± 4.2	43.0 ± 3.5	37.3	4.2 × 10^–8	83.5 ± 2.7	53.8 ± 2.5	29.8	<0.0001
	chr13: 74671661	85.0 ± 2.1	50.3 ± 0.3	34.7	4.0 × 10^–8	84.8 ± 2.2	58.1 ± 2.3	26.6	<0.0001
	chr13: 74671734	84.7 ± 2.0	47.0 ± 1.0	37.7	1.0 × 10^–9	86.9 ± 4.7	57.5 ± 2.8	29.4	<0.0001
	chr13: 74671739	88.3 ± 1.5	51.0 ± 2.9	37.3	9.3 × 10^–11	82.9 ± 2.8	54.3 ± 3.2	28.6	<0.0001
Cmtm1	chr8: 106833817	67.7 ± 6.2	30.3 ± 4.1	37.4	1.2 × 10^–6	94.5 ± 1.9	84 ± 3.3	10.5	<0.0001
	chr8: 106833863	74.3 ± 1.8	43.7 ± 2.8	30.6	5.1 × 10^–4	91.1 ± 3.1	77.1 ± 2.8	14.0	<0.0001
	chr8: 106833866	79.3 ± 3.2	37.7 ± 3.9	41.6	8.7 × 10^–8	95.0 ± 2.7	84.5 ± 5.0	10.5	<0.0001
	chr8: 106833874	81.0 ± 4.0	41.7 ± 2.8	39.3	4.3 × 10^–6	80.1 ± 1.7	65.6 ± 4.5	14.5	<0.0001
	chr8: 106833881	81.7 ± 4.7	44.0 ± 4.2	37.7	5.5 × 10^–6	77.3 ± 2.6	65.0 ± 4.6	12.3	<0.0001

Brain region-specific methylation differences

In the Methyl-Seq experiment, a total of 14,519 DMCs were identified in the comparison between cortex and hippocampus. Our validation work described above examined regions in 9 genes that showed differences in methylation between these 2 brain regions. Significant differences by bisulfite pyrosequencing were observed in 7 of these. More methylation was seen in cortex than hippocampus for CpGs in 5 genes: Prtg, Apcdd1, Gm13345, and Zfp825, and Cmtm1. Conversely, more methylation was seen in hippocampus than cortex for CpGs in 2 genes: Cacbn4 and Pot1a.

Discussion

Mice are commonly used to model different paradigms of disease. To this end, Methyl-Seq baits were developed to enable a targeted approach to studying patterns of the mouse methylome. To determine the validity of this assay, we investigated DNAm patterns between tissues that are functionally and structurally different from each other, with the expectation that such distinct tissues would maximize the observed methylation differences. Indeed, data from Methyl-Seq did reveal about 10 times more DMCs and a larger magnitude of difference between the brain and liver tissues compared to differences between the cortex and hippocampus, which are more similar. Strikingly, a disproportionate number of these DMCs reside in intronic and intergenic regions. Similar observations have been reported in previous genome-wide DNAm studies that investigated cell-type and tissue-specific DNAm differences.^4,18 Similarly, DMCs were also observed predominantly outside of CpG islands in CpG shores, shelves and open sea, as has been similarly observed in other genome-wide DNAm studies.^4,18,20 Together, these observations suggest that changes to DNAm patterns in introns, intergenic regions, and in regions outside of CpG islands may play an important role in tissue-specific expression of genes. Importantly, this study highlights that the current mouse Methyl-Seq assay can adequately capture regions where DNAm differences are present between tissues.

Although we did not investigate the function of these non-coding regions, many other studies have clearly demonstrated the involvement of non-coding regions in mediating tissue-specific gene expression. For example, an intron present in a gene ∼1 Mb upstream of the sonic hedgehog gene was shown to mediate limb bud expression of sonic hedgehog.²⁴ In an independent study, a conserved non-coding element ∼42 kb upstream of the galanin gene was observed to drive promoter expression in the paraventricular nucleus, arcuate nucleus and amygdala where the galanin peptide is expressed.²⁵ Our assay was designed to capture such genomic regions with putative regulatory function, and the finding that many of them contain tissue-specific DMCs suggests that DNAm may influence tissue-specific expression of their corresponding target genes. Indeed, the DMRs within which validated DMCs are located contain histone modifications and DNase I hypersensitivity sites, which are hallmarks of transcriptional regulatory features.²³ Furthermore, liver-specific enhancers, where regions are enriched for histone marks of active enhancers H3K4me1 and H3K27ac in the liver but depleted in the cortex, are highly methylated in the cortex. Likewise, cortex-specific enhancers are highly methylated in the liver. Interestingly, enhancers that are commonly shared between tissues also displayed tissue-specific methylation differences. This suggests that tissue-specific differences may not only be mediated by the suppression of tissue-specific enhancers, but may also involve varying methylation patterns within shared enhancers.

To demonstrate the utility of the new Methyl-Seq assay, we performed tests of gene ontology (GO) annotation enrichment on differentially methylated genes to determine whether they were enriched for terms that reflect the cellular identity of brain and liver tissues. This was indeed observed. Genes that are hypomethylated in the brain are distinctly enriched in neuronal components. In contrast, genes that are hypomethylated in the liver are enriched in cellular components prevalent in the liver, such as the microbody, or that play a role in the metabolic function of the liver. Similarly, pathway analyses also show that genes that are hypomethylated in the brain are enriched in signaling pathways involved in neuronal function. For example, the Wnt signaling and neurotrophin pathways have been shown to play a role in neurogenesis and arborization of neurites, which are important cellular processes that affect behavior and cognition.^26,27 Likewise, the Erb signaling pathway has been observed to mediate specification of adult neural stem cells into oligodendrocyte lineage.²⁸ In contrast, genes that are hypomethylated in the liver are enriched in pathways involved in various liver functions including protein metabolism, haematopoietic cell lineage and complement and coagulation cascades. The reliability of Methyl-Seq at detecting DMCs was demonstrated by its high correlation with the same DMCs detected by bisulfite-pyrosequencing in a separate cohort of mice. In addition, genes that were validated as differentially methylated also generally showed decreased gene expression in tissues where they were more highly methylated. There were exceptions, suggesting that, for some of these genes, other factors, such as histone modifications, may be involved in controlling their expression. Together, these observations underscore the reliability of the assay at detecting differentially methylated cytosines.

As a separate goal, we aimed to extend prior observations about brain-specific differences.^4,29 Previous genome-wide DNAm studies have observed a number of tissue-specific differentially methylated genes in the brain involved in cellular mechanisms and signaling pathways, such as neurotransmission, neurotrophins and Wnt signaling, that have been implicated in neuropsychiatric disorders.^4,29,30 In this study, we have identified genes with a role in neurodevelopment that have not previously been reported to be differentially methylated across brain regions. For example, the Adenomatous Polyposis Coli Downregulated 1 (Apcdd1) gene is observed to be hypomethylated in the mouse hippocampus relative to the mouse cortex. This gene is a negative regulator of the Wnt signaling pathway and has a migratory role in astrogliogenesis.³⁰ Similarly, the msh homeobox1 (Msx1) gene is hypomethylated in the cortex relative to the hippocampus. Although this gene is known to be involved in orofacial development,³¹ disruption of this gene has also been observed to impair brain development resulting in mice that have hydrocephalus at birth among other neuroanatomical defects.³² Due to their involvement in brain development, it may be of interest to explore whether methylation changes in human to orthologous regions of the DMRs of these genes may be involved in neurodevelopmental disorders.

The present study has several limitations. First, while our Methyl-Seq data showed good consistency with bisulfite-pyrosequencing validation results, they were not in perfect agreement. The magnitude of difference in DNAm varied between methods. This likely reflects variation in methylation detection between platforms. It is also worth noting that the Methyl-Seq assay is designed to detect DNAm levels on the negative strand, which is a default design parameter of Agilent's SureSelect target enrichment system. Although every attempt is made to validate DNAm changes on the negative strand by bisulfite-pyrosequencing, this is not always possible due to limitations in primer designs with sequences of low complexity. Thus, it is possible that discrepancies between the 2 platforms may be due to the effect of hemi-methylation. None of the genes reported are, however, currently known to be hemi-methylated. Second, while our Methyl-Seq assay covers about 4% of the mouse genome with annotated regulatory features or with regions where differential methylation is known to occur, there remain many more regions where relevant DNAm variation between tissues might exist. These will likely only be covered when whole genome sequencing for methylation status becomes affordable. Third, due to the increasing awareness of methylation outside of the CpG context, an external control should be spiked into the Methyl-Seq library preps to provide accurate detection of bisulfite conversion error rates. This is presently not performed using the standard Agilent Methyl-Seq protocol. In this study, bisulfite conversion error rates were computed using methylation levels detected from random cytosines in the non-CpG context within the sample. Although a slightly elevated bisulfite error rate was observed in brain tissue compared to liver tissue (1.5% vs 0.8%), which is in line with elevated methylation levels of the non-CpG context in the brain,³³ this had little impact on our ability to accurately detect tissue-specific DMCs and DMRs, as observed from our robust validation result.

In summary, we have demonstrated the robustness of this platform at detecting DNAm differences between varying mouse tissues. We have also shown the ability of this platform to identify novel DMCs and DMRs. The platform may be a valuable tool in future studies for examining DNAm changes in various mouse models of disease including those involving exposure to environmental stimuli such as stress, medications and substance abuse.

Material and Methods

Tissue extraction from mice

Adult C57/Bl6 mice (>12 weeks) were purchased from Jackson Laboratories. Mice were sacrificed followed by immediate extraction of the brain and liver. Cortices and hippocampi were microdissected from the brain and retained for the study. For Methyl-Seq experiments 3 biological replicates of these tissues were used. For subsequent bisulfite-pyrosequencing validation, 8 biological replicates of these tissues were used from a separate mouse cohort.

DNA extraction

DNA extraction was performed using the MasterPure^TM DNA extraction kit (Cat. No. MCD85201 Epicenter, an Illumina company, USA) as described in the manufacturer's protocol. Briefly, tissue was lysed in 600 μl of tissue and cell lysis solution with 2 μl of proteinase K (50 μg/μl). The reaction was briefly vortexed and incubated at 65°C for 30 min. The reaction was then incubated with 2 μl of RNase A (5 μg/μl) at 37°C for 30 min. After incubation at 4°C for 5 min, reaction was mixed with 350 μl of MPC protein precipitation reagent and centrifuged at 4°C for 10 min at 14,000 rpm. The supernatant was mixed with an equal volume of isopropanol and centrifuged at 4°C for 10 min at 14,000 rpm. The pellet was washed in 70% ethanol twice before air-drying for 10 min to remove ethanol residue prior to being dissolved in Tris-EDTA buffer.

Design of the custom RNA hybridization bait library

Recent studies have highlighted the importance of the non-coding genome for regulating gene expression.³⁴ We have designed baits to not only capture RefSeq genes, but to also include putative regulatory sequences, as DNA methylation patterns in both types of genomic regions may reveal important functionally relevant information. The commercially available Agilent SurePrint® oligo manufacturing and SureSelect® hybridization target enrichment technology uses biotinylated 120 base pair RNA baits to selectively hybridize with complimentary DNA sequences that have been randomly fragmented. Agilent offers a fully customizable bait design option through their eArray® online utility (https://earray.chem.agilent.com/earray/). Using this utility, we designed a custom hybridization array that targets 99 Mb of CpG islands, shores, shelves and regulatory elements in the mouse genome based on criteria described below. This is the most comprehensive murine methylome hybridization capture sequencing array currently reported.

Using the publicly available databases from the National Center for Biotechnology Information (NCBI-NIH) (http://www.ncbi.nlm.nih.gov), the University of California Santa Cruz (UCSC) Genome Browser (http://genome.ucsc.edu) and Ensembl (European Bioinformatics Institute/Wellcome Trust Sanger Institute) (http://www.ensembl.org), genomic coordinates for all known CpG islands, shores and shelves in the mouse genome (mm9) were obtained. Verification of each coordinate set was carried out against previously published reports of the predicted CpG islands in the mouse genome.^35-41 Additionally, we included tissue-specific differentially methylated regions (DMRs),²⁰ DNase I hypersensitive sites, and regulatory elements such as promoters, enhancers and transcription factor binding sites (Table 1). The combined total mouse methylome design included ∼201,788 targets and ∼110 Mb of bases covered by 297,773 baits. Due to some overlap between target regions, the final product was 99 Mb. Table 1 shows the composition of each of the categories included in the target design. Of the 297,773 baits that were designed, ∼75% of them successfully isolated targeted regions for sequencing. Of these, ∼51% were present in gene bodies covering 17,190 genes with an average of 80 CpGs per gene. This capture library can be obtained from Agilent (Cat. No. 931052). An annotated table of each bait probe can be found in Table S1. Table S2 contains a list of probe baits that did not generate usable data.

Methyl-Seq library preparation

Methyl-Seq was performed using Agilent SureSelect Methyl-Seq kit (Cat. No. G9651A) as described in the manufacturer's protocol (Methyl-Seq protocol, version A.1, April 2012). Briefly, 3 μg of DNA was used for the library prep. From this, >550 ng of adaptor ligated DNA was used for the hybridization capture and the final concentration of indexed library was 8–15 nM. The following changes to minimize loss of DNA in the process of enzymatic reactions. AMPure XP beads (Cat. No. A63880, Beckman Coulter, USA) were incubated with DNA reaction mix for 10 min at room temperature prior to pelleting by magnetization. AMPure XP beads were then washed twice with 80% ethanol and dried at 37°C for 5 min. DNA was dissolved in water for 10 min at 37°C. The amount of water used for dissolving is as described in the manufacturer's protocol. AMPure XP beads were retained in solution for adenylation and end-repair reaction. After these reactions, an equal volume of binding buffer (20% PEG 8000/2.5M NaCl) was added to the reaction mix (1:1 ratio) to enable the AMPure XP beads to rebind to DNA.⁴² The reaction was incubated at room temperature for 10 min. DNA was subsequently purified as described above. AMPure XP beads were not retained after this step. Concentration and sizes of DNA fragments after shearing, after adaptation ligation, and after indexing, were analyzed by Agilent Bioanalyzer 2100, and were observed to be within the range advised by the manufacturer's protocol. Greater than 550 ng of adaptor ligated DNA was used for the hybridization capture as described in the manufacturer's protocol. This protocol is described in detail in the supplementary material and methods. The Agilent protocol used in this study can be found here: http://www.medicine.uiowa.edu/psychiatry/SureSelect%20Methyl-Seq%20protocol%20version%20A.1.pdf

Bisulfite conversion

Bisulfite conversion was performed using EZ DNA Methylation-Gold^TM (Zymo Research, USA) as described in the manufacturer's instructions. DNA was eluted in 12 μl of elution buffer.

Next-generation sequencing

Sequencing was performed by 100 bp paired-end sequencing by Illumina HiSeq2000 as described in the manufacturer's protocol. Quality and concentration of library preparation were determined by bioanalyzer using a high sensitivity DNA kit (Cat. No. 50674626, Agilent, USA) and a KAPA real-time PCR assay (Cat. No. #KK4835, KAPA biosystem, USA). For HiSeq, 3 samples were multiplexed in a lane to a final concentration of 12 pmole for sequencing. This provided a cluster density between 600–700 K/mm². To overcome color imbalance inherent to low complexity in a bisulfite-converted genome, 1% of phiX genome was spiked into the reaction. In addition, a control lane containing regular genome that was not bisulfite converted was processed in the same flow cell as bisulfite converted libraries. Q30 scores of bases from HiSeq reactions were within the threshold recommended by the manufacturer.

Nested PCR for validation experiment

Validation of differentially methylated cytosines as detected by Methyl-Seq was performed by bisulfite pyrosequencing. To accomplish this we first designed nested PCR primer sets using Primer 3⁴³ and Pyromark assay design SW 2.0 (Cat. No. 9019077, Qiagen, USA). Due to the low complexity of bisulfite converted DNA, nested PCR was performed to promote specific amplification of targeted sequence. Primers were purchased from Integrated DNA Technology, USA. For inner PCR primer sets, one of the inner primers was biotinylated and HPLC purified as indicated in Table S3. Outer PCR amplification was performed with 1 μl of bisulfite converted DNA using either Taq polymerase with Thermopol® buffer (Cat. No. M0267S, New England Biolabs, USA) or the PyroMark PCR kit (Cat. No. 978703, Qiagen, USA). Inner PCR amplification was performed using 2 μl of outer PCR product. Thermal cycling conditions were performed as described in the manufacturers’ instructions. For thermopol PCR, annealing temperatures for outer and inner primer sets were 57°C and 53°C, respectively. For Pyromark PCR, both inner and primer sets were annealed at 56°C. The extension times for both PCR kits were 45 sec and 30 sec for inner and outer primer sets, respectively. The number of thermal cycles was 30 and 40 for inner and outer primer sets, respectively. The PCR products were electrophoresed on a 1% agarose gel to check for product specificity. Table S3 describes the primer sets used for the amplification.

Pyrosequencing

Pyrosequencing was performed using capillary dispensing tips on the PyroMark Q96 MD (Qiagen, USA) as described in the manufacturer's instructions. PyroMark CpG software (Cat. No. 9019067, Qiagen, USA) was used in this assay. Bisulfite conversion was shown to be efficient for all samples as the fluorescence signal by cytosine in a non-CpG context was ≤ 1% of the signal produced by thymine. Sequence specificity was also ascertained by introducing a base that is not in the target sequence. As an additional quality control, we investigated DNA methylation levels at the promoter of the glial fibrillary acidic protein (GFAP) gene in NeuN+ and NeuN- cells. NeuN+ cells showed higher methylation in the GFAP promoter than NeuN- cells (Fig. S3). This experiment demonstrated the reliability of the nested PCR bisulfite-pyrosequencing assay.

Data analysis

Bisulfite sequencing data from HiSeq were trimmed using Trim Galore! (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) to remove standard Illumina paired-end adaptors using the settings –paired -t -q 30.⁴⁴ Data were subsequently aligned to the mouse genome mm9 using Bismark and Bowtie2 with default settings and deduplicated using a Bismark perl script.⁴⁵ Mapping efficiency was ∼84%. Bisulfite error rate was estimated from the non-CpG context giving an average of ∼0.9% in liver and 1.5% in cortex and hippocampus. These estimates are comparable to the bisulfite conversion error rate of ∼1% previously reported in other studies.⁴⁶ A higher bisulfite conversion error rate from brain tissues is expected as methylation in the non-CpG context has been observed to be elevated in the brain.^47,48 Data were subsequently analyzed using the R package MethylKit.⁴⁹ Only cytosines with a minimal coverage of 10× and phred quality score of Q30 in the CpG context were analyzed. The 5-methylcytosine maps and lists of tissue-specific differentially methylated cytosines can be found at http://www.medicine.uiowa.edu/Psychiatry/Cortex_Hippocampus_Liver_DNA_methylation_map.7z and http://www.medicine.uiowa.edu/Psychiatry/Tissue_specific_DMCs.7z, respectively.

Analysis of differential methylation

Differentially methylated cytosines (DMCs) were considered as cytosines whose difference in methylation level between 2 tissues was > 27% (see below for rationale in selecting this threshold), with the q-value for the difference ≤ 0.05. The initial P-values of DMCs were calculated using logistic regression by comparing the number of methylated reads (Cs) and unmethylated reads (Ts) between the different tissues. The depth of coverage (Cs and Ts) per CpG was ≥ 10× with a median of ∼40–50×, values that were consistent between tissue comparisons. P-values were converted to q-values by MethylKit using the sliding linear model method to correct for multiple comparisons. Differentially methylated regions (DMRs) were determined as targeted regions with > 1 DMCs. For DMR analysis, we began with the default approach used in MethylKit, in which the number of Cs and Ts for all CpGs in a targeted region were combined, with Ns again reflecting the number of reads per CpG. Counts (Cs and Ts) were then divided by the total number of CpGs in the targeted region to produce an average number of Cs and Ts for the region. The median depth of coverage for the averaged region was ∼40–50×, values that were consistent between tissue comparisons. P-values were then calculated by comparing the number of Cs and Ts between tissues using logistic regression as described above and corrected for multiple testing by false discovery rate (FDR). DNAm differences identified by bisulfite pyrosequencing were analyzed by either Student's t-test or one-way ANOVA with Sidak's post-hoc multiple comparison test using Graphpad prism. For comparisons of absolute methylation levels between different tissues, the Mann-Whitney test with Bonferroni correction was employed.

Annotation

Differentially methylated cytosines (DMCs) were mapped to different CpG contexts (i.e., CpG islands, CpG shores, CpG shelves and CpG open sea), genomic contexts (i.e., promoters, exons, introns, and intergenic regions) and enhancers by overlapping DMCs to targeted capture baits mapped to the CpG contexts, genomic contexts or enhancers. Coordinates of CpG islands from the mm9 genome were extracted from the UCSC Genome Browser using Galaxy.^50-52 CpG shores were identified as 2 kb regions flanking CpG islands as previously defined.²⁰ Likewise, CpG shelves were considered as regions between 2–4 kb flanking CpG shores,⁵³ while CpG open sea was considered as regions beyond CpG shelves.¹⁸ The coordinates of promoters for the mm9 genome were extracted from the UCSC Genome Browser. Promoter coordinates were considered as 2 kb regions upstream of annotated transcription start sites of RefSeq genes with annotated 5′ untranslated regions. Baits that do not overlap promoter coordinates were mapped to exons, introns and intergenic regions using Peak Annotation and Visualization (PAVIS).¹⁹ Coordinates of cortex and liver enhancers were retrieved from the UCSC Genome Browser defined by peaks of histone 3 lysine 4 mono-methylation (H3K4me1) and histone 3 lysine 27 acetylation (H3K27ac) obtained from chromatin immunoprecipitation (ChIP)-seq of the cortex and liver of adult (8 weeks) C57/Bl6 mice. We note that the sex of the mice and collection procedures used were not available. As H3K4me1 and H3K27ac partially overlap promoters, only H3K4me1 and H3K27ac peaks which were >2 kb from the transcriptional start site were used. The expected number of differentially methylated CpGs in each of the CpG contexts, genomic categories, and enhancer regions was calculated as follows. First, Methyl-Seq baits were mapped to the different CpG contexts, genomic categories, and enhancer regions. Sequences for the baits were retrieved using Galaxy. The number of CpGs mapped to each CpG context/ genomic category was calculated and expressed as a ratio of the total number of CpGs in all CpG contexts, genomic categories, or enhancer regions. This ratio was then multiplied by the total number of DMCs observed to determine the expected number of DMCs for each CpG context, genomic category, and enhancer region.

Since most previous studies have focused on methylation in the CpG context, we chose to evaluate tissue-specific differentially methylated cytosines (DMCs) in the CpG context for this study. Additional data on CHG and CHH contexts (where H is A, T or C) was generated. Given the unique challenges in analyzing and interpreting those data, we chose to consider them in a subsequent report.

Gene Ontology and Pathway Analysis

Gene ontology enrichment and pathway analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was performed using the functional annotation tools of DAVID bioinformatics resources version 6.7.^54,55 Reported P-values have been corrected for multiple testing using the Benjamini-Hochberg method.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Supplemental Material

Supplemental data for this article can be accessed on the publisher's website.