Genome-wide hydroxymethylcytosine pattern changes in response to oxidative stress

The TET enzymes convert methylcytosine to the newly discovered base hydroxymethylcytosine. While recent reports suggest that TETs may play a role in response to oxidative stress, this role remains uncertain, and results lack in vivo models. Here we show a global decrease of hydroxymethylcytosine in cells treated with buthionine sulfoximine, and in mice depleted for the major antioxidant enzymes GPx1 and 2. Furthermore, genome-wide profiling revealed differentially hydroxymethylated regions in coding genes, and intriguingly in microRNA genes, both involved in response to oxidative stress. These results thus suggest a profound effect of in vivo oxidative stress on the global hydroxymethylome.

Despite the above-mentioned advances, the impact of oxidative stress on TET-mediated hydroxymethylation is therefore unclear. We have addressed this issue using two experimental models, namely SY5Y neuroblastoma cells treated with buthionine sulfoximine (BSO), one of the most commonly used drug to induce oxidative stress, and Glutathione peroxidase 1 and 2 (GPx1/2) double-knockout mice. Our results show, for the first time, that both exogenous and in vivo oxidative assaults deeply modify the hydroxymethylome. Of particular interest, we unveiled an unexpected link between oxidative-stress-induced hydroxymethylation pattern changes, a set of microRNAs, and oxidative-stress-related genes.

