Genome-wide impact of hydrogen peroxide on maintenance DNA methylation in replicating cells

Background Environmental factors, such as oxidative stress, have the potential to modify the epigenetic landscape of cells. We have previously shown that DNA methyltransferase (DNMT) activity can be inhibited by sublethal doses of hydrogen peroxide (H2O2). However, site-specific changes in DNA methylation and the reversibility of any changes have not been explored. Using bead chip array technology, differential methylation was assessed in Jurkat T-lymphoma cells following exposure to H2O2. Results Sublethal H2O2 exposure was associated with an initial genome-wide decrease in DNA methylation in replicating cells, which was largely corrected 72 h later. However, some alterations were conserved through subsequent cycles of cell division. Significant changes to the variability of DNA methylation were also observed both globally and at the site-specific level. Conclusions This research indicates that increased exposure to H2O2 can result in long-term alterations to DNA methylation patterns, providing a mechanism for environmental factors to have prolonged impact on gene expression. Supplementary Information The online version contains supplementary material available at 10.1186/s13072-021-00388-6.


Background
Epigenetic modification of chromatin provides a mechanistic basis through which environmental stimuli can modulate gene expression. Epigenetic marks, such as methylation of cytosine and histone modifications, are utilized during development to generate a wide range of cellular phenotypes from the same DNA code. While the majority of epigenetic changes in adult organisms are due to the continued differentiation of stem cells, they can also be modified during ageing [1,2] and by interactions with environmental stimuli such as cigarette smoke, pollutants, and alcohol consumption [3][4][5][6][7][8][9]. However, there is a lack of mechanistic studies that investigate how environmental factors impact DNA methylation.
Several potential environmental modifiers are associated with increased oxidative stress, including aging and chronic inflammation [10,11]. Mitochondrial dysfunction is linked to increased superoxide and hydrogen peroxide (H 2 O 2 ) production in a range of pathologies [12]. Ageing, smoking, radiation and xenobiotic metabolism can also result in increased H 2 O 2 in biological systems [13][14][15][16][17][18][19]. Immune cells such as neutrophils and macrophages generate a range of reactive oxygen species following NADPH oxidase (NOX) activation, many of which are cell permeable, and contribute to oxidative stress during inflammation [11,20,21]. H 2 O 2 and myeloperoxidase-derived oxidants produced by immune cells can oxidize redox-sensitive thiol proteins and methionine [22][23][24]. Importantly, methionine is a precursor of the methyl donor S-adenosyl-methionine (SAM).
We have previously reported that H 2 O 2 and chloramines generated by phagocytic immune cells are able to inhibit DNA methyltransferase (DNMT) activity during DNA replication [25]. These results suggested that this inhibition was due to the oxidation of DNMT1's catalytic site cysteine and occurred at concentrations that did not have major impacts on cell proliferation. DNMTs copy DNA methylation to the nascent DNA strand during DNA replication, a process which ensures that the correct pattern of methylation is maintained throughout cellular division. However, in a cell culture system, only the chloramines had a significant impact on global DNA methylation because they also depleted SAM [25].
In this study, we have revisited the impact of H 2 O 2 on maintenance DNA methylation in cultured cells by using a more sensitive genomic approach that individually monitors 850,000 CpG sites across the human genome. We have also investigated the reversible nature of any modifications analysing changes seen at early and late timepoints. We reveal that sublethal doses of H 2 O 2 caused a significant disturbance to DNA methylation patterns, increasing overall variability and site-specific methylation in the 1-2 h immediately following exposure to this oxidant. Many of the changes were lost 72 h after exposure, but not all. Collectively, these findings describe a potential mechanism for various environmental factors and pathologies that promote oxidative stress to impact the methylation pattern of proliferating cells.

