Intraindividual Variation and Short-term Temporal Trend in DNA Methylation of Human Blood

This article is featured in Highlights of This Issue, p. 481

Epigenetic modifications such as DNA methylation of cytosine residues at CpG dinucleotides have been extensively examined for their potential association with disease and aging in humans. DNA methylation marks are stable indicators of tissue lineage; thus, different tissues have fundamentally different methylation patterns. DNA methylation profile is established at neonatal stage and undergoes rapid change at the early age and further modifications throughout lifespan. Changes of methylation patterns have been seen in normal human aging process, and aberrant methylation changes in tissues such as blood have been linked to a number of diseases, including cancer (1–5). A link has also been established between age-related and cancer-related methylation changes. Important cancer-related genes become hypermethylated during aging, including key developmental genes and those encoding the estrogen receptor, insulin growth factor, and E-cadherin (1, 6). Furthermore, the acquisition of methylation changes during aging could underlie the development of age-related pathologic conditions (7).

Research into the role of DNA methylation in disease is affected by variability in measured methylation levels, which can be attributed to a number of sources: differences among study participants due to demographic, environmental and genetic factors; random and systematic variation between samples taken from the same individual, which may result from differences in tissue sample composition or timing of sample collection; technical variability due to measurement error and limited precision of analytic tools. The relative proportion of between- and within-person variability in the total variability is typically described by the intraclass correlation coefficient (ICC). The ICC varies between 0 and 1, with high (low) proportion of within-person variation represented by values close to 0 (respectively,1). Low values of ICC signify lower data reliability and may affect validity of analysis results. A number of studies (8–11) have analyzed within-person variability of different blood biomarkers; however, we are unaware of any reports on that in DNA methylation, although the importance of this aspect for epigenome-wide association studies has been recognized (12).

In the present study of leukocyte DNA methylation among 24 healthy women, we examined CpG-site–specific intraindividual variation in DNA methylation and analyzed short-term temporal changes in DNA methylation using Illumina HumanMethylation450 BeadChip (13). This array has been reported to have low technical variability and high reproducibility (13), although a number of issues have been pointed out (14). The BeadChip contains 485,577 loci that include 482,421 (99.35%) CpG dinucleotides, 3,091 (0.64%) CNG targets, and 65 (0.01%) SNP sites. The array covers 96% of known CpG islands (CGI) and 99% of NCBI Reference Sequence (RefSeq) genes, with an average of 17 CpG sites per gene distributed across the upstream of the transcription start sites (TSS)1500, TSS200, 5′ untranslated regions (UTR), first exon, gene body, and the 3′UTR. Because the function of DNA methylation seems to vary by functional location and by proximity to CGIs (15), we also sought to describe the distribution of high–low variability CpG sites across functional locations in the genome, and identify whether such changes are more likely to occur in specific types of genomic regions. To assess a potential effect of X chromosome inactivation on methylation variability, separate analyses were conducted for X chromosome CpG sites.

Study participants were drawn from the Shanghai Physical Activity (SPA) study subcohort of the Shanghai Women's Health Study (SWHS). SWHS, a prospective cohort study of 74,943 women ages 40 to 70 years who were recruited from seven communities in Shanghai, China, between 1997 and 2000 (16). The SPA subcohort of SWHS comprises 300 women from two communities (17), who provided blood samples at study entry and after 9 months. Twenty-four women from the SPA subcohort were randomly selected for the present analysis.

For each consented participant, a 10-mL blood sample was drawn into an EDTA vacutainer tube and kept in a portable Styrofoam box with ice packs at 0°C to 4°C. Within 6 hours of collection, samples were processed and separated into 2-mL aliquots of plasma (4), buffy coat (2), and red blood cells (2). Samples were then stored at −80°C until DNA extraction.

DNA was extracted from buffy coats using a QIAamp DNA Blood Mini Kit (Qiagen Inc.). The quality and quantity of extracted DNA was examined by NanoDrop-2000 spectrometer (Thermo Scientific). Bisulfite conversion of 500 ng of DNA was performed on each sample according to the manufacturer's recommendations for the HumanMethylation450 BeadChip using the EZ DNA Methylation Kit (Zymo Research). The treatment protocol included 16 cycles of denaturing at 95°C for 30 seconds, incubation at 50°C for 60 minutes, and holding at 4°C for a minimum of 2 hours.