Results
BSO-treated SY5Y cells and GPx1/2 double-knockout mice have a decreased level of hydroxymethylcytosine. To see whether hmC levels are affected in cells exposed to exogenous oxidative stress, we treated SY5Y cells with BSO and quantified hmC by dot-blot experiments and mass spectrometry (LC/MS/MS) (Fig. 1A,B and supplementary Fig. S1C). SY5Y is a human cell line commonly used as a model for studying oxidative stress and oxidative-stress-related diseases such as Alzheimer's and Parkinson's 11,12 . Additionally, neuroblastomas are derived from neural crest cells, and neuronal cell types have high levels of hmC 13 , making those cells a good model for our study. BSO inhibits the synthesis of the antioxidant glutathione (GSH) and causes oxidative stress by generating H 2 O 2 ( Fig. 1A and supplementary Fig. S1B) 14 . As shown in Fig. 1B, BSO-treated cells were found to have a significant decrease of the global hmC level. Further, this observation appeared not to be attributable to an increased apoptosis or necrosis, as BSO treatment did not increase the proportion of apoptotic or necrotic cells as compared to mock-treated cells (supplementary Fig. S1A).
Several studies have shown that BSO-treated SY5Y cells display a significant decrease in proliferation in response to oxidative stress 15,16 . We therefore measured proliferation with the xCELLigence technology, which records changes in conductivity that are proportional to the number of cells attached to the incubation chamber. Upon treatment with BSO, the proliferation decrease was greater in TET1-knockdown cells than in TET2-knockdown cells or in control cells (supplementary Fig. S1D). It is worth noting that TET3 does not seem to be expressed in SY5Y cells, and that the effects observed here cannot be attributed to an off-target effect of TET1-RNAi on TET2 (supplementary Fig. S1D, and data not shown).
To confirm and extend the above findings, we used an in vivo model of oxidative stress: double-knockout mice lacking the genes encoding glutathione peroxidases 1 and 2 (called hereafter GPx1/2 Dko). The GPx enzymes are known to reduce H 2 O 2 and are viewed as the major cellular antioxidant enzymes, especially in intestinal epithelial cells where reactive oxygen species accumulate upon their depletion (Fig. 1A) 10 . Dot-blot and mass spectrometry experiments were performed on GPx1/2 Dko and wild-type (wt) colon epithelia, and in agreement with our data on SY5Y cells, the global hmC level was found to be lower in GPx1/2 Dko mice than in their wt counterparts (Fig. 1C).
Our results thus suggest that hydroxymethylcytosine levels are decreased upon exogenous and in vivo oxidative assaults, and that SY5Y cells with reduced TET1 expression are more sensitive to oxidative stress.
hmC deep-sequencing profiles of BSO-treated SY5Y cells highlight pathways involved in the oxidative stress response. The global decrease in hmC seen upon treatment of SY5Y cells with BSO ( Fig. 1B) led us to interrogate its genome-wide distribution. For this we used the previously described hmC-selective chemical labeling technique hMe-seal to selectively isolate hydroxymethylated DNA fragments 17 , and subjected these to Illumina deep sequencing (referred here as "hmC-seq"). As previously described (reviewed in 18 ), gene bodies appeared most highly represented among the captured fragments, exons being more enriched than introns (supplementary Fig. S2A). In agreement with the above dot-blot and mass spectrometry data, BSO-treated SY5Y cells displayed a significant global decrease in hmC ( Fig. 2A, left panel).
Despite the observed global decrease of hmC, we next looked at differentially hydroxymethylated genes (termed "dhMGs") to see if certain genes might locally loose or even gain hmC upon oxidative stress. We found 2846 dhMGs (supplementary table I), 53% of which displayed a local decrease in hmC and 47% a local increase ( Fig. 2A, right panel). Remarkably, Ingenuity gene ontology analysis applied to these dhMGs revealed significant over-representation of toxicogenomic pathways associated with oxidative stress response, such as mitochondrial dysfunction, decreased polarization of mitochondria, and cytochrome P450 response (Fig. 2B). Of note, the most highly over-represented pathways were different according to whether a differentially hydroxymethylated gene showed a gain or a loss of hmC: the mitochondrial dysfunction pathway in the former case and pathways related to the physiopathology of the heart, liver and kidney in the latter (supplementary Fig. S2B).
Interestingly, genes such as BCL2 or the BCL2-related MCL1 gene, found to be differentially hydroxymethylated ( Fig. 2C; see also sequencing tracks on Fig. 2D), are known to exert important functions during the oxidative stress response: Oxidative stress is attenuated in mice overexpressing BCL2, and the hepatitis B HBx protein sensitizes liver cells to H 2 O 2 -induced oxidative stress by reducing MCL1 expression 19,20 .
Our results thus show that BSO treatment affects the hmC patterns both globally and locally, notably in genes important for a protective response to oxidative stress.  Mice lacking the glutathione peroxidases 1 and 2 display an altered hydroxymethylation on genes involved in the oxidative stress response. As for the above, hmC-seq experiments were performed on GPx1/2 wt and Dko colon epithelia, and confirmed the results obtained with SY5Y cells: mice lacking the GPx genes showed a significant overall decrease in hmC (Fig. 3A, left panel); reads most often corresponding to a gene body location and more frequently to an exon than to an intron location (supplementary Fig. S3A). We next identified 475 dhMGs (supplementary table II), 81% of which showed a local decrease in hmC and 19%, a local increase (Fig. 3A, right panel). Interestingly, a significant proportion of dhMGs were found to be involved in toxicogenomic pathways associated exclusively with oxidative stress, and notably with GSH depletion, precisely emphasizing GPx deletion that occurs in the Dko mice (Fig. 3B,C. In supplementary figure S3B, results are shown separately for genes showing a gain or loss of hmC).
As for BSO-treated SY5Y cells, we again found key genes involved in the oxidative stress response. Examples include Nfe2l1/Nrf1 and Prdx6 (Fig. 3C,D), whose products are known, respectively, to initiate transcription of antioxidant genes and to participate in the redox regulation of the cell by reducing H 2 O 2 21,22 . It is also worth noting that our pathway analyses did not highlight the same pathways in BSO-treated SY5Y cells and GPx1/2 Dko mice (Figs 2B and 3B). This probably reflects the fact that the cultured cells were treated for only 48 h, whereas mice lacking the GPx genes were under constant oxidative assault from embryogenesis to the day of sacrifice (1 month).
In summary, these results confirm restructuration of the hydroxymethylome, here in response to an in vivo oxidative stress, and uncover a set of genes/pathways that reflect the need for an appropriate cellular response.

