Data pre-processing

Determination of technical versus biological variation

First to assess the sensitivity to detect biological variability between co-twins versus technical variation, we performed replicate hybridizations of three MZ twin pairs both at birth and 18 months. We compared the level of variation between co-twins (biological variation) to the level of variation between each technical replicate sample (technical variation). Biological variation (twin 1 versus twin 2) consistently exceeded technical variation (twin 1 versus twin 1; twin 2 versus twin 2) for each twin pair (Figure 1). We determined the average level of differential methylation between all biological and technical replicate arrays using a moderated paired t-test with false discovery rate correction. Precisely 230,340/330,155 probes were differentially methylated (adjusted P-value <0.05) across all biological replicates, whilst 858/330,155 probes were found to vary (adjusted P-value <0.05) across all technical replicates of twin pairs.

Determining relationships between samples

Unsupervised hierarchical clustering of the entire dataset (Figure S1A in Additional file 1) revealed that most samples cluster according to age. The majority of co-twins also cluster together: 7/9 (78%) MZ co-twins cluster at birth and 6/6 (100%) at 18 months, while with DZ co-twins, 4/5 (80%) cluster at birth and 2/4 (50%) cluster at 18 months. To explore the variation in this dataset attributable to the effect of sequence variation on methylation values via cis genetic effects or probe hybridization, we performed hierarchical clustering selectively for probes overlapping known SNPs as defined by the HM450 SNP manifest (version 3, 103,148 probes). The results compared well to the full dataset: 7/9 MZ co-twins cluster at birth and 6/6 cluster at 18 months, while 4/5 DZ co-twins cluster at birth and 3/4 cluster at 18 months (Figure S1B in Additional file 1). Restricting this analysis to probes with reported SNPs at the CpG site assayed by the probe (2,527 probes in this data set) resulted in 8/9 MZ co-twins clustering at birth and 6/6 at 18 months (Figure S1C in Additional file 1). Interestingly, on average for this set of probes, DZ twins did not cluster with their co-twin; rather, DZ twins at birth clustered with their matched samples at 18 months. Thus, data for such probes are likely to reflect the genotype of the individual rather than representing purely methylation levels. A random sampling of the same number of SNP-associated probes did not reproduce this clustering (data not shown), indicating this effect did not represent a sampling bias. These results suggest that SNP-containing probes account for little variation in the overall data set, with the exception of probes with SNPs at the CpG site assayed.

Identification of age-associated differentially methylated probes