Bisulfite-converted DNA (4 μL) was hybridized onto and analyzed using HumanMethylation450 BeadChip (Illumina). The hybridization protocol consisted of a whole-genome amplification step followed by enzymatic endpoint fragmentation, precipitation, and resuspension. Resuspended samples were hybridized onto the BeadChip for 16 hours at 48°C. After hybridization, unhybridized and nonspecifically hybridized DNA was washed away, followed by single-nucleotide extension using the hybridized bisulfite-treated DNA as a template. The Illumina iScan SQ scanner was used to create images of the single arrays. Image intensities were extracted using GenomeStudio (v.2011.1) Methylation module (v.1.9.0) software.

Data normalization was performed in GenomeStudio using “Background Subtraction” and “Normalization to Internal Controls” methods and has been described elsewhere (18). Briefly, background subtraction values were derived from built-in negative control bead signals and subtracted from probe intensities. Normalization was performed using internal control probe pairs designed to target the same region within housekeeping genes. The methylation score for each CpG site was represented by a β value, calculated according to normalized probe fluorescence intensity ratios between methylated and unmethylated signals and varying between 0 (fully unmethylated) and 1 (fully methylated). Data quality control (QC) analyses performed on β values included principal component analysis, to assess potential batch effects, and sample histograms for signal distributions (data not shown). Probes with detection P value above 0.05 were excluded.

Statistical analyses were performed using methylation levels represented by M values, computed as logit of β values: M = log₂(β/(1−β)). Because the distribution of M values is closer to normality, they are widely used as a measure of DNA methylation in association studies (19). For every CpG locus, we estimated variance components that correspond to within-person |$(\hat \sigma _w^{2})$| and between-person |$({\hat\sigma}_{b}^{2})$| variation using a mixed model with study participants as random effects:

ICCs were then computed as |${\rm ICC} = {\hat\sigma}_{b}^{2}/\left({\hat\sigma}_{b}^{2} + {\hat\sigma}_{w}^{2}\right)$|⁠. To summarize patterns of variability by functional genomic location, by type of probe and for known differentially methylated regions (DMR), we computed the percentage of CpG loci with high (>0.8), midrange (0.5–0.8), and low (<0.5) ICCs for each of these groups. The estimated between- and within-person variance components were converted to the β-scale as follows:

where |$\hat\mu$| is the mean M value across all study subjects. For 14,248 CpG loci with very low mean methylation, ICCs could not be estimated and were set to 0, as the worst possible case. The short-term temporal trend in DNA methylation across all CpG loci was examined using a mixed model with time since study entry (years) as a fixed effect, and study participants as random effects:

To assess the influence of quantile normalization on ICC and trend, all models were fit with untransformed and quantile normalized M values. All analyses were adjusted for cell-type composition using the cell mixture deconvolution method of Houseman and colleagues (20), which establishes DNA methylation signature panels for each cell type and uses constrained optimization to map methylation profiles of interest onto the signature panel and to predict proportions of cell types in a blood sample.

A number of recent reports define meaningful difference (14) in methylation as Δβ value (a measure of average difference in methylation level at a particular CpG site between study subjects) being above a certain threshold. A threshold of 0.2 on the β-scale is commonly used. To assess the impact of intraindividual variation on such high-variance CpG loci, we summarized ICC and trend by genomic location and probe type for high-variance (Δβ ≥ 0.2), mid-variance (0.1 ≤ Δβ < 0.2), and low-variance loci (Δβ < 0.1). We estimated the Δβ value for every CpG site as twice the between-person SE: |$\Delta \beta = 2\hat \sigma_{b}$|⁠. All analyses were conducted using R 3.0.3 software and Bioconductor package nlme 3.1–117.

The median age at study entry among the 24 study participants was 54.5 years (range, 46.5–68.8 years). Across the 483,880 CpG loci that passed data QC, the ICCs ranged from 0 to 0.999, with a median of 0.37. For the majority of loci, the within-person variance component tended to be smaller than and uncorrelated with the between-person component (Fig. 1). The ranges of |${\hat\sigma}_{b,\beta}$| and |$\hat \sigma_{w,\beta}$| were 0 to 0.42 and 0 to 0.36, respectively.

Table 1 lists proportions of low-, mid-, and high-ICC loci by genomic location, probe type, and locus variability for nonsex chromosomes. The median ICC across all CpG loci was 0.36 [interquartile range (IQR), 0.13–0.63]. More than 64% loci had low ICC, whereas 23% had mid-range ICC, and 13% had high ICC. Among the low-variability loci, the proportions of low- and mid-range ICC were comparable, whereas only 8.1% loci exhibited high ICC. On the other hand, more than 90% of moderate-variability loci (0.1 ≤ Δβ < 0.2) and more than 99% of high-variability loci (Δβ ≥ 0.2) had high ICC, whereas <1% loci in these two groups had low ICC. This latter result was observed across all genomic locations.