MicroRNAs encoded by differentially hydroxymethylated genes are predicted to target transcripts involved in the oxidative stress response. During our analysis of dhMGs in SY5Y cells
and GPx1/2 wt and Dko mice, we observed an unexpected high proportion of differentially hydroxymethylated microRNA-encoding sequences (Fig. 4A). MicroRNAs (miRs) are non-coding RNAs of 18 to 24 nucleotides long that hybridize to target mRNAs and, depending on the level of complementarity, cause their degradation or translational repression 23 . As depicted in Fig. 4A, 25% of the dhMGs identified in SY5Y cells and 21% of those identified in GPx1/2 Dko mice were found to correspond to miR-encoding sequences (see also supplementary table III).
As many miR targets are conserved among mammalian species 24 , we examined more closely the miR-encoding sequences showing altered hydroxymethylation in BSO-treated SY5Y cells or in GPx1/2 Dko mice. Of these miRs, 20 showed a robust increase or decrease in hmC in both SY5Y cells and in GPx1/2 Dko mice ( Fig. 4B and supplementary Fig. S4). Focusing on this set of 20 miRs, we used the miRWalk database, a bioinformatic tool that exploits several prediction programs 25 , to build a comprehensive list of predicted mRNAs targets (Fig. 4B, and supplementary table IV). Surprisingly, Ingenuity pathway analysis applied to the corresponding human and mouse mRNAs revealed over-representation of toxicogenomic pathways involved in oxidative stress, as observed earlier ( Fig. 4C; see also Figs 2B and 3B, and supplementary Figs S2B and S3B).
These results thus suggest that hydroxymethylation pattern changes in DNA regions that do not solely code for proteins (here, microRNA genes), but that may contribute to determine the cellular response to oxidative stress.

Discussion
In this report, we have used different approaches to investigate the links between oxidative stress, the global hmC level, and the local distribution of this epigenetic mark. In vitro and in vivo data on SY5Y cells and in our mouse model show that hmC is significantly decreased upon oxidative stress. In addition, we show that TET1 depletion sensitizes SY5Y cells to oxidative stress induced by BSO, the most commonly used drug to increase intracellular levels of reactive oxygen species (ROS). This suggests that TET1, at least in part, may play a role in protecting cells against oxidative stress.
Of note, our genome-wide hmC profiling data shows a global decrease of hmC in both SY5Y cells and in GPx1/2 Dko mice, but locally, some of the differentially hydroxymethylated regions display an increase of hmC in response to oxidative stress. This increase at specific genomic loci might be caused by recruitment of TET enzymes to these regions in order to initiate an appropriate transcriptional program to respond to oxidant assaults.
Surprisingly, around 25% of the dhMGs identified here in either SY5Y or mice lie in microRNA-encoding sequences. In silico analyses have revealed that these miRs could target important genes involved in ROS detoxification. MiRs have been shown to target the transcripts of numerous genes that are associated to life-threatening diseases such as neurodegenerative disorders or cancer, and this set of miRs might be exploited, in the future, in biomarker discovery for oxidative stress related diseases 26,27 . The present work could thus be the starting point in the development of exciting biomedical applications, once we have gained better understanding of the links between oxidative stress and hydroxymethylcytosine patterns, and between miR-genes hydroxymethylation and the miR-genes transcript levels.
Of note, our study is the first that evaluated the impact of oxidative stress in vivo. On the one hand, prolonged oxidative stress is associated with inflammation, which is thought to precede cancer development 28 . In the other hand, a decrease of hmC is now widely accepted as a hallmark of many cancers, linking the TET proteins to cancerogenesis 29 . It is thus tempting to speculate that the decrease of hmC observed during cancer development could be explained, at least partly, by a decrease of hmC emerging during oxidative stress and inflammation.
In conclusion, we propose, on the basis of our results, a model linking oxidative stress to hydroxymethylcytosine patterns. Our results notably highlight the unexpected potential role of certain microRNAs in determining how cells respond to oxidative stress. In addition, we expect that these results, in association A set of 20 dhMG-associated miR genes, common to human and mouse, were identified and their targets predicted. Next, the identified mRNA targets were processed with the Ingenuity software. (C) IPA-Tox applied to the predicted S5Y5 and mouse mRNA targets reveals significant over-representation of oxidative-stress-related pathways (marked in brown). In each graph, the x axis represents the log(p-value) and the dashed line, the significance threshold for pathway over-representation. See also supplementary Fig. S4. with other studies, will lift the veil on molecular aspects that could explain the global decrease of hmC in cancers. However, additional investigation on how the TETs are recruited on genes that gain hmC, as well as the transcriptional effects of differential hydroxymethylation on oxidative stress-related genes and microRNAs, still need to be addressed and should be the subject of future studies.

