Characterization of the DNA methylome and its interindividual variation in human peripheral blood monocytes

Subjects

The study was approved by the University of Missouri-Kansas City (MO, USA) Institutional Review Board, and each participant signed an informed consent document before entering the study. A total of 18 unrelated healthy postmenopausal Caucasian females, aged 50–56 years, were recruited through the Kansas City Osteoporosis Study, a repertoire of approximately 6000 subjects collected for genomic and epigenomic studies of complex diseases/traits. All subjects were living in Kansas City (MO, USA) and its surrounding areas, and were self-identified as being of European origin. An extensive set of exclusion criteria (Supplementary Material, see online at www.futuremedicine.com/doi/suppl/10.2217/EPI.13.18) was adopted for this study in order to exclude subjects with diseases/conditions that may potentially lead to altered DNA methylation and/or gene expression patterns in PBMs. The basic characteristics of the study subjects are summarized in Supplementary Table 1.

Sample preparation

Peripheral blood mononuclear cells (PBMCs) were isolated from 60 ml of freshly collected peripheral blood from each subject, using Lymphoprep™ (Axis-Shield, Oslo, Norway). PBMs were then isolated from the PBMCs with the Dynabeads^® Untouched™ Human Monocytes kit (Life Technologies, CA, USA) using a previously established and routinely performed protocol [19,23–25]. The kit depleted unwanted cells (i.e., T and B cells) from the PBMCs, leaving PBMs free of the surface-bound antibody and beads with minimum disturbance. The isolated PBMs were visually checked for purity and counted under a microscope. Genomic DNA and total RNA were extracted from the freshly isolated PBMs with the AllPrep^® DNA/RNA/Protein Mini Kit (Qiagen, CA, USA).

MeDIP-seq

MeDIP-seq was performed by Arraystar Inc. (MD, USA) using a protocol based on Down et al. with minor modifications [26]. In brief, genomic DNA was sonicated to approximately 100–500 bp with a Bioruptor^® sonicator (Diagenode, NJ, USA). A total of 1 µg of sonicated DNA was end repaired, A-tailed and ligated to single-end adapters following the standard Illumina (CA, USA) genomic DNA protocol. Methylated DNAs were enriched by immunoprecipitation with an anti-5-methyl-cytosine monoclonal antibody (MAb-081-100, Diagenode). Precipitated DNA fragments were purified using Qiagen MinElute^® columns and PCR amplified using the single-end Illumina PCR primers. The PCR products were separated on 2% agarose gel to select fragments in the approximately 200–300 bp size range. The completed libraries were quantified and checked for quality with Agilent (CA, USA) DNA Bioanalyzer 2100. In addition, an aliquot of each library was analyzed by real-time PCR for a specific methylated site of H19 locus to confirm the enrichment of the methylated region (Supplementary Table 2).

The libraries were denatured with 0.1 M NaOH to generate single-stranded DNA molecules, loaded onto Illumina flow cells, amplified in situ using a standard cluster generation kit v2 (Illumina), and subsequently sequenced on an Illumina Genome Analyzer IIx (with single-end, 36-bp reads). Raw sequencing image analysis and base calling were performed using Solexa (CA, USA) pipeline v1.6. The MeDIP-seq data has been uploaded onto the Gene Expression Omnibus under accession number GSE39604 [101].

Sequence alignment & data analyses

Sequencing reads that passed through the Solexa CHASTITY quality filter were aligned to the human genome reference sequence (GRCh37/hg19) using NOVOALIGN (v2.05; Novocraft, Selangor, Malaysia). DNA methylation profilies were inferred from the uniquely aligned reads by using the MEDIPS analysis package [27]. Briefly, each 36 bp read was extended to 250 bp along the plus or minus direction, depending on the strand information of the reads. The number of extended reads in every 50-bp bin throughout the genome was counted as the raw MeDIP-seq signals, which were then transformed into a reads per million format in order to assure that the methylation profiles derived from samples with different amounts of total sequencing reads were comparable. To correct for the influence of MeDIP enrichment by local CpG density [26,28], MeDIPS incorporated the coupling factor for each bin (representing the local CpG densities) into the transformed raw MeDIP-seq signals (reads per million). The corrected methylation signals were termed relative methylation scores. For any specified genomic region, the mean relative methylation score values of the 50-bp bins that fall into the specified region can be calculated and further corrected for the relative CpG density of the region to produce the absolute methylation score (AMS), which is an absolute methylation estimate allowing for comparisons of methylation profiles between genomic regions with different CpG densities.