Among the CpG loci in known DMR (13), ICCs tended to be higher than across the entire genome, with median 0.49 (IQR, 0.24–0.72), and 51.5%, 31.8%, and 16.7% loci exhibiting low-, mid-range and high ICCs, respectively. Across loci from CGIs, shores and shelves, the median ICC ranged from 0.34 to 0.43, with the proportion of high-ICC loci highest in CpG shores (overall and among low-variability sites) and CGIs (moderate-variability sites). The ICC profiles across different functional locations were similar, with median ICC between 0.29 and 0.36, and the percentage of high-ICC loci 7.8% to 12.4% overall, 5.2% to 8.2% among low-variability sites, and 92.0% to 96.2% among moderate-variability sites. There was also little difference in ICC profiles by probe type. CpG loci with methylation measured by Infinium I probes, compared with Infinium II probes, had lower proportion of high-ICC loci overall (11.3% vs. 13.5%), among low-variability sites (6.6% vs. 8.7%), and among moderate-variability sites (92.7% vs. 94.7%).

The distribution of ICCs on the X chromosome (Table 2) differed substantially from that on nonsex chromosomes, with the median ICC of 0.80 (IQR, 0.60–0.89) overall and nearly half of the X chromosome CpG loci having high ICC. The main difference from nonsex chromosomes in the distribution of ICCs was found among low-variability sites, where 27.8% to 51.1% loci had high ICC. Of all genomic locations, lineage-defining DMRs appeared the most stable, with 61.6% loci having high ICC. CGIs and shores and functional locations on the promoter side of gene coding regions contained more high-ICC loci than CGI shelves, gene body, and 3′UTR locations. The proportion of high-ICC loci measured by Infinium II probes was somewhat higher than that among Infinium I probes. ICC estimates without cell-type composition adjustment were very similar (Supplementary Tables S1 and S2).

A summary of temporal trend estimates by genomic location, probe type, and locus variability is presented in Table 3. Less than 4% CpG loci exhibited a statistically significant trend overall and for most functional locations, with the exception of 5.3% loci with trend for 3′UTR on the X chromosome. There was little difference in the proportion of CpG sites with temporal trend between low-, mid-, and high-variability locus groups (data not shown). The proportion of sites with negative trend was somewhat higher in most functional locations. The proportion of loci with positive trend among Infinium II sites on the X chromosome was about 1.5 times that for Infinium I sites. Without adjustment for cell-type composition, the proportion of loci with significant negative trend was several times higher across functional locations, 8% to 17.5% on nonsex chromosomes and 5.2% to 9.4% on the X chromosome (Supplementary Table S3). We estimated the effect of quantile normalization transformation on ICC profiles and trend estimates of CpG loci. Across all genomic locations and probe types, the application of quantile normalization to DNA methylation data had no noticeable effect on the proportion of high-ICC loci, but lowered the proportion of mid-ICC loci and, correspondingly, increased the proportion of low-ICC loci (Fig. 2A). This effect was observed in both nonsex chromosomes and the X chromosome (data not shown). In temporal trend analysis, quantile normalization had the effect of equalizing the proportions with negative and positive trend, and an overall increase in the proportion of loci with trend (Fig. 2B).

Figure 3 presents a comparison of temporal trend estimates from untransformed and quantile-normalized models, with and without cell-type adjustment. The application of quantile normalization alone resulted in the reduction of the temporal effect magnitude, consistent across most loci (Fig. 3A). Cell-type adjustment alone left temporal effects nearly unaffected (Fig. 3B). Compared with cell-type adjusted model, quantile normalization resulted in inconsistent changes in the temporal effect magnitude, with larger changes for some loci and smaller changes for others (Fig. 3C and D). Neither cell-type adjustment nor quantile normalization, alone or combined, changed the direction of the observed temporal effects.

To our knowledge, this is the first study to look at within- and between-person variability profiles of CpG sites from the HumanMethylation450 BeadChip. Our results show that on average, ICC < 0.5 is common across the genome, and that ICC profiles are similar across all functional locations and probe types. This implies that genomic location at or near CGIs or functional regions has little effect on inter- and intraindividual variability of a CpG locus. We also found that within-person variation tended not to exceed a certain threshold (0.36 for |${\hat\sigma}_{w}$|⁠). Thus, larger between-person variation generally implies a higher ICC. In particular, among CpG loci with between-person variation above the threshold of meaningfulness (Δβ ≥ 0.2), more than 99% loci had high ICC. Therefore, for differentially methylated loci with sizeable differences in methylation across subjects, that is, loci of interest to most association studies, one measurement may be sufficient for a reliable estimate of methylation level.