Results
We propose that oxidative stress is able to inhibit the maintenance methylation that occurs following DNA replication in proliferating cells. Cultured Jurkat T-lymphoma cells were synchronized with a thymidine block to maximize the number undergoing DNA synthesis at the time of exposure to oxidative stress. Flow cytometry was used to monitor cell cycle transitions and determine the period of maximal DNA replication during the S-phase (Fig. 1a). The asynchronous population of Jurkat cells shows that just before cell division the majority of cells were in the G 0 /G 1 phase, ~ 20% were in S-phase and ~ 14% of cells were in G 2 /M phase. The period of 2-5 h post-release was when the greatest number of cells entered the S-phase (Fig. 1b) [26,27]. Cells were therefore exposed to two doses of H 2 O 2 , at 2 and 3 h post-release, to maximize the potential for inhibition of maintenance methylation during DNA replication. To determine the highest concentration of H 2 O 2 that could be used without significant cell death, synchronized cells were exposed to a range of concentrations (Fig. 2a). Treatment with 15 µM H 2 O 2 resulted in a large decrease in viability after 24 h, while 5 µM and 10 µM H 2 O 2 treatments only had a small effect. With 10 µM H 2 O 2 , which was selected for subsequent experiments, there was no significant difference in viability between the treatment and control groups at 4 or 72 h (p = 0.25) (Fig. 2b). Assessment of cell growth indicated that only the control cells doubled in cell density after 24 h, whereas the H 2 O 2 -treated cells had significantly lower percentage growth (p = 0.02) (Fig. 2c). Cells were diluted to a concentration of 0.25 × 10 6 cells/ mL at 24 h, and the decrease in fold growth for the treatment samples compared to control was still observed at 48 and 72 h, although not statistically significant (p = 0.1 and 0.6). This slowing of growth rate was confirmed by the analysis of CFSE dilution by flow cytometry (Fig. 3). Collectively this data shows that both control and treated cells underwent at least one round of division by the 72-h time point.

DNA methylation analyses
DNA methylation profiles were assessed using the Illumina Infinium MethylationEPIC 850K array at 4 h and 72 h after the addition of cytidine to release cells from thymidine block, which was 2 and 70 h from the first H 2 O 2 exposure, respectively. This technology assesses DNA methylation at over 850,000 CpG sites across the human genome, with nucleotide-level resolution.

Principal component analysis (PCA) of DNA methylation and hierarchical clustering
Multidimensional scaling was used to assess the top contributing factors in data variability, which corresponded with treatment, and to a lesser extent time (Fig. 4a). Overall, there was a distinct separation between samples from the control and H 2 O 2 treatment groups. Samples from the control group formed a relatively tight cluster with less spatial separation, and as expected there was no separation between time points for the control samples. The H 2 O 2 treatment samples from the 72-h time point more closely represented the control samples than the H 2 O 2 treatment samples from control at the 4-h time point.
Hierarchical clustering was performed on a distance matrix calculated using the β-values for all probes, regardless of significance. Except for one sample, which incorrectly grouped within the treatment, this approach successfully differentiated between the control and H 2 O 2 treatment groups (Fig. 4b).

Relative average DNA methylation change
At the 4-h time point, the mean methylation level of the H 2 O 2 treatment group was 3.33% less than the control group (p = 0.013, t (5) = 3.72) (Fig. 5a). However, at 72 h this significant difference was no longer observed (p = 0.10, t (4) = 2.2). There was no significant difference between the 4-h and 72-h time points within each group, although the H 2 O 2 treatment group demonstrated a minor increase and had substantially reduced variability (Fig. 5a).
Plotting the individual probe-wise comparisons of significance versus effect size change confirmed there was a substantial decrease in methylation after H 2 O 2 treatment , with several probes from the 4-h time point demonstrating a significant log 2 FC greater than ~ 1. There was a comparative reduction in the effect size and significance at the 72-h time point, which indicates that the H 2 O 2 treatment and control samples were more similar to each other than at the 4-h time point (Fig. 5b, c). Together these observations suggested that on a global scale, there were significant decreases in DNA methylation in the treatment group compared to the control group at 4 h, which were largely restored by 72 h (Fig. 5).

Analysis of variability in DNA methylation after H 2 O 2 treatment
Principal component analysis and global methylation change suggested that there was a higher level of heterogeneity in the DNA methylation profiles after H 2 O 2 treatment. To investigate the individual probes that contributed towards this observation, we assessed the differential variability at individual CpG positions Forty-six DVPs were common between the two time points, however, only 18% (7 DVPs) showed a consistent direction of change (Fig. 6).