To evaluate genome-wide DNA methylation profiles, AMS values were calculated for 1-kb overlapping sliding windows (step size: 0.5 kb) throughout the genome. In addition, the methylation profiles (AMS) were estimated for various classes of genomic features, such as CGIs and promoters. Specifically, promoters were defined as the −0.5 to +2-kb genomic regions around the transcription start sites (TSSs) of all RefSeq genes [29], and subdivided into three subclasses based on their CpG content [30]: high CpG promoters (HCPs; promoters containing a 500-bp interval with a CpG observed/expected ratio ≥0.6 and GC content ≥55%); low CpG promoters (LCPs; promoters containing no 500-bp interval with CpG observed/expected ratio ≥0.4); and intermediate CpG promoters (ICPs; promoters that are not HCPs or LCPs). CGIs were identified as regions of ≥200 bp in length with a GC content ≥50% and with a CpG observed/expected ratio of ≥0.6 [31]. The genomic locations of CGIs were obtained from the University of California (CA, USA), Santa Cruz Genome Browser [32]. CGIs were further grouped into three subclasses on the basis of their positions relative to the closest RefSeq genes, including 5′-end CGIs (from −1 to +0.3 kb of a RefSeq gene TSS), intragenic CGIs (from +0.3 kb of TSS to −0.3 kb of the gene transcription end site) and intergenic CGIs, which are CGIs that do not fall into the other two categories. CGI shores were defined as 2-kb regions on either side of the CGIs. Annotations of other genomic features were obtained from the University of California, Santa Cruz Genome Browser.

To determine the common methylation patterns across individual PBM methylomes, the AMS values for a specific class of genomic features (or genome-wide 1-kb windows) were transformed into (0–1000) intervals in each individual sample, with 0 representing the lowest methylation level and 1000 representing the highest methylation level within all examined regions [27]. Based on the bimodal distribution of transformed AMS (trans-AMS) values observed in our samples (data not shown), the methylation status of a specific region can be subsequently defined as follows: trans-AMS ≥600: highly methylated; 460< trans-AMS <600: partially methylated; and trans-AMS ≤460: lowly methylated.

MassARRAY^® EpiTYPER™ assay

Targeted methylation analysis was performed by using MassARRAY^® EpiTYPER™ assays (Sequenom Inc., CA, USA) designed for selected genomic regions. The EpiTYPER assay is a MALDI-TOF mass spectrometry-based method for quantitatively measuring the mean DNA methylation level within each short DNA fragment. The EpiTYPER assays were carried out by the Genomics Shared Resources Facility of Roswell Park Cancer Institute (NY, USA). Briefly, 1 µg of genomic DNA was modified by bisulfite treatment and regions of interest were amplified by PCR with specific primers, by which a T7-promoter tag was incorporated into the 5′-end of the reverse primers. Next, amplified fragments were transcribed in vitro into RNAs. The synthesized RNAs were treated with RNase A for base-specific cleavage. The cleavage products differ in molecular weight depending on the methylation state of the original genomic DNA, and the differences can be quantitatively analyzed by the MassARRAY analyzer system. DNA methylation standards (0, 20, 40, 60, 80 and 100% methylated) were used to control for PCR amplification bias.

Gene ontology analysis

The program GOEAST [33,102] was employed to identify significantly enriched gene ontology terms among list of genes with specific characteristics (e.g., genes with highly methylated promoters across the 18 subjects). The p-values were calculated by hypergeometric tests and adjusted for multiple comparisons by stringent Yekutieli (false discovery rate under dependency) adjustment [34].

Interindividual variation analysis

To identify genomic regions showing substantial interindividual variability in DNA methylation, we calculated median absolute deviation, a robust estimator of scale, for each 1-kb sliding window across the genome. Specifically, for each 1-kb window (i):

• Median absolute deviation_i = median (| AMS_ij − median [AMS_i]|)