We note that CpG loci within DMRs had better ICC profiles: even among those with low between-person variation, the percentage of high-ICC loci was higher than in other genomic locations. This observation is expected, as methylation of tissue-specific lineage-defining regions is stable in blood DNA.

The low ICC values may result in substantial attenuation (reduction in magnitude compared with the true value) in estimated parameters, so repeated measurements would be required per study participant to keep attenuation within some acceptable level. The average ICC across all CpG sites in our study was under 0.5, with the largest proportion of low- and mid-range ICC values observed among low-variability CpG loci. More than one measurement would be needed for these loci to adequately assess their methylation and to limit attenuation in the estimates; however, because their variability is well below the threshold of meaningful difference, they will likely not be primary targets of an association study.

It should be noted that the low-variability CpG locus group would include most of unmethylated loci with mean methylation close to 0, as such averages can only be attained with low variation. Because of the nature of the logit transformation, this translates to negative M values of very large magnitude, whereby even small differences in the β values may translate to M values that are wide apart. As a result, any estimates for such loci are very unstable and should be treated cautiously. For these loci, it may be preferable to use β values rather than M values, with appropriate adjustments in analytic methods. For example, beta regression techniques (21) could be used for comparison of DNA methylation between two groups.

We have found that the percentage of high-ICC loci is much larger on the X chromosome, implying more stability in methylation patterns. Prior studies have shown that X chromosome inactivation is accompanied by methylation increase at CGIs and at the promoters of genes silenced by X chromosome inactivation (22). Furthermore, methylated promoter CGIs are usually associated with genes in a stable long-term silenced state (15). Although we were unable to distinguish between active and inactive X chromosomes, our observation of more stable methylation patterns at CGIs and near the start of gene coding regions on the X chromosome appears consistent with the effect of gene silencing in X chromosome inactivation.

In our study, we have found that most CpG sites did not exhibit temporal trend that is sufficiently strong to be detected in a small window of <1 year. Among a rather small percentage of sites that did exhibit a trend, negative trend was somewhat more common than positive trend. This finding is in agreement with several prior studies of age-related methylation changes. Using HumanMethylation450 platform on 421 individuals ages 14 to 94 years, Johansson and colleagues (23) reported 29% of CpG sites significantly associated with age, of which 60.5% exhibited decrease, and 39.5% increase in methylation. Heyn and colleagues (24) showed that most of the genome undergoes age-associated hypomethylation. At the same time, a number of studies found that high CpG density promoters of key developmental genes tend to exhibit age-associated hypermethylation, during both early and late stages of life (6, 25, 26). In two longitudinal studies of newborns followed for 1.5 to 5 years, Martino and colleagues (26, 27) found clear distinction in methylation profile between samples collected at birth and at subsequent clinic visits. They also observed increase in methylation across all classes of annotated genomic regions, with intergenic regions most likely to undergo such changes. A unique aspect of our study is its longitudinal nature combined with a focus on short-term temporal changes. We have been able to show that in a number of CpG loci, temporal changes in methylation that occur in mature adulthood are detectable over a 1-year period of time.

In addressing intraindividual variability in DNA methylation, adjustment for known sources of such variability could improve the ICC. Temporal change may be one such source of within-person variability for CpG loci that exhibit methylation changes over time. Removing a temporal trend may potentially reduce intraindividual variation and improve statistical power. However, because age-related methylation change may be an important contributor to carcinogenesis or other pathogenic process, temporal trend removal may sometimes obscure a true association between DNA methylation and disease risk.