Differentially methylated CpGs at the 4-h time point
To identify changes that associated with H 2 O 2 treatment we assessed DNA methylation levels at each individual probe using the "topTreat" algorithm, within the R package, Limma. To summarize the results, each probe was mapped to the corresponding genomic region and  unadjusted significance was plotted across the whole genome as a Manhattan plot (Fig. 7). This analysis demonstrated that the significant changes were distributed across the genome, and did not appear to be concentrated at specific genomic locations. There were 1162 probes that demonstrated a significant decrease in methylation using a log 2 FC (M-values) cutoff of 1, and 89 probes that demonstrated a significant increase in methylation (adj. p < 0.05). Increasing the log 2 FC cutoff to 1.2 identified 87 probes that demonstrated a significant decrease in methylation and three probes that demonstrated a significant increase (Additional file 1: Table S1). The top most significant probes had an initial methylation value of 20-25% in the control group, and a value of 10-15% in the H 2 O 2 treatment group (Fig. 8a). Therefore, the raw data for the top six most significant probes that did not conform to this general trend are presented in Fig. 8 (including the three probes that demonstrated increased methylation). The top 20 significantly differentially methylated positions with a log 2 FC > 1.2 in order of largest fold change are presented in Table 1.

Differentially methylated CpGs at the 72-h time point
At the 72-h time point, 19 probes demonstrated a significant decrease in methylation and one probe demonstrated a significant increase in methylation (adj. p < 0.05), when assessed using a log 2 FC cutoff of 1 (Additional file 1: Table S2). The magnitude of these changes was less than 10%, therefore there were no significant probes detected at a log 2 FC cutoff of 1.2 at the 72-h time point.
Seventeen of the 20 significant probes detected at the 72-h time point were also detected in the 4-h time point when a log 2 FC of > 1 was used (Table 2) and the raw data from a subset of these probes are presented in Fig. 9. Two probes were also detected in the 4-h time point when a log 2 FC > 1.2 was used (cg05577994 mapping to gene: EMX2OS at chr10:119,254,968, adj. p < 0.05 and cg18363176 mapping to gene: FLI1 at chr11:128,606,110, adj. p < 0.04). Both probes had a smaller log 2 FC at 72 h than at 4 h. However, all other significant probes showed a substantially larger log 2 FC at 72 h than at 4 h ( Table 2).

Top differentially methylated gene regions
We next investigated if multiple probes mapping to the same genomic loci demonstrated consistent methylation changes. We limited this analysis to differentially methylated regions (DMRs) that contained five or more adjacent significant CpGs, with an average change in methylation larger than 5%. Under these parameters, we detected 70 significant DMRs at the 4-h time point and five significant DMRs at the 72-h time point (Additional file 1: Tables S3 and S4). Two significant DMRs were consistent between the 4-h and 72-h time points (Table 3 and Fig. 10).
The genes corresponding to the top most significant CpG sites at 4 h, were assessed for significant biological relevance by pathway analysis within the GO and KEGG databases. Neither comparison identified significantly enriched pathways, after correction for multiple testing (Additional file 1: Table S5).