where AMS_ij denotes methylation level (AMS value) of individual j at region i, and median (AMS_i) is the median of the AMS values at region i across the subjects.

Gene expression analysis

Genome-wide gene expression profiles were obtained for seven of the 18 subjects by using 1.0 Human Exon ST arrays (Affymetrix, CA, USA). Microarray data were processed using the robust multiarray analysis algorithm implemented in the R package of ‘oligo’ [35]. Gene level expression estimates were derived using the ‘core’ metaprobe list annotation release 21.

Correlation of the methylation profile with gene expression level

To examine the role of promoter methylation in influencing transcription levels across different genes, we first calculated the average promoter methylation levels and average transcription levels for each gene, and then compared the methylation levels between genes of different transcription levels. In addition to comparison between individual genes, we also divided the assessed genes into quintiles based on the ranking of the mean expression levels of individual genes and compared between the gene groups for the DNA methylation profiles in a ±1.5-kb region surrounding the TSSs. Similar comparisons between gene groups of different expression levels were also carried out separately in genes associated with HCPs, ICPs and LCPs.

To evaluate the correlation between promoter methylation and gene expression levels across individuals, we computed the Pearson’s correlation coefficient between promoter methylation levels and gene expression levels for each gene across all individuals. We defined a conservative threshold (Pearson’s r ≤ −0.678) for determining significant negative correlation. This threshold was based on the fifth percentile of the empirical distribution of the correlation coefficient under the null hypothesis of no association between DNA methylation and gene expression levels across individuals, which was obtained from the permutations of the data sets from the seven subjects with both DNA methylation and gene expression data collected.

Data generation & quality assessment

Using MeDIP-seq, we generated a total of approximately 553 million 36-bp single-end sequence reads in the 18 PBM samples. Of these raw reads, approximately 283 million reads passed the Solexa CHASTITY quality filter and were uniquely aligned to the human genome (GRCh37/hg19), resulting in an average of approximately 16 million uniquely aligned high-quality reads per sample (Supplementary Table 1). Based on the high-quality mapped reads, we estimated coverage at genome-wide 50-bp bins using the MEDIPS package [27]. For each individual sample, on average 64% (range: 52–67%) of the 28.2 million CpG sites in the human genome were covered at least once (Supplementary Table 3). In addition, saturation analysis indicated that we had sufficient reads to obtain reproducible genome-wide methylation profiles (Supplementary Figure 1).

In order to confirm the validity of our data and the reproducibility of the analyses, we first checked the genome-wide DNA methylation patterns for several previously reported features. In agreement with previous reports [36], we observed that subtelomeric regions were frequently found to be highly methylated and genomic regions associated with lightly stained Giemsa bands often showed higher methylation levels than darkly stained regions (Supplementary Figure 2). We examined methylation signals at specific genomic regions, such as the region harboring the GAPDH gene. GAPDH is a ubiquitously expressed housekeeping gene, and thus the promoter region of this gene was expected to be unmethylated. Indeed, we consistently observed high levels of GAPDH expression (ranked in the top 1% most highly expressed genes among the 20,384 genes assayed on the Affymetrix exon arrays) and a lack of methylation at the GAPDH promoter in all samples (Supplementary Figure 3), whereas other regions surrounding the gene showed variable methylation status. Similar results were also observed at several other housekeeping gene regions (data not shown). Additionally, we further validated our MeDIP-seq data by confirming the DNA methylation status at a number of selected genomic regions (five promoters and three CGIs) using MassARRAY EpiTYPER assays, which demonstrated a high concordance rate (Pearson’s r = 0.74–0.92) between the two approaches (Supplementary Table 4). Taken together, we have successfully generated comprehensive and reproducible genome-wide DNA methylation profiles in our PBM samples.

Characteristics of the PBM methylome

Based on the MeDIP-seq data, we characterized DNA methylation patterns at 20 distinct genomic features (Figure 1A). The overall DNA methylation patterns in PBMs were similar to previously reported patterns observed in other cell types [36–38]. For instance, promoters tend to be lowly methylated, whereas gene bodies, particularly protein coding regions, were more frequently found to be highly methylated (Figure 1A). Interestingly, we found that approximately 24% of CGIs were highly methylated. This number may at first seem to be higher than expected since it has been widely accepted that CGIs are largely unmethylated [39]. However, it has become increasingly evident that a significant fraction of CGIs were methylated across a variety of different tissues/cells [36,40–42].