Statistical analysis. Unless otherwise indicated, all experiments included technical replicates, and
were repeated at least three independent times. Data and graphs are presented as averages ± standard deviations. Data were compared by means of two-tailed unpaired t-tests for comparison, and statistical significance was accepted at p-values ≤ 0.05. * and ** represent p-values ≤ 0.05 and ≤ 0.01, respectively. Statistics relative to hmC deep-sequencing analysis are further described in supplementary methods. Ingenuity p-values were calculated online with Fisher's exact t-test.
Cell culture and treatments. SY5Y cells (ATCC #CRL-2266) were maintained in 1:1 Dubelcco's modified Eagle's medium (DMEM) and F12 medium, then supplemented with 10% fetal bovine serum and grown at 37 °C under 5% CO2. The cells were then treated with 500 μ M BSO for 48 h. The mock-treatment control contained water instead of BSO.
GPx1/2-knockout mice. In our in vivo experiments, purified genomic DNA extracts obtained from the colon epithelia of 129 mice with combined disruption of the GPx1 and GPx2 genes were generously provided by Dr. F. F. Chu, and are fully described in 10 . The methods used on mice were carried out in accordance with the approved NIH guidelines, which authorize to perform all of the animal studies described in the related manuscript. In addition, all experimental protocols were approved by the COH Research Animal Care Committee (Duarte, California).
Five wild-type and five double-knockout (Dko) mice were used for hmC quantification in dot-blot and mass spectrometry experiments. Of these, two wild-type and two Dko mice were used for subsequent deep-sequencing analysis.
Dot blot for 5-hydroxymethylcytosine quantification. Dot-blot experiments were done as previously described, with some modifications 30 . For complete procedure, please refer to supplementary methods.
Analysis of global hmC levels by mass spectrometry (LC/MS/MS). 500 ng of genomic DNA was incubated with 5 U of DNA Degradase Plus (Zymo Research) at 37 °C for 3 h. The resulting mixture of 2'-deoxynucleosides was analysed on a Triple Quad 6500 mass spectrometer (AB Sciex) fitted with an Infinity 1290 LC system (Agilent) and an Acquity UPLC HSS T3 column (Waters), using a gradient of water and acetonitrile with 0.1% formic acid. External calibration was performed using synthetic standards, and for accurate quantification, all samples and standards were spiked with isotopically labeled nucleosides. HmC levels are expressed as a percentage of total cytosines. Hydroxymethylated DNA fragment affinity purification (hMe-seal). At least 500 ng of genomic DNA was diluted in ultra-pure water to 35 ng/μ L and then sonicated in cold water with a Bioruptor sonicator (Diagenode) to obtain fragments averaging 300 bp in size. The fragmented DNA was used in combination with the hydroxymethyl collector (Active Motif) following to the manufacturer's protocol. Briefly, a glucose moiety that contains a reactive azide group was enzymatically linked to hydroxymethylcytosine in DNA, creating glucosyl-hydroxymethylcytosine. Next, a biotin conjugate was chemically attached to the modified glucose via a "click reaction", and magnetic streptavidin beads were used to capture the biotinylated-hmC DNA fragments. After extensive washing steps and chemical elution, the hydroxymethylated DNA fragments released from the beads were used in sequencing experiments. Ingenuity software. Ingenuity IPA software was used for toxicogenomic analysis ("IPA-TOX"). The genes and hmC fold changes were loaded into the Ingenuity database and then core analyses were done, using default parameters and specifying the relevant species ("Human" for SY5Y cells and "Mouse" for mice samples). The tissues used in the experiments were also indicated.
For further details on experimental procedures, including, RNA interference and RT-qPCR analysis, dot blot for 5-hydroxymethylcytosine quantification, Reactive Oxygen Species measurements, measures of apoptosis and necrosis, proliferation experiments, library preparation, deep sequencing workflow, and statistical analysis, please refer to supplementary methods.