Discussion
We used methylation arrays to investigate if exposure to H 2 O 2 during DNA replication could interfere with the transfer of DNA methylation patterns in proliferating cells. H 2 O 2 exposure during the S-phase of DNA replication disrupted DNA methylation of the nascent strand in proliferating cells 2 h after the first exposure to H 2 O 2 , with most but not all of these changes reversed by 72 h. By 72 h, all cells had undergone at least one round of division, indicating that these changes were representative of the epigenome of daughter cells. Sublethal H 2 O 2 is known to alter gene expression in the short-term through impacting the activity of redox-sensitive transcription factors. Our results provide a mechanism by which H 2 O 2 may have long-term impacts on gene expression and cell function.
Because aliquots of the treatment and control samples were split from the same source immediately prior to treatment, methylome changes relating to cell cycle were controlled for in the study design. The design of these experiments was such that DNA methylation was assessed 1-2 h after treatment, thus comparing the efficiency of nascent DNA methylation between the control and treatment samples. DNA methylation was then assessed again at 72 h to establish if the observed changes were corrected, or were inherited into the epigenome of subsequent cell generations. Cell viability and live cell counts indicated that approximately 72 h was required before cellular replication in the treatment group resembled that of the control group. This is likely to reflect the length of time required for repair of DNA damage, and re-establishment of the epigenome, and therefore DNA methylation was not assessed at the 24-and 48-h time points. An unsupervised assessment of data variability indicated that the treatment samples from the 4-h time point grouped distinctively from all other samples. There was no noticeable separation between the control samples, which more closely resembled the 72-h treatment samples.
Analysis of the average relative DNA methylation levels between groups indicated that H 2 O 2 exposure corresponds with a significant decrease of ~ 3.3% in average methylation between treatment and control at the 4-h time point, however, there was no significant change between the 72-h time points. This change was driven by the combined effect of many probes with a log 2 FC approximately equal to one, however, at 72 h the effect sizes (and therefore significance) were substantially reduced. Hierarchical clustering was relatively successful at distinguishing the control samples from the treatment samples, and also identified that samples from the same biological replicate were more closely related, an effect which was subsequently controlled for in the statistical design.
These findings build on our previous work where we showed that while H 2 O 2 inhibited DNMT, there was no detectable change in global methylation levels [25]. The main experimental difference our previous work was that H 2 O 2 exposure occurred at the time of cytidine release. We modified our previous protocol by treating the cells twice with 10 µM H 2 O 2 during the S-phase of DNA replication, at 2 and 3 h after the replication block  was released, in order to maximize H 2 O 2 exposure during the time period where DNMTs are most active [28]. Furthermore, while our previous study utilized a mass spectrometry approach which analysed the methylation of incorporated heavy-labelled deoxycytidine in the nascent strand, the current study used array technology that increased both the magnitude of the effect size and sensitivity, enabling detection of very small changes that were not previously possible using traditional methods. Our current work also highlighted the importance of nucleotide resolution, which enabled the detection of both increased and decreased methylation.
Intriguingly, there was a substantial increase in differential variability at both global and site-specific measurements. Differential variability of methylation at for the 4-h time point was significantly increased for the H 2 O 2 treated samples compared to controls. Theoretically, his could either be due to increased variability at specific sites or random sites across the genome. We found that 1035 probes had significantly increased variability associated with H 2 O 2 treatment at 4 h, which was reduced (41 probes) by 72 h. Therefore, site-specific increases in variability were contributing to the global variability in our system. It is important to note that there was also low consistency in the direction of change between the two time points, where 82% of the probes that were significant at both time points showed contrasting results.
The role of stochasticity in biological systems has received increasing attention in multiple contexts [29]. Whether cell-intrinsic or extrinsic, stochasticity can impact on a wide range of cellular processes including epigenetics and cell differentiation. Multicellular organisms have a remarkable capacity to produce consistent phenotypic outcomes even when challenged by variable conditions [30,31]. This "buffering capacity" of normal cells is dysregulated in cancer and can be unmasked by treatment with a toxic stress such as chemotherapy. In practical terms, epigenetic heterogeneity in response to an environmental stress enables a range of transcriptional states across a population, some of which might confer resistance. This has been observed in acute myeloid leukaemia where epigenetic heterogeneity has been correlated with poorer outcomes and a shorter time to relapse [32]. Our finding that treatment with H 2 O 2 increased the variability of methylation at the early time point is consistent with this and potentially points to a mechanism that allows a proportion of cancer cells across a population to adapt to oxidative stress. We also investigated significant differentially methylated sites that demonstrated the largest fold change after H 2 O 2 treatment. At the 4-h time point, 1162 probes showed a significant reduction in DNA methylation (log 2 FC > 1) in the treatment samples compared to control samples. In addition to this, 90 probes had an effect size larger than log 2 FC of 1.2, which suggests that in these instances the methylation level of the parental DNA strand may have been influenced [33]. With the exception of two probes (mapping to EMX2OS and FLI1), these 90 probes did not show a significant change at the 72-h time point. This suggests that the large changes observed at the 4-h time point were eventually corrected, or occurred in cells that eventually died. The majority of the significant probes at the 72-h time had a smaller effect size at 4 h (log 2 FC ~ 1), which suggests that they represent a delayed effect. Although no significant enrichment of gene pathways were observed, several of the top most differentially methylated genes are disrupted in pathologies such as diabetes and cancer [34][35][36][37][38][39][40][41][42][43][44][45].
Two gene regions corresponding to the SPATA16 gene and LOC100506102 (also overlapping with the small nuclear RNA: snoU13.410-201, and the miRNA: RP11-322J23.1-001) demonstrated significant differential methylation at both 4 h and 72 h. The potential results of these changes are unclear; however, the SPATA16 gene is related with sperm count and motility. Extensive research has suggested that reactive oxygen species are a contributing factor to male infertility [46], with oxidative stress identified as a common pathology in approximately half of all infertile men [47]. Deregulation of snoRNA has been observed in several different cancer types including chronic lymphocytic leukaemia [48], hepatocellular carcinoma [49], colorectal cancer [50] and endometrial cancer [51], and also metabolic stress disorder [52]. Their roles as diagnostic and prognostic biomarkers have also been extensively studied [50,53,54].
Oxidants can directly damage and modify DNA, leading to mutations and genetic instability, but they can also cause interference in epigenetic pathways. One of the most common products of cellular oxidative stress is oxidized guanine or 8-hydroxy-2′-deoxyguanosine (8-OHdG) [55][56][57]. Increased frequencies of 8-OHdG in hepatocellular carcinoma patient DNA have been associated with increased oxidative stress, and in vitro studies have revealed that H 2 O 2 exposure (50-250 µM) increased the occurrence of 8-OHdG and corresponded with altered methylation patterns in promoter regions of tumour suppressor genes [58]. When 8-OHdG is formed in close proximity to CpG sites, it has been shown to interfere with DNMT1 activity by impeding its access to the nascent strand, resulting in impaired maintenance methylation [59]. There is also evidence that 8-OHdG, through the recruitment of 8-oxoguanine DNA glycosylase-1, activates demethylation by ten-eleven translocation methylcytosine dioxygenase 1 (TET1) [60]. However, one study has suggested that the cumulative occurrence of 8-OHdG is low [61], and no reports have demonstrated specific enrichment at CpG sites. Therefore, 8-OHdG formation due to H 2 O 2 exposure, particularly at low concentrations, would be unlikely to account for the large genome-wide decreases in methylation observed in our study. Another potential mechanism for regulation of methylation levels by oxidants is decreased TET activity. TET enzymes are sensitive to redox regulation, as their iron-dependent catalytic regions require ascorbate for optimal activity [62][63][64][65]. Oxidants such as H 2 O 2 , chloramines and superoxide, ubiquitous in an environment of oxidative stress, have been shown to deplete ascorbate levels, and could contribute to reduced TET activity [66,67]. Alternatively, since there have been several reports that H 2 O 2 can directly inhibit DNMT1 activity, it is more likely that the transient inhibition of DNMT1 due to H 2 O 2 exposure is largely responsible for the decreased methylomic changes observed [25,68]. However, the methylation changes from the 4-h time point were spread across the genome, with no strong evidence for enrichment of specific gene pathways, and yet sequence-specific effects on DNA methylation were observed across four independent experiments. This result would not be anticipated if H 2 O 2 were acting through inhibition of DNMT1 activity alone. One possibility is that since many of the changes were consistent across multiple adjacent probes, the effects of H 2 O 2 treatment relate to the positioning of DNMT1 at the time of exposure. The cells were released from synchronization at the same time, therefore they would progress through the cell cycle at a similar rate. This may indicate that there are undiscovered mechanisms that combine with inhibition of DNMT activity to yield sequence-specific alterations in DNA methylation.
Most of the large epigenetic changes reported in the literature are associated with long-term exposure to environmental stressors such as smoking and chronic inflammatory disease states [5,[69][70][71][72], it was therefore unexpected to observe a large number of methylation changes with a relatively short oxidant exposure period. However, our findings have a substantial effect size, and the oxidant concentration used in this analysis is similar to in vivo concentrations observed during chronic infection [73][74][75][76]. Although it is unlikely  that a large population of cells would exist in S-phase at any one time point in vivo, in the situation of prolonged inflammation, a large proportion of cells will be exposed to a range of oxidants during cellular replication. Therefore, it is conceivable that H 2 O 2 -induced inhibition of DNMT could also occur in an asynchronous population of replicating cells in an inflammatory environment. Furthermore, although the current study was limited by the utilization of a single cell type, DNMT inhibition due to H 2 O 2 treatment has been reported in other cell lines [68], and therefore methylation changes as a result of oxidative stress should be futher explored in other cell lines as well as primary cells.
These data provide supporting evidence for a mechanistic link between environmental oxidative stress and DNA methylation. We propose that oxidative stress can influence DNMT1 activity, causing longer term changes in cytosine methylation patterns. Although it is accepted that environmental factors can influence the epigenome, few studies have investigated the pathways through which this can occur and the effect in subsequent cell generations. These data demonstrate that exposure to oxidative stress during cellular replication corresponds with substantial decreases in DNA methylation, and while the majority of these changes are subsequently corrected, a smaller number of sites remain demethylated. Furthermore, exposure to H 2 O 2 resulted in increased variability of methylation in a site-specific manner. This is of great interest for cancer biology where metabolic reprogramming and the tumour microenvironment combine to create a milieu of oxidative stress that can promote altered epigenetic profiles and modify cell behaviour of both neoplastic and normal cells.

Conclusions
The results of this research indicate that the human Jurkat T-lymphoma epigenome is susceptible to sustained changes in DNA methylation after exposure to oxidative stress during cellular replication. These studies need to be expanded to other proliferating cells. H 2 O 2 exposure was associated with genome-wide decreases in DNA methylation, which were largely corrected during subsequent cycles of cell division. The magnitude of the observed changes generally decrease with time, however, a subset of changes that displayed a relatively small initial effect size at 4 h demonstrated an increased effect size at 72 h. These results suggest that oxidantinduced changes in DNA methylation can be inherited by subsequent cell generations. Furthermore, we show that treatment with H 2 O 2 affects variability of methylation in a global and site-specific manner.

Cell synchronization
Cellular replication of E6.1 Jurkat cells was synchronized using a single thymidine block protocol to arrest cells between late G 1 and early S-phase. Briefly, cells (1 million cells/ml) were treated with thymidine (2′-deoxythymidine, Sigma) at a final concentration of 1 mM and incubated for 18 h. To promote release into the S-phase, the block was alleviated by two media changes and addition of cytidine (2′-deoxycytidine, Sigma) to a final concentration of 50 µM. Cell cycle progression was assessed by flow cytometry every two hours over a 12-h period by analysing cellular DNA content following cell fixing in ethanol and staining with propidium iodide (PI) (Sigma). This approach was used to detect cells in each of the three major phases of the cycle (G 0 /G 1 , S and G 2 /M), and also allowed for detection of apoptotic cells with minimal DNA content.

Cell culture procedure
Each of four replicates originated from frozen cells of differing passage lineages, and the following procedure was performed for each replicate, on independent days. Cellular synchronization was performed on 200 × 10 6 cells in 200 mL as described above. Immediately after release, the cells were split into two equal volumes at 0.5 × 10 6 / mL for treatment and control samples. At two hours post-release, cells in media were treated with H 2 O 2 at a final concentration of 10 µM. This was repeated one hour later, meaning the cells received two H 2 O 2 bolus doses totaling 20 µM. Cells were split to a concentration of 250,000 cells/ mL after 24 h. Samples containing approximately five million cells, were collected at four hours and 72-h post-release for DNA extraction. Flow cytometry for cell viability (using PI), cell cycle analyses, and live cell counts were performed at pre-block, pre-release, and 4, 24, 48 and 72-h post-release time points.

Cell viability, proliferation and cell cycle progression
The percentage of viable cells at each major time point was assessed by the exclusion of PI. Cell cycle analyses were performed by fixation with ethanol followed by incubation with PI for 30 min. Cell proliferation was observed by labelling cells with carboxyfluorescein diacetate succinimidyl ester (CFSE; Invitrogen, New Zealand), as described previously [77].

DNA extraction and bisulfite conversion
DNA was extracted from cell samples using a Gene-Jet DNA extraction kit (ThermoFisher Scientific, USA) according to the manufacturer's instructions. Sodium bisulfite treatment of genomic DNA was performed on 1000 ng of DNA using the Zymo Research EZ DNA methylation kit, according to the manufacturer's recommended protocol for use on the Illumina Infinium MethylationEPIC 850K array. Cells harvested at 4 h and 72 h after cell cycle release were assessed for DNA methylation profiles using the Illumina methylationEPIC bead chip array. Methylated and non-methylated values were determined using the minfi pipeline, within R. To control for DNA methylation changes that occur naturally during cellular replication, the block-release protocol was performed on all cells together.

Bioinformatics analyses
DNA methylation profiles of the control samples were subtracted from the treatment samples. Data analyses were performed using the Minfi and Limma Bioconductor software packages within the R statistical program (www.R-proje ct. org). All packages, programs and pipelines used in these analyses are freely available, and the workflow was based upon published scripts [78]. Any additional scripts are available from the authors on request. All samples passed quality control together using Subset-quantile Within Array Normalization (SWAN) [79] and variance-stabilizing transformation [80,81]. Probes containing detection p-values > 0.05 for 1% or more samples, were excluded from further analysis. Because the cell samples were of the same sex, all chromosomes were retained, however, probes identified as having polymorphic hybridizing potential and homology to common SNPs [82] were removed. The final data frame contained 809,911 probes available for analysis. Multidimensional scaling of the top 1000 methylation values was performed using pairwise distance method for gene selection, and normalized, filtered methylation values. Hierarchical clustering was performed on a "minkowski" distance matrix calculated using the β-values for all probes, regardless of significance [83].

Average relative DNA methylation changes
Average relative DNA methylation changes were calculated using mean β-values from all probes for each sample in the four groups. The absolute mean differences between control and treatment groups at each time point were calculated using paired t-tests.

Differential variably
Differentially variable positions (DVPs) were identified using the DiffVar algorithm [84] within the missMethyl package [85]. This analysis used a similar linear model as below, however, biological replicates were incorporated into the statistical design as covariates.

Identification of differentially methylated CpGs
The methylation status of each probe was calculated using normalized probe signals represented as methylation values (M-values) and β-values. M-values were generated within Minfi as the log 2 ratio of the signal intensities of methylated probe divided by the unmethylated probe. Unless stated otherwise, statistical analyses were performed using the M-values. β-values (average DNA methylation level for each probe) were used for data visualization and range from 0 (unmethylated) to 1 (methylated). β-values were generated by dividing the methylated probe signal with the sum of the methylated and unmethylated probe signals [86]. Normalized β-values and M-values were manually assessed for fit (Additional file 1: Figure S1).
Differentially methylated positions which correlated with treatment were identified using a linear regression model within the Limma package, with adjustment for multiple testing. This was calculated by comparing the methylation measurements of the control samples with the treatment samples, for each of the two time points (4-h post-release, and 72-h post-release). Samples that originated from the same cell passage were treated as a biological replicate and the correlation coefficient between these samples was incorporated into the statistical design using Limma's duppcorr function. The top most significant, differentially methylated CpG positions were identified by log fold change (logFC) weighted functions using Limma's topTreat algorithm. Adjustment for multiple correction was performed using the "Benjamini, Hochberg" method within Limma and Quantile-Quantile plots for the linear regression model were assessed at each time point using the observed against expected -log 10 (p-value) (Additional file 1: Fig. S2). The observed data show a minor inflation of p-values smaller than expected by chance. Therefore, we substantially reduced our criteria for significance by excluding probes that displayed less than a 10% mean difference in methylation between treatment and control samples. This also removed significant changes observed at sites that displayed methylation values less than 10% or greater than 90%, as changes within these ranges are less likely to be of biological relevance.
Pathway analysis was performed by comparison with the Kyoto Encyclopedia of Genes and Genomes (KEGG) [87] and gene ontology databases (GO) [88] using the missMethyl R package [85], with correction for probe bias.

Differentially methylated regions
Differentially methylated regions were interrogated within the Minfi package using the statistical package DMRcate [89]. A methylation differential cutoff of 10 was used, and unless stated otherwise significance was determined using a false discovery rate (FDR) of 0.05 in conjunction with a p-value cutoff of 0.05.