To identify specific sites of differential DNA methylation associated with age, we used an empirical Bayes method [37] to compare birth samples with matched 18-month samples in all individuals and performed a probe-wise moderated paired t-test for differential methylation. Using this approach we found that 30.1% (99,198) probes changed significantly over time (adjusted P-value <0.05). These age-associated differentially methylated probes (aDMPs) changed by a mean β of 0.031 (3.1%) per year. Adding a further stringent cutoff of >20% absolute change over time to minimize technical effects [38] resulted in 0.8% (2,632) probes classified as stringent aDMPs (Table S1 in Additional file 2). Of these aDMPs, 87% showed a gain in DNA methylation over time whereas 13% showed reduced methylation (Figure 2a). We selected candidate aDMPs for validation based on their ranked change in methylation β value from birth to 18 months. The Sequenom MassArray Epityper platform was used to provide an independent measure of DNA methylation at aDMPs and confirmed the validity of the HM450 dataset. Using this approach, we confirmed that aDMPs identified by HM450 analysis are also representative of methylation at surrounding CpG sites (Figure 2b). Ontology and pathway analyses of the aDMP-associated genes showed an over-representation of cell development, morphogenesis (especially neuronal cells), and GTPase signaling pathways (Table 2; Tables S2 and S3 in Additional file 2). In order to determine whether aDMPs were more likely to occur at specific regions in the genome, we calculated the observed/expected frequency (enrichment) of genomic locations annotated in the HM450 manifest and assigned P-values with hypergeometric means tests. Intergenic regions were most likely to show changes in DNA methylation from birth to 18 months (Figure 3, grey bars; enrichment = 6.0×), followed by enhancers (2.5×) and 'open sea' regions >4 kb distant from CpG islands (1.7×). Promoters and CpG islands, but not their flanking shores and shelves, were less likely to show changed methylation over time (Figure 3; relative enrichment of 0.23× and 0.39×, respectively

Identification of age-associated differentially methylated regions

In order to identify larger regions of coordinated methylation change over time, we adopted a recently published differentially methylated region (DMR)-finding method [39]. This 'bump hunting' method identifies genomic regions in which clusters of consecutive CpG sites exhibit change over time in the same direction. Estimates were obtained for aDMRs by computing group medians and obtaining a value for the smoothed estimate that exceeds a t-statistic cutoff of 0.995. Using these criteria, we defined 897 aDMRs consisting of 4 or more consecutive probes changing in methylation between birth and 18 months (aDMPs). These aDMRs ranged in size from 33 to 1,698 bp. Twelve of these regions contained ten or more consecutive probes within approximately 1 kb of each other (Table S4 in Additional file 2 with an example shown in Figure 4). Of all aDMRs, 44% are located within 5 kb of a transcriptional start site, compared to 29% for aDMPs (Figure S2 in Additional file 1 and Table S4 in Additional file 2). Ontology analysis indicated that the aDMRs were significantly enriched for biological processes associated with cellular and organ development and in DNA binding (FDR <0.05; Table 3). As the sixth largest DMR (Table S4 in Additional file 2) and a representative of a number of the top age-associated ontologies, DNA methylation was validated at the cytoplasmic FMR1 interacting protein 1 (CYFIP1) gene in all samples using the Sequenom MassArray Epityper platform (Figure 2b).

Epigenetic discordance within twin pairs at birth and 18 months

Within-pair epigenetic discordance resulting from non-shared environmental factors has been postulated to underscore variation in phenotypic traits [40, 41]. We examined discordance in DNA methylation profile within twin pairs at birth and at 18 months of age. We calculated twin discordance as the absolute difference in β methylation values within pairs in birth samples and separately for 18-month samples. We ranked all probes according to average within-pair discordance at each age and performed 'ranked-list' ontology, which differs from 'gene-list' ontology in that there is no requirement for a predefined cutoff. All probes on the array were ranked by their scores for average within-pair discordance at each age (most discordant to least discordant), and the ranked list of probes was analyzed by the GOrilla bioinformatics tool [42] to identify ontology terms over-represented at the top of the list, compared with the bottom. We found that the most discordant genes at birth were consistently enriched for ontology terms associated with RNA metabolism, including spliceosome components and transcription factors (Table S5 in Additional file 2). At 18 months of age, the most discordant genes were associated with a similar set of gene ontologies as seen at birth (Table S6 in Additional file 2). The genes with discordant probes at both time points include a wide array of spliceosome components (for example, WDR83 and CWC22), zinc finger proteins (for example, ZNF267, ZBTB1, ZNF10), ribosomal proteins (for example, RPS26, RPL15, RPL12) and transcription factors (for example, MAML1, HOXB13).

We next investigated the distribution of DNA methylation discordance across genomic regions to determine whether discordance is more likely to occur at specific genomic locations. We have shown previously, using HM27 arrays, that median within-pair methylation discordance increased with increasing distance from CpG islands in three tissues (cord blood mononuclear cells, human umbilical vein endothelial cells and placenta) in both MZ and DZ twins at birth [7]. As the HM27 array focuses primarily on gene promoters and CpG islands, we repeated this analysis taking advantage of the diversity of genomic locations contained within the HM450 arrays. We calculated absolute within-pair discordance as before, and plotted probe discordance across genomic location at birth and at 18 months. The distribution of discordance values was consistent across all genomic annotations targeted on the array, with no evidence of regional enrichment (Figure S3A in Additional file 1). Similar results were observed selecting the top 10,000 most variable probes, or alternatively when the analysis was performed separately at birth and 18 months separately for both MZ and DZ twins (data not shown). We then filtered the dataset to include only probes present on the HM27 arrays and found evidence of higher levels of discordance around shores and shelves of CpG islands (Figure S3B in Additional file 1), which is consistent with our previously published observation with this platform [7].

Level of change in epigenetic discordance (drift versus convergence) over the first 18 months is a pair-specific phenomenon

Since previous cross-sectional studies suggest that epigenetic discordance in twins increases with age [30], we next investigated the degree of epigenetic drift from birth to 18 months of age within our twin pairs. The probe-wise level of within-pair discordance for CpG sites exhibiting a β-discordance value of greater than 0.2 (>20% discordant) was visualized at each age on scatterplots (Figure 5a, points in red). In contrast to the anticipated drift associated with age, we observed that the degree of within-pair discordance over time varies in a pair-specific manner (Figure 5a), with some pairs becoming more discordant in 18 month samples compared to birth samples (that is, epigenetic drift), some pairs becoming less discordant in 18 month samples, which we termed 'convergence', and others similarly discordant at both ages ('stable'). This was supported by Euclidean distance measures of twin discordance [7] (Figure 5b). These phenomena were not associated with zygosity or chorionicity, nor influenced by the effects of probes targeting SNPs on the array (Figure S4 in Additional file 1). We further calculated the change in discordance with age (delta discordance), as the difference in twin discordance (absolute values) from birth to 18 months. The distribution of delta discordance values was strongly centered about zero, with no evidence for overall skewing with age (Figure 5c). A comparison of the magnitude of the absolute values of differences in within-pair discordance over time (delta discordance) compared to the absolute values for methylation change over time (18 months - birth) indicated that age-related changes are far greater on average than changes to within-pair discordance (Figure 5d).

Epigenetic discordance at birth and 18 months of age

Very few genome-wide studies of methylation or expression have been performed on buccal cells. One such study, of buccal cells collected from 20 twins aged 13 to 14 years using a low resolution CpG island array, found no significant methylation differences within pairs [44]. However, a study of smoking-induced differential gene expression in buccal cells identified a differentially expressed network of genes with, at the hub, transcription factors REL and CREB [58], which are among the top 10% most discordant genes at birth and 18 months in our data (Tables S5 and S6 in Additional file 2). Despite the extensive longitudinal changes in DNA methylation described above, we found that, in general, probes located within genes associated with RNA metabolism (for example, spliceosome components) and control of gene expression (for example, transcription factors) were consistently more discordant within twin pairs at both birth and 18 months of age. Of interest, this class of genes has previously been shown to have altered levels of transcription in buccal cancer [59].

Epigenetic drift and convergence

In the current study, we found that a summed value (Euclidean distance) of epigenetic discordance across hundreds of thousands of loci can vary between and within pairs and can increase or decrease over time. In accordance with our genome-scale findings, a longitudinal study of DNA methylation at seven imprinted gene loci in buccal cells between birth and one year of age showed that inter-individual variation similarly increased, decreased or remained similar in singletons and that the direction of change could differ between individuals [45]. A longitudinal study of DNA methylation at three genes in buccal cells in 46 MZ and 45 DZ twin pairs found that methylation drifted in some pairs and converged in others over time [8]. Similar results were found for MZ and DZ twins and a role for genetic, shared and non-shared environmental factors, dependent on genomic location, in these longitudinal changes was postulated [8]. For MZ pairs, changes in within-pair discordance must be influenced solely by stochastic and non-shared environmental factors. Evidence for the latter comes from our previous studies of methylation in newborn twins [6, 7, 60] and from a cross-sectional study of DNA methylation in seven genes in whole blood from >200 MZ twin pairs aged 18 to 89 years [61]. Data from a longitudinal, genome-scale study of DNA methylation (using HM450 arrays) in whole blood from an independent cohort of young adults (aged 22 to 32 years) also provides evidence of genome-scale methylation drift and convergence defined by changes in Euclidean distance over time (Figure S5 in Additional file 1).

Epigenetic drift has been postulated to arise from the cumulative effects of (non-shared) environment and stochastic events [30, 62, 63], the latter influenced by epigenetic events such as promoter occupancy by transcription factors [64] and by errors made during the maintenance of DNA methylation profile following DNA replication [30, 63]. Recent studies suggest that epigenetic drift may also reflect differing rates of change of methylation among the population [65]. Furthermore, others have argued that epigenetic variability (or noise) is itself genetically programmed and has evolved to mediate some degree of plasticity (via canalization) [66]. In contrast, we suggest that 'convergence' may involve sites of methylation equalization between co-twins, possibly reflecting regression to the mean as a contributing factor. Regression to the mean is a phenomenon in which it is a statistical certainty that individual phenotypes, such as growth patterns [67], shift to the population mean over time [68]. This explains why twins with birth weight discordance become more similar over time [69] and can be understood in terms of twin-specific uterine-specific restrictions being replaced postnatally by a greater degree of shared environment [69–71]. Indeed, the twins in the current study had a median weight discordance [(Weight of the heavier twin - Weight of the lighter twin)/Weight of the heavier twin] of 13.3% at birth and 2.8% at 18 months. Although caution is needed with interpretations from a small sample size, we note that 'converging' pairs were more likely to start with a higher within-pair discordance (mean Euclidean distance = 375) than the drifting pairs (mean Euclidean distance = 295) (Figure 5b), although this difference did not reach significance (P = 0.11). Clearly, larger longitudinal twin-based studies are needed to further investigate factors contributing to epigenetic drifting and convergence over time.

Subjects, tissues and DNA extraction

Sample collection from twins at the time of delivery was carried out with appropriate human ethics approval from the Royal Women's Hospital (project number 06/21), Mercy Hospital for Women (project number R06/30), and Monash Medical Centre (project number 06117C), Melbourne and the study was conducted according to the Declaration of Helsinki principles. The twin pairs chosen for methylation array analysis are shown in Table 1. The 10 MZ pairs and 5 DZ pairs shared a similar sex ratio, gestational age and birth weight to the full group of 250 pairs. Buccal cells were collected with Catch-all Sample Collection Swabs (EPICENTRE Biotechnologies, Madison, WI, USA) and were stored at -20°C until DNA extraction, which was performed as previously described [60].

Infinium HumanMethylation450 BeadChip data acquisition and processing

DNA samples (1 μg) were bisulfite converted using the Methyl EasyXceed bisulfite modification kit (Human Genetic Signatures, North Ryde, Australia), according to the manufacturer's instructions. Conversion efficiency was assessed by bisulfite-specific PCR. DNA samples were hybridized to Illumina Infinium Human Methylation450 (HM450) BeadChip arrays according to the manufacturer's instructions. Raw intensity data (IDAT) files were imported into the R environment (version 2.14.1) [74] and processed using the minfi package [75]. All analyses were performed in R using packages available from the Bioconductor project [76]. Data quality was assessed in minfi using plots derived from various control probes on the array. Poor performing probes defined as those with an average detection P-value >0.001 in one or more samples were removed from the analysis (n = 132,113). Data from five samples with an average detection P-value >0.05 and with evidence of poor bisulfite conversion efficiency were removed completely. Probes on the × and Y chromosomes were also discarded from all samples. The resulting data were pre-processed using the Illumina method within minfi and subset-quantile within-array normalization was performed [36] for combined normalization of Infinium type I and type II probes. The log2 ratio of methylated probe intensity to unmethylated probe intensity were calculated in minfi and the resulting M-values [77, 78] were quantile normalized between arrays using the limma package [79]. Sample quality was further assessed using hierarchical clustering plots available in minfi and lumi [77] packages. Following this, three additional samples were removed as outliers constituting a final data set of 330,168 probes and 53 samples.

Statistical analysis

Exploratory analysis of sample relationships was performed using unsupervised hierarchical clustering analysis with the Euclidean distance and complete linkage algorithm, and dendrogram was created using gplots [80]. Differential methylation analysis was performed on M-values using the limma package using a cutoff of FDR-corrected P-values <0.05 [81] and delta beta values >0.2. To study discordance among co-twins at the probe-level, a linear model was fitted to the M-values with twin-pair as a predictive factor to model the twin relationship. The level of discordance among co-twins was interpreted as the residual measurement for each CpG from the model-fit. For enrichment analysis, gene sets were populated with probe IDs using the annotated regions provided in the Illumina HM450 manifest file (version 1.1). Annotations used were classified as gene-related (TSS1500 and TSS200, regions from -1500 to -200 and -200 to the transcriptional start site respectively, 5' UTRs, first exons, gene bodies, 3' UTR and intergenic (no gene annotation)); CpG island-related (islands (also split into intragenic and intergenic)), shores (0 to 2 kb flanking islands), shelves (2 to 4 kb flanking islands) and open sea (>4 kb from islands) [82]); DMRs (associated with cancer (CDMRs) and induced pluripotent stem cell reprogramming (RDMRs); [83] and regulatory regions (promoters, enhancers and DNAse hypersensitivity sites, likely to be a mixture of promoters and enhancers [84, 85]). Boxplots were produced to graph each category by discordance score. The 'bump-hunting' methods described by Jaffe and colleagues [39] were implemented using the charm package available in Bioconductor [86]. We used the 'dmrFinder' algorithm without covariate adjustment, using the default SPAN settings and specifying a minimum four probes, and a t-statistic cutoff to identify probes as being in a DMR at 0.995. For gene ontologies the GOrilla bioinformatics tool [42] was used to perform ranked-list ontology using the entire array content ranked by scores for discordance. Gene-list ontology enrichment was performed on significant gene lists (FDR <0.05) using the DAVID bioinformatics tool under the default settings [87]. Pathway analysis data were analyzed through the use of Ingenuity Pathway Analysis (Ingenuity Systems, Redwood City, CA, USA). The analysis tool GREAT (Genomic Regions of Annotations Tool) [88] was used to analyze the functional significance of aDMRs using the single nearest gene association rule within a 100 kb window.

Sequenom MassArray target validation

Target validation was performed using the Sequenom MassArray EpiTYPER (Sequenom, San Diego, CA, USA) performed as previously described [18, 60]. Amplicons were designed using Sequenom EpiDesigner software. Primers are listed in Table S7 in Additional file 2. In brief, amplification was performed after bisulfite conversion of genomic DNA with the MethylEasy Xceed bisulphite conversion kit (Human Genetic Signatures, North Ryde, Australia). All PCR amplifications and downstream processing were carried out at least in duplicate and the mean methylation level at specific CpG sites determined. Raw data obtained from MassArray EpiTYPING were cleaned systematically using an R-script to remove samples that failed to generate data for more than 70% of CpG sites tested [60]. Also, technical replicates showing ≥10% absolute difference from the median value of the technical replicates were removed and only samples with at least two successful technical replicates were analyzed. Samples were compared across each analyzable CpG site in the amplicon, as well as the mean across the whole amplicon.

Data availability

Array data described in this manuscript have been submitted to the Gene Expression Omnibus public repository and are freely available under the accession number GSE42700.