HCPs and LCPs have been shown to exhibit distinct DNA methylation patterns in other human somatic and germline cells [38,43,44]. Accordingly, we performed further DNA methylation analyses for promoters by subdividing the 25,814 promoters into three classes, namely HCP (n = 12,416), ICP (n = 5510) and LCP (n = 7888), based on their CpG content (see ‘Materials & methods’ section). Consistent with previous findings in other cell types [38,43,44], we observed distinct methylation patterns between HCPs and LCPs (Figure 1B). As such, on average, 84% of HCPs (range: 60–97% in individual samples) were lowly methylated and less than 0.5% of HCPs (range: 0–2.3% in individual samples) were highly methylated. By contrast, on average, 28% of LCPs (range: 9–51%) were highly methylated and an additional 55% of LCPs (range: 39–68%) showed partial methylation. Consequently, highly methylated promoters were almost exclusively (average: 95%; range: 91–98%) detected in LCPs, whereas lowly methylated promoters were predominantly (average: 70%; range: 66–73%) observed in HCPs (Figure 1B).

Similarly, we took a closer examination of the CGIs by subdividing the CGIs into three groups, specifically, 5′-end CGIs (n = 9426), intragenic CGIs (n = 7078) and intergenic CGIs (n = 11,033), on the basis of their positions relative to the nearest RefSeq genes (see ‘Materials & methods’ section). CGIs located near the 5′-end of RefSeq genes were much less frequently found to be highly methylated than CGIs located in gene bodies and intergenic regions. For instance, on average, only 7% of 5′-end CGIs were found to be highly methylated, in contrast to 32–34% of intra- and inter-genic CGIs. Accordingly, highly methylated CGIs were most often to be found at the intra- and inter-genic regions (average 89%; range 83–96%) than at the 5′-end regulatory regions (average 11%; range 4 –17% ) (Figure 1C).

Next, we investigated the highly/lowly methylated regions (1-kb sliding windows) that were common to all 18 PBM methylomes (Figure 2A; profiles are saved as bedGraph files, which can be downloaded from the Supplementary Material and visualized on the University of California, Santa Cruz genome browser). Regions that were constitutively highly methylated (n = 128,982) or lowly methylated (n = 246,136) in PBMs were distributed throughout the genome. In general, highly methylated regions have lower CpG density (CpG observed/expected ratios) than lowly methylated regions (Figure 2B; p < 2.2 × 10⁻¹⁶). This is more evident at some long-stretch segments enriched for highly methylated or lowly methylated windows (Figure 2C–D). For instance, a 350-kb segment on chromosome 9q33.1 was found to be significantly enriched for highly methylated windows and a close examination of the sequence in this area revealed that this region has relatively low CpG density and is largely depleted for CGIs (Figure 2C). By contrast, a 500-kb CpG-rich region on chromosome 7q22.1 was largely enriched for constitutively lowly methylated windows (Figure 2D). Out of the 25,814 promoters, we identified 276 promoters that were consistently highly methylated in all 18 subjects, including three ICPs and 273 LCPs (Supplementary Table 5). Gene ontology analysis revealed that genes associated with these constitutively highly methylated promoters were significantly enriched in a number of biological processes (Supplementary Table 6), such as sensory perception (p = 5.27 × 10⁻¹⁰) and keratinization (p = 1.21 × 10⁻⁷). As many of these biological processes are not related to PBM function, constitutive silencing of these genes with highly methylated promoters may not be unexpected.

Interindividual variation in the PBM methylome

Previous studies suggested that there is substantial interindividual variation in DNA methylation and this variability may confer important functional implications for gene regulation and contribute to phenotypic variation in humans. Therefore, in addition to the qualitative characterization of common PBM methylomes, we also assessed the interindividual variability of DNA methylation in normal PBMs in a quantitative manner. Specifically, we measured the interindividual variability with median absolute deviation for each 1-kb sliding window throughout the genome and selected the top 0.5 and 1% most-variable windows for further characterization. These highly variable regions were not evenly distributed along the chromosomes (Figure 3A). Chromosomes 13 and 21 had a significantly higher density of variable regions than predicted from their length, whereas chromosomes 1 and 3 had the lowest density of variable regions across the genome. Furthermore, by annotation of the variable 1-kb windows to their closest genomic features, we found that relative to their sizes in the reference genome, transposable elements, such as satellite and long terminal repeats, were enriched in these highly variable regions; whereas genic components, such as promoters and exons, tend to be under-represented in these highly variable regions (Figure 3B).