Different cell mixture composition between samples taken from the same individual may also confound measured DNA methylation and contribute to the observed within-person variability. Reinius and colleagues (28) established that CpG methylation differs between cell types, such as mononuclear cells, granulocytes, natural killer cells, B cells, and T cells. Jacoby and colleagues (29) analyzed methylation of 58 CpG sites and observed differences in interindividual variability across blood cell types. We adjusted our analyses for cell-type composition using the method of Houseman and colleagues (20). Koestler and colleagues (30) further tested this algorithm on data from the Gene Expression Omnibus database, reporting moderate to high agreement between predicted cell-type composition and that from complete blood cell counts. In our analyses, there was little difference in ICCs before and after cell-type adjustment. It has also been suggested that changes in cell composition of human blood across a person's lifespan may largely explain age-associated methylation change (1, 31). In our study, we found that a significant temporal trend disappears after cell-type adjustment in the majority of CpG sites with such trend (Table 3, Supplementary Table S3). We also observed an age-related decrease in the proportion of CD8⁺ T cells (3.5%/year; P = 0.015) and B cells (1.2%/year; P = 0.056) and a corresponding increase (4.7%/year; P = 0.044) in the proportion of granulocytes (data not shown). Thus, our results support the hypothesis that cell composition change largely accounts for temporal changes in DNA methylation. Although two of the aforementioned studies of methylation change did not account for cell-type composition of whole blood (23, 25), others either considered cell-type composition (27) or focused on specific cell types (6, 24). Thus, associations reported by this latter group of studies are independent of age-related changes in cell composition.

Quantile normalization is a widely used technique to correct for batch effects in DNA methylation and gene expression data. As this technique changes the data, it can affect variability patterns and trends in methylation levels. In our study, we have observed that, although quantile normalization has little effect on the ICC of differentially methylated CpG sites, it often lowers the ICC of low-variability loci. In addition, it changes the magnitude of the temporal trend estimates for CpG loci with significant trend, but does not change its direction. These observations suggest that for most likely targets of a DNA methylation study, such as differentially methylated loci or those with significant trend, applying quantile normalization as a batch correction technique is unlikely to significantly alter analysis results. However, for studies interested in the magnitude of temporal effects in methylation, this normalization technique can have a profound effect on the results and should be used with caution. Other batch correction techniques should also be considered (32).

Limitations of the present study include its modest sample size, which, nonetheless, is not unusual for epigenome-wide DNA methylation studies to date due to the cost constraint. Second, the absence of samples from men may restrict generalizability of our results. Although a recent study by Lam and colleagues (33) found only subtle sex differences in methylation of a small subset of CpG loci, one cannot discount the possibility that variability and temporal changes in some CpG loci may differ between men and women. Besides, there is conflicting evidence on the effect of sex on age-related changes in the methylome outside of sex chromosomes (23, 34). Third, in the absence of duplicate contemporaneous samples from the same subject and of technical replicates, it is unclear how much of the detected within-person variation is due to the temporal trend or other nontemporal factors, such as cell fraction differences. Also for this reason, we could not separate technical variability as part of the total variability in DNA methylation. Despite these limitations, we were able to examine variability patterns and detect short-term temporal methylation changes in a substantial number of CpG sites.

In summary, for CpG loci with differences in methylation between people, methylation levels can be reliably assessed with one blood sample; however, more samples and possibly special statistical methods are needed for low-variability and unmethylated loci. The X chromosome exhibits more stable methylation patterns, especially in CGIs and gene promoters, which is consistent with the effects of X chromosome inactivation. Although short-term temporal changes are largely driven by changes in the cell type composition of blood, trend unrelated to cell type was also detected in a small fraction of CpG sites. Further studies are needed to examine whether CpG loci with short-term temporal trend undergo similar methylation changes throughout lifespan and to what extent such changes are related to the onset or progression of disease.

No potential conflicts of interest were disclosed.

Conception and design: Y.B. Shvetsov, M.-A. Song, X.-O. Shu, H. Yu

Development of methodology: Y.B. Shvetsov, M. Tiirikainen

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Q. Cai, M. Tiirikainen, Y.-B. Xiang, H. Yu

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y.B. Shvetsov, M. Tiirikainen, H. Yu

Writing, review, and/or revision of the manuscript: Y.B. Shvetsov, M.-A. Song, Q. Cai, M. Tiirikainen, Y.-B. Xiang, X.-O. Shu, H. Yu

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y.-B. Xiang, X.-O. Shu

Study supervision: M. Tiirikainen, X.-O. Shu, H. Yu

Other (experiments: Illumina array, data quality check): M.-A. Song

The authors thank Wei Zheng for his leadership and support, Dake Liu, Gong Yang, and Charles Matthew for their contribution to the SPA, Regina Courtney and Jie Wu for assistance in DNA extraction, and Hui Cai for data preparation. Sample preparation was performed at the Survey and Biospecimen Shared Resource, Vanderbilt-Ingram Cancer Center.

This study was supported by NIH grants and contracts P30 CA068485 (Vanderbilt-Ingram Cancer Center), R37 CA070867 (to W. Zheng), and NO2-CP11010-66 (to X.-O. Shu).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.