Correlation between DNA methylation & gene expression

The modification of DNA by methylation is an important epigenetic mechanism that affects the spatial and temporal regulation of gene expression. Among the 18 samples used in this study, seven have been previously tested for genome-wide gene expression profile data in PBMs. Subsequently, we investigated the relationship between DNA methylation and gene expression for the 20,684 unique genes that were assessed on the 1.0 Human Exon ST arrays.

First, we compared the promoter DNA methylation and gene expression levels across individual genes (Figure 4A). On the one hand, genes with highly methylated promoters exhibited low expression levels and highly expressed genes showed low promoter methylation, supporting the notion that promoter methylation is correlated with gene silencing. On the other hand, low promoter methylation is not necessarily translated into high gene expression, as many silenced/lowly expressed genes also contained lowly methylated promoters, indicating that various other mechanisms may regulate gene expression independent of the effects of promoter methylation [44,45]. The negative correlation between DNA methylation and gene expression levels was also evident when we split the 20,684 assessed genes into five groups based on the ranking of mean expression levels of individual genes and examined the DNA methylation levels in a 1.5 kb region up- and down-stream of the TSSs (Figure 4B). As such, the top 20% most highly expressed genes exhibited the lowest methylation levels in the TSS-surrounding region and groups of genes with increasingly lower expression levels showed increasingly higher methylation levels. Additionally, we observed a gradual decrease in DNA methylation levels when approaching the TSS from both upstream and downstream, with the bottom of the methylation ‘valley’ at around ±450 bp from the TSS (Figure 4B). A similar pattern of DNA methylation around the TSS has been observed in human embryonic stem cells [46], lymphocytes [4] and PBMCs [37]. Since previous studies suggested that the relationship between promoter methylation and gene expression is related to the CpG density of the promoter [44], we further examined the methylation patterns at the TSS-surrounding region, and the correlation between DNA methylation and gene expression in genes associated with HCPs, ICPs and LCPs, respectively (Figure 4C). Similar to what we observed for the combined data, we found that DNA methylation around the TSS region was negatively correlated with gene expression within each gene group associated with a distinct promoter class, even for genes associated with LCPs. However, the DNA methylation ‘valley’ centered at the TSS was only observed for genes associated with HCPs and ICPs, whereas the methylation levels were relatively constant along the TSS-surrounding regions for genes associated with LCPs (Figure 4C). Subsequently, we explored the pattern of CpG density in the TSS-surrounding region in order to determine whether differential CpG density along the region can explain the observed differences in their DNA methylation patterns between genes associated with HCPs/ICPs and LCPs. Indeed, we observed distinct CpG density profiles surrounding TSSs when comparing genes associated with HCPs/ICPs to those associated with LCPs (Figure 4D). In the TSS-surrounding regions of genes associated with HCPs and ICPs, we observed a gradual increase in CpG density when approaching the TSS, with a peak center at the TSS. By contrast, CpG density was relatively constant throughout the TSS-surrounding region for genes associated with LCPs.

To explore how differential promoter methylation was correlated with variation in gene expression across individuals, we calculated the Pearson’s correlation coefficient between methylation and gene expression levels for each gene across all individuals. We observed a modest but significant excess (p < 2.2 × 10⁻¹⁶) of genes showing negative correlation between promoter methylation and transcription levels across individuals, particularly for genes associated with HCPs (Figure 5). We selected genes with correlation coefficients of r ≤ −0.678 (equivalent to the fifth percentile of the empirical distribution of the correlation coefficient) as significant genes. Gene ontology analysis for these significant genes revealed an enrichment of biological processes that are closely related to PBM function (Table 1 & Supplementary Table 7), suggesting that the interindividual variation in expression levels for genes involved in these functions may be under strong influence from individual differences in DNA methylation levels.