ABSTRACT
Epigenetic mechanisms determine the access of regulatory factors to DNA during events such as transcription and the DNA damage response. However, the global response of histone modifications and chromatin accessibility to UV exposure remains poorly understood. Here, we report that UV exposure results in a genome-wide reduction in chromatin accessibility, while the distribution of the active regulatory mark H3K27ac undergoes massive reorganization. Genomic loci subjected to epigenetic reprogramming upon UV exposure represent target sites for sequence-specific transcription factors. Most of these are distal regulatory regions, highlighting their importance in the cellular response to UV exposure. Furthermore, UV exposure results in an extensive reorganization of super-enhancers, accompanied by expression changes of associated genes, which may in part contribute to the stress response. Taken together, our study provides the first comprehensive resource for genome-wide chromatin changes upon UV irradiation in relation to gene expression and elucidates new aspects of this relationship.
INTRODUCTION
Maintenance of genome integrity is essential for cell survival and reproduction as DNA is continuously challenged, for example, by environmental factors like solar UV light. UV light primarily induces DNA lesions, such as cyclobutane pyrimidine dimers (CPDs), pyrimidine 6-4 pyrimidone photoproducts (6-4PPs) and their Dewar isomers, but it can also cause DNA double-strand breaks (Rastogi et al., 2010). Many DNA repair pathways have evolved to correct these diverse types of damages in order to prevent DNA lesions causing pathologies such as skin cancer (Iyama and Wilson, 2013; Lo and Fisher, 2014). The DNA damage response involves recognition of damaged sites and transduction of the signal to effector molecules, which further regulate transcriptional changes, DNA repair, cell cycle arrest or, if the damage is too severe, cell death (Zhou and Elledge, 2000).
Consequences of DNA lesions, either as a result of UV irradiation or other genotoxic agents, to cell physiology are widely governed by changes in gene expression regulated on various levels (Andrade-Lima et al., 2015; Bowden et al., 2015; Dawes et al., 2014; Sesto et al., 2002). UVB has been shown to inhibit initiation of transcription by impeding the binding of RNA polymerase II to the promoters of many transcribed genes (Gyenis et al., 2014). UV can further modulate the action of polymerases through activation of transcription factors such as the tumor suppressor protein TP53 and activator protein 1 (AP1) resulting in gene expression changes (Engelberg et al., 1994; Laptenko and Prives, 2006). Recent studies have revealed that epigenetic mechanisms also play a role in the transcriptional response following UV irradiation. For example, the components of the ATP-dependent SWI–SNF chromatin remodeling complex BRG1 and BRM (also known as SMARCA4 and SMARCA2, respectively) have been implicated in the transcriptional regulation of UV-induced genes (Hassan et al., 2014; Zhang et al., 2014). Furthermore, many posttranslational histone modifications have been shown to be modulated during the cellular response to DNA damage (Corpet and Almouzni, 2009; Kumar et al., 2012; Li et al., 2014; Tjeertes et al., 2009). For example, following UV exposure, the histone acetyltransferase p300 (also known as EP300) is recruited close to the MMP-1 promoter, regulating MMP-1 transcription by its catalytic activity (Kim et al., 2009). UV exposure induces phosphorylation of histone H3 serine 28 (H3S28p) at promoters of stress-response genes and causes dissociation of histone deacetylase (HDAC) co-repressor complexes that further accompanies enhanced levels of histone acetylation and transcriptional induction of these genes (Keum et al., 2013; Sawicka et al., 2014). H3S28 phosphorylation is furthermore involved in the regulation of RNA-polymerase-III-dependent transcription (Zhang et al., 2011). Given that many long non-coding RNAs (lncRNAs) are deregulated upon X-ray irradiation and doxorubicin treatment in mammalian fibroblasts (Rashi-Elkeles et al., 2014; Younger et al., 2015), and several lncRNAs have been implicated in the DNA damage response (Huarte et al., 2010; Liu and Lu, 2012; Negishi et al., 2014; Wan et al., 2013), it is likely that lncRNAs also play a crucial role in modulating gene expression upon UV irradiation.
In addition to promoters, the importance of distal regulatory regions in regulating gene expression is increasingly being appreciated. Enhancers are among such distal regulatory regions that function to augment the transcription of associated genes (Shlyueva et al., 2014). Their active state is characterized by high levels of histone 3 acetylated at lysine 27 (H3K27ac) and chromatin accessibility (Creyghton et al., 2010; Rada-Iglesias et al., 2011; Wang et al., 2008b). Recent findings suggest that enhancer regions might also play an essential role in the DNA damage response by recruiting factors such as TP53 (Younger et al., 2015). Recently, clusters of putative enhancers (named ‘super-enhancers’) have been shown to regulate the expression of genes defining cell identity (Whyte et al., 2013). They have also been found to be deregulated in diseases such as cancer and inflammation (Brown et al., 2014; Chapuy et al., 2013; Hnisz et al., 2013; Mansour et al., 2014; Schmidt et al., 2015; Vahedi et al., 2015), but have not yet been studied in context of UV exposure.
Chromatin dynamics constitute a crucial part of the DNA damage response and many studies have shown recruitment of several ATP-dependent chromatin-remodeling complexes as well as histone-modifying enzymes to damage sites (Corpet and Almouzni, 2009; Lans et al., 2012; Price and D'Andrea, 2013). An early and transient nucleosome destabilization at the site of damage has been shown (Rubbi and Milner, 2003; Smerdon and Lieberman, 1978, 1980; Smerdon et al., 1982), which allows binding of the repair machinery to DNA lesions before chromatin architecture is restored after the repair (the ‘prime-repair-restore model’) (Soria et al., 2012). In this direction, a recent study has shown that, upon UV irradiation, chromatin transiently undergoes de-condensation to allow binding of several repressive factors such as the polycomb complex component BMI1, the nucleosome-remodeling deacetylase complex (NuRD) as well as the heterochromatin protein HP1 and the HP1-binding protein 3 (HP1BP3) (Izhar et al., 2015). Another study proposes that during the response to UV irradiation, recruitment of repressive complexes induces condensation of chromatin to protect DNA from further potential damage (Burgess et al., 2014). Until now, the function and the resultant chromatin changes following recruitment of these complexes remains controversial (Papamichos-Chronakis and Peterson, 2013). Therefore, there is a need to perform a comprehensive investigation into the dynamics of chromatin state following UV irradiation.
Here, we have generated genome-wide datasets to study transcriptome and epigenome changes in response to UV irradiation. We assessed genome-wide changes in the transcriptome by performing RNA-seq, chromatin accessibility by performing formaldehyde-assisted isolation of regulatory elements (FAIRE)-seq and epigenetic landscape by performing H3K27ac chromatin immunoprecipitation (ChIP)-seq at 6 h after UV irradiation of murine fibroblast cells (NIH3T3). Strikingly, we observed a genome-wide loss of chromatin accessibility affecting all genomic regions. The active histone mark H3K27ac also showed massive reorganization throughout the genome. Many potential regulatory regions harboring open chromatin and H3K27ac undergo changes upon UV exposure and serve as binding sites for sequence-specific transcription factors. Importantly, we also observed a dramatic remodeling of super-enhancers that often accompanied expression changes of associated genes. In general, a large fraction of the UV-induced gene expression changes could be explained by the observed chromatin changes. Our observations also reveal that the chromatin status prior to UV damage might influence the expression and epigenetic changes following UV exposure. Overall, these findings provide the first comprehensive resource revealing global changes in the epigenetic state as well as chromatin accessibility following UV damage and their relation to the UV-induced expression changes.
RESULTS
UV-induced gene expression changes primarily occur at expressed genes exhibiting active chromatin
To identify global transcriptome changes in response to UV, we irradiated murine fibroblast cells (NIH3T3) with UVC and performed genome-wide expression profiling (RNA-seq) on RNA collected after 6 h (Fig. S1A,B). Computational analysis revealed 832 and 1236 genes that were significantly up- or down-regulated in response to UV, respectively (Fig. 1A). Interestingly, both the up- and down-regulated genes were expressed in the untreated condition indicating that active genes are prime responders to UV irradiation (Fig. 1A, inset, Fig. S1C). Importantly, not only the count of down-regulated genes, but also the magnitude of down-regulation was significantly greater compared to that of the up-regulation (Fig. 1B). Gene ontology (GO) enrichment analysis of the down-regulated genes showed a strong enrichment for terms such as regulation of transcription (including genes encoding general transcription factors, mediators and zinc finger proteins) and chromatin organization (Fig. 1C; Table S1). In contrast, the up-regulated genes were involved, for example, in oxidation reduction, cell death and the stress response (Fig. 1C; Table S1). Furthermore, genes in pathways related to translation were found to be up-regulated as exemplified by the up-regulation of ribosomal proteins (e.g. Rpl18 and Rpl17), translation initiation factors (e.g. Eif5b and Eif3g) and elongation factors (e.g. Eef1g and Eef1b2). To further consolidate our findings we also performed Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) and cellular component enrichment analysis for up- and down-regulated genes. These results supported the indication of an enhanced translation process in response to UV irradiation (Fig. S1D,E; Table S1). Although we did not observe any increase in the total protein amount 6 h after UV (Fig. S1F), the potential increase in translation of certain proteins might until that time be employed to compensate for the observed transcriptional down-regulation. A number of genes from the AP1 complex (e.g. Atf3, Fosb and Junb), which is involved in cellular stress response, were found to be strongly up-regulated in response to UV irradiation, whereas, for example, Cbx2, a component of the polycomb multiprotein complex involved in chromatin organization, was down-regulated (Fig. 1D) (Chaum et al., 2009; Darwiche et al., 2005; Tanos et al., 2005).
Active genes are prime responders to UV exposure. (A) MA plot showing expression changes of protein-coding genes 6 h after UV treatment of NIH3T3 cells revealed by RNA-seq analysis (n=3). The inset shows the number of differentially expressed genes (y-axis) in four bins of log2 expression (x-axis). (B) Box plots showing distribution of the absolute values of the log2 expression ratio of differentially expressed genes. The box represents the 25–75th percentiles, and the median is indicated. The whiskers show 1.5 times the interquartile range (IQR) added to the 75th percentile (upper whisker) or subtracted from the 25th percentile (lower whisker). The notches represent median±1.57×IQR/(n0.5). **P<2.2×10−6 (Wilcoxon rank-sum test). (C) Bar plot showing enrichment of biological processes for up-regulated (right panel) and down-regulated (left panel) genes. Bar length represents the number of genes (primary x-axis), and the P-value is shown by line graph (upper x-axis). (D) Independent validation of detected gene expression changes 6 h after UV treatment in NIH3T3 cells by RT-qPCR analysis. Mean±s.e.m. expression levels detected by RNA sequencing are plotted as normalized read counts (left y-axis; n=3) and mean±s.e.m. mRNA abundance measured by RT-qPCR are plotted normalized to Ctcf (ΔCT) (right y-axis; n=4). ***P<0.001; ****P<0.0001 (unpaired t-test for RT-qPCR results and as determined by DESeq for the RNA-seq data). (E) Scatterplots showing gene expression changes for genes harboring different epigenetic features at their promoters in the untreated condition with highlighted differential expressed genes identified by DESeq (blue, up-regulated; red, down-regulated). Genes whose promoters show no enrichment for FAIRE or H3K27ac (double negative, first panel), only enrichment for H3K27ac (H3K27ac positive, second panel), only enrichment for FAIRE (FAIRE positive, third panel) and enrichment for FAIRE and H3K27ac (double positive, fourth panel) are displayed. The table provides the total number of genes and the percentage of differentially expressed genes in each category.
Active genes are prime responders to UV exposure. (A) MA plot showing expression changes of protein-coding genes 6 h after UV treatment of NIH3T3 cells revealed by RNA-seq analysis (n=3). The inset shows the number of differentially expressed genes (y-axis) in four bins of log2 expression (x-axis). (B) Box plots showing distribution of the absolute values of the log2 expression ratio of differentially expressed genes. The box represents the 25–75th percentiles, and the median is indicated. The whiskers show 1.5 times the interquartile range (IQR) added to the 75th percentile (upper whisker) or subtracted from the 25th percentile (lower whisker). The notches represent median±1.57×IQR/(n0.5). **P<2.2×10−6 (Wilcoxon rank-sum test). (C) Bar plot showing enrichment of biological processes for up-regulated (right panel) and down-regulated (left panel) genes. Bar length represents the number of genes (primary x-axis), and the P-value is shown by line graph (upper x-axis). (D) Independent validation of detected gene expression changes 6 h after UV treatment in NIH3T3 cells by RT-qPCR analysis. Mean±s.e.m. expression levels detected by RNA sequencing are plotted as normalized read counts (left y-axis; n=3) and mean±s.e.m. mRNA abundance measured by RT-qPCR are plotted normalized to Ctcf (ΔCT) (right y-axis; n=4). ***P<0.001; ****P<0.0001 (unpaired t-test for RT-qPCR results and as determined by DESeq for the RNA-seq data). (E) Scatterplots showing gene expression changes for genes harboring different epigenetic features at their promoters in the untreated condition with highlighted differential expressed genes identified by DESeq (blue, up-regulated; red, down-regulated). Genes whose promoters show no enrichment for FAIRE or H3K27ac (double negative, first panel), only enrichment for H3K27ac (H3K27ac positive, second panel), only enrichment for FAIRE (FAIRE positive, third panel) and enrichment for FAIRE and H3K27ac (double positive, fourth panel) are displayed. The table provides the total number of genes and the percentage of differentially expressed genes in each category.
Although long intergenic non-coding RNAs (lincRNA) are crucial regulators of gene expression, to our knowledge there have not been any genome-wide studies exploring the role of lincRNAs in response to UV. Therefore, we decided to explore differentially expressed lincRNAs from our RNA-seq data. Surprisingly, in contrast to the coding genes, differentially expressed lincRNAs were predominantly up-regulated (n=86) and only a small proportion showed down-regulation (n=12) (Fig. S1G). Interestingly, unlike coding genes, up-regulated lincRNAs were mostly either not expressed or only expressed at low levels in the untreated condition. Previous studies have linked lincRNAs to the transcriptional regulation of nearby genes (Kornienko et al., 2013; Luo et al., 2013; Vance and Ponting, 2014). Correlating the expression changes of lincRNAs to the expression change of the closest gene revealed that genes near to down-regulated lincRNAs were almost exclusively reduced in expression, whereas genes next to up-regulated lincRNAs were also often differentially expressed, but were either up- or down-regulated (Fig. S1H). Interestingly, some up-regulated lincRNAs occurred next to established stress response genes such as Btg2, Egr3, Ier2 and Ier3 (Fig. S1I) (Arlt and Schäfer, 2011; Gallitano-Mendel et al., 2008; Imran and Lim, 2013; Landau et al., 2012).
We next investigated whether the chromatin state of gene promoters and distal regulatory elements could predict the transcriptional behavior of genes in response to UV irradiation. H3K27ac is an established mark of active promoters and enhancers. Therefore, we performed ChIP assays in NIH3T3 cells using H3K27ac-specific antibody followed by next-generation sequencing (ChIP-seq) of recovered genomic DNA. To further assess whether chromatin accessibility predicts changes in gene expression, we also performed a FAIRE assay and sequenced the isolated genomic DNA (FAIRE-seq) (Giresi et al., 2007). Computational analysis revealed that the genes that changed expression after UV exposure were mainly expressed genes whose promoters harbored H3K27ac and/or that were accessible [denoted double-positive (promoters enriched for both H3K27ac and FAIRE), or as H3K27ac-positive or FAIRE-positive] compared to genes whose promoters had none of these features (denoted double-negative) (Fig. 1E). In general, as expected, genes that showed neither H3K27ac nor were accessible were mostly not expressed; however, genes in this category who changed their expression upon UV were expressed in the untreated condition (Fig. 1E, left panel). This finding is consistent with our above observation that active genes are primarily differentially expressed upon UV irradiation (Fig. 1A). Taken together, these observations suggest that gene expression status and the prior chromatin state at promoters determine the UV-induced gene expression response.
UV irradiation results in a genome-wide reduction of chromatin accessibility
We next investigated whether UV treatment influences chromatin accessibility. Towards this, we performed FAIRE-seq 6 h after UV treatment in biological duplicates and compared it to FAIRE-seq performed on untreated cells to reveal genome-wide changes in chromatin accessibility following UV irradiation. Quality control analyses revealed a good correlation between replicates and also confirmed that analysis with individual replicates showed the same biological interpretations (Fig. S2A–H).
Analyzing the total accessible genome showed an approximately one-third reduction upon UV irradiation (Fig. 2A). The distribution of accessible sites in the genome revealed a loss in accessibility following UV treatment at all genomic features; however, the reduction was much stronger in non-promoter regions (Fig. 2B; Fig. S2I). Furthermore, enrichments normalized to the genomic feature size revealed that promoters were more enriched for accessible sites compared to other genomic regions (Fig. S2J,K). Our analysis revealed that 16,782 genomic regions lost accessibility following UV treatment (denoted untreated unique), whereas only 2206 acquired open chromatin (denoted UV unique) (Fig. 2C). To further investigate the chromatin compaction, we next analyzed the FAIRE peak width in both conditions, as this may reflect the size of the accessible DNA available for the binding of proteins such as transcription factors. Interestingly, regions that gained accessibility were of relatively smaller width compared to the accessible sites that were lost following UV treatment, with the exception of exonic sites (Fig. 2C,D). We found ∼40-fold more regions losing accessibility compared to regions gaining accessibility, reflecting the finding that chromatin compaction was much more pronounced upon UV exposure (Fig. 2E, left panel). Similar observations were made when analyzing promoters and intergenic regions separately (Fig. 2E, right panels). In order to validate FAIRE-seq-detected chromatin accessibility changes at promoters as well as intergenic regions, we used independent methods to assess chromatin accessibility changes. First, accessibility changes measured by FAIRE-quantitative PCR (FAIRE-qPCR) at earlier time points after UV irradiation revealed dynamic changes of chromatin accessibility over time arguing against any technical bias due to UV-induced crosslinking or that DNA damage is preventing effective isolation of open regions (Fig. S2L). Furthermore, we also employed an assay for transposase-accessible chromatin (ATAC) on untreated and UV-irradiated cells as an independent measure of chromatin accessibility, revealing that the FAIRE-seq results were largely reproducible (Fig. S2M) (Buenrostro et al., 2013, 2015; Davie et al., 2015). To further augment these results, nucleosome occupancy was investigated by assessing total H3 enrichments at such sites by ChIP-qPCR. As expected, the results showed the opposite dynamics to the accessibility changes at the same loci (Fig. 2F; Fig. S2L). This argues that indeed the observed global loss of accessibility upon UV treatment led to chromatin condensation at these loci with an increased amount of nucleosomes.
UV irradiation results in a genome-wide chromatin compaction. (A) Bar plot showing total open genome size in Mbp of untreated and UV-treated NIH3T3 cells. Total open genome size was calculated as the sum of the widths of all enriched FAIRE peaks. The percentage displays the loss of accessible chromatin upon UV. (B) Same as in A, but for individual genome features. (C) Box plots showing distribution of peak width for unique open regions in untreated and UV-treated cells. The box represents the 25–75th percentiles, and the median is indicated. The whiskers show 1.5 times the interquartile range (IQR) added to the 75th percentile (upper whisker) or subtracted from the 25th percentile (lower whisker). The notches represent median±1.57×IQR/(n0.5). The numbers above each box plot display the counts of unique peaks for each condition. **P<2.2×10−16 (between peak width in untreated and UV condition; Wilcoxon rank-sum test with a cut-off of 0.01). (D) Same as in C, but for different genomic regions. Significant changes upon UV are indicated above the blue box plots. **P<0.01 (promoters, 1.581×10−9; exons, 0.774; introns, 1.762×10−15; intergenic regions, 8.188×10−7). (E) Scatter plot showing the dynamics of FAIRE enrichment changes between untreated and UV conditions for all genomic regions together (left panel), promoter regions (right, upper panel) and intergenic regions (right, lower panel). Peaks with at least 1.5-fold enrichment change following UV irradiation are highlighted (gain in blue, loss in red). The numbers of up- and down-regulated peaks are displayed in the plots. (F) H3 ChIP-qPCR results for promoter (left) and intergenic (right) regions that showed gain or loss in the FAIRE-seq analysis (indicated below). The results are plotted as immunoprecipitated DNA/input DNA (IP/input) (n=4, mean±s.e.m.). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test). (G) Scatter plot showing relationship between chromatin accessibility at promoters and expression of the corresponding genes. Quadrants indicated by dashed lines show changes in expression (y-axis) and accessibility (x-axis) greater than 1.5-fold and numbers indicate the amount of genes falling in each quadrant. (H) Boxplots summarizing the absolute values of expression changes for genes belonging to the green and red categories in G. Significance was determined using a Wilcoxon rank-sum test for non-normal distributions with a cut-off of 0.01. Up and down indicate up- or down-regulation of the gene expression 6 h following UV exposure, FAIRE− indicates loss of FAIRE-seq enrichment 6 h after UV treatment. (I) Genome browser plots for genes showing either concomitant increase or decrease in expression and FAIRE enrichment at their promoters. Each track is shown from 0 to the indicated height. (J) Validation of RNA-seq and FAIRE-seq results by independent RT-qPCR and FAIRE-qPCR for genes selected from the red and blue categories in G (indicated by arrows). Mean±s.e.m. expression measured by RNA-seq analysis is plotted as normalized read counts on the left y-axis (n=3) and the mean±s.e.m. mRNA abundance determined by RT-qPCR is plotted as ΔCT on the right y-axis, normalized to Ctcf (n=3). As determined by qPCR, FAIRE enrichment above input is plotted normalized to FAIRE enrichment at the Ctcf promoter at the right y-axis (n=4, error bars represent s.e.m.), whereas FAIRE-seq enrichment normalized to input is plotted at the left y-axis (n=2). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test for RT-qPCR, FAIRE-qPCR and FAIRE-seq, and as determined by DESeq for RNA-seq data). (K) Relative expression determined by RT-qPCR is plotted as the fold change for NIH3T3 cells 30 min and 6 h after UV treatment relative to the untreated condition (ΔΔCT normalized to Ctcf, n=3, error bars represent s.e.m.). FAIRE enrichment above input determined by FAIRE-qPCR at same time points was normalized to FAIRE enrichment at the Ctcf promoter and then plotted as fold change to untreated condition (n=4, error bars represent s.e.m.). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test). (L) Same plot as in G, but FAIRE enrichments are shown for peaks falling in intergenic regions and plotted against the expression changes of the closest gene. (M) Same plot as in H, but for the red and green quadrants in L.
UV irradiation results in a genome-wide chromatin compaction. (A) Bar plot showing total open genome size in Mbp of untreated and UV-treated NIH3T3 cells. Total open genome size was calculated as the sum of the widths of all enriched FAIRE peaks. The percentage displays the loss of accessible chromatin upon UV. (B) Same as in A, but for individual genome features. (C) Box plots showing distribution of peak width for unique open regions in untreated and UV-treated cells. The box represents the 25–75th percentiles, and the median is indicated. The whiskers show 1.5 times the interquartile range (IQR) added to the 75th percentile (upper whisker) or subtracted from the 25th percentile (lower whisker). The notches represent median±1.57×IQR/(n0.5). The numbers above each box plot display the counts of unique peaks for each condition. **P<2.2×10−16 (between peak width in untreated and UV condition; Wilcoxon rank-sum test with a cut-off of 0.01). (D) Same as in C, but for different genomic regions. Significant changes upon UV are indicated above the blue box plots. **P<0.01 (promoters, 1.581×10−9; exons, 0.774; introns, 1.762×10−15; intergenic regions, 8.188×10−7). (E) Scatter plot showing the dynamics of FAIRE enrichment changes between untreated and UV conditions for all genomic regions together (left panel), promoter regions (right, upper panel) and intergenic regions (right, lower panel). Peaks with at least 1.5-fold enrichment change following UV irradiation are highlighted (gain in blue, loss in red). The numbers of up- and down-regulated peaks are displayed in the plots. (F) H3 ChIP-qPCR results for promoter (left) and intergenic (right) regions that showed gain or loss in the FAIRE-seq analysis (indicated below). The results are plotted as immunoprecipitated DNA/input DNA (IP/input) (n=4, mean±s.e.m.). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test). (G) Scatter plot showing relationship between chromatin accessibility at promoters and expression of the corresponding genes. Quadrants indicated by dashed lines show changes in expression (y-axis) and accessibility (x-axis) greater than 1.5-fold and numbers indicate the amount of genes falling in each quadrant. (H) Boxplots summarizing the absolute values of expression changes for genes belonging to the green and red categories in G. Significance was determined using a Wilcoxon rank-sum test for non-normal distributions with a cut-off of 0.01. Up and down indicate up- or down-regulation of the gene expression 6 h following UV exposure, FAIRE− indicates loss of FAIRE-seq enrichment 6 h after UV treatment. (I) Genome browser plots for genes showing either concomitant increase or decrease in expression and FAIRE enrichment at their promoters. Each track is shown from 0 to the indicated height. (J) Validation of RNA-seq and FAIRE-seq results by independent RT-qPCR and FAIRE-qPCR for genes selected from the red and blue categories in G (indicated by arrows). Mean±s.e.m. expression measured by RNA-seq analysis is plotted as normalized read counts on the left y-axis (n=3) and the mean±s.e.m. mRNA abundance determined by RT-qPCR is plotted as ΔCT on the right y-axis, normalized to Ctcf (n=3). As determined by qPCR, FAIRE enrichment above input is plotted normalized to FAIRE enrichment at the Ctcf promoter at the right y-axis (n=4, error bars represent s.e.m.), whereas FAIRE-seq enrichment normalized to input is plotted at the left y-axis (n=2). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test for RT-qPCR, FAIRE-qPCR and FAIRE-seq, and as determined by DESeq for RNA-seq data). (K) Relative expression determined by RT-qPCR is plotted as the fold change for NIH3T3 cells 30 min and 6 h after UV treatment relative to the untreated condition (ΔΔCT normalized to Ctcf, n=3, error bars represent s.e.m.). FAIRE enrichment above input determined by FAIRE-qPCR at same time points was normalized to FAIRE enrichment at the Ctcf promoter and then plotted as fold change to untreated condition (n=4, error bars represent s.e.m.). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test). (L) Same plot as in G, but FAIRE enrichments are shown for peaks falling in intergenic regions and plotted against the expression changes of the closest gene. (M) Same plot as in H, but for the red and green quadrants in L.
Intrigued by the largely reproducible chromatin accessibility response upon UV treatment, we next asked whether these changes were caused by the DNA-damage response occurring at these loci. Because it has been reported that UV damage induces γH2AX at damaged sites within 6 h in human fibroblasts (Oh et al., 2011) and we also observed a global gain of γH2AX 6 h following UV treatment (Fig. S2N), we monitored damaged DNA regions by performing γH2AX ChIP in untreated and UV-irradiated NIH3T3 cells. We measured γH2AX enrichments at regions that were non-accessible or accessible in the untreated condition and either changed or did not change their accessibility at 6 h after UV exposure. These results indicate that UV-induced γH2AX occurrence is independent of the accessibility status of the regions before UV treatment and, furthermore, chromatin accessibility changes seem to be independent from the γH2AX occurrence (Fig. S2O,P). We also investigated whether the sites responding or non-responding in accessibility upon UV irradiation had different chromatin features prior to the treatment by using published ChIP-seq datasets for different histone modifications, histone variants and phosphorylated RNA polymerase generated from ChIP experiments performed on untreated NIH3T3 cells. The analyses revealed that certain chromatin features are significantly different between these classes indicating that the chromatin status prior to UV irradiation impacts on the chromatin response to UV with respect to changes in accessibility (Fig. S2Q).
The observation of a major chromatin accessibility loss is further consistent with our findings that UV-induced transcriptional down-regulation is much more pronounced compared to up-regulation (Fig. 1B). Prompted by these findings, we next investigated the correlation of chromatin accessibility at promoters with the expression level of the genes in the untreated and UV conditions and observed a good positive correlation (Fig. S2R). We then asked whether changes in the accessibility of promoters and distal regions are in accordance with alterations in expression of the associated genes following UV treatment. This analysis showed that the majority of genes whose promoters lost accessibility were also down-regulated (n=513, red) (Fig. 2G–I). These genes were predominantly those involved in regulation of transcription and chromatin organization (Fig. S2S; Table S1). We observed another set of genes that lost accessibility at promoters but slightly gained transcription (n=295, green) (Fig. 2G,H). There were four genes found to gain accessibility as well as expression (Cdkn1a, Fosl1, Marcksl1 and Olfr1226). Among them were typical DNA-damage response genes like the proliferation inhibitor Cdkn1a and the AP1 component Fosl1, whose changes in expression and promoter accessibility were validated by independent reverse transcription quantitative real-time PCR (RT-qPCR) and FAIRE-qPCR (Fig. 2G,J). Moreover, analysis of expression and accessibility changes that had occurred by 30 min after UV irradiation revealed that accessibility changes either precede expression changes or occur at the same time (Fig. 2K). A comparison of accessibility changes at potential distal regulatory regions with the expression changes of the closest gene revealed patterns comparable to those seen for promoters. Most genes whose distal regions lost accessibility were down-regulated (n=339, red) following UV treatment (Fig. 2L,M). Another set lost accessibility at intergenic regions but showed a minor gain in expression (n=212, green) (Fig. 2L,M).
Overall, these findings suggest that UV exposure causes a genome-wide loss of chromatin accessibility that largely accompanies transcriptional down-regulation of the associated genes. However, despite a genome-wide reduction in transcriptional competence, distinct genomic loci associated with the stress response are selectively kept accessible to allow their expression.
UV exposure causes a global reprogramming of the H3K27ac mark
Following our observations of a global change in chromatin accessibility, we next investigated the effect of UV on a specific histone modification that defines distinct transcriptional states. To this end, we selected H3K27ac, an established marker of active promoters and enhancers, and performed ChIP-seq for H3K27ac in untreated and UV-treated NIH3T3 cells. Quality control analyses with respect to antibody specificity, chromatin shearing, alignment statistics and replicate correlation were performed (Fig. S3A–I). Whereas the total fraction of the genome and of various genomic regions enriched for H3K27ac stayed largely unchanged following UV treatment (Fig. 3A,B; Fig. S3J), the enrichment for H3K27ac at these genomic loci underwent a dramatic reorganization. Sites that gained H3K27ac (n=12,241) showed a significantly higher peak width compared to regions that lost H3K27ac following UV irradiation (n=11,441) (Fig. 3C). Furthermore, this observation pertains also to regions distributed throughout various genomic locations (Fig. S3K). The enrichment for H3K27ac at different genomic regions mainly showed an increase in H3K27ac enrichment after UV (Fig. 3D; Fig. S3L). We observed massive changes with respect to H3K27ac enrichment after UV treatment with nearly equal numbers of up- and down-regulated sites in response to UV irradiation (22,869 up and 20,103 down) (Fig. 3E, left panel). This observation further holds true when looking at the enrichment level changes either at promoters or at peaks falling in intergenic regions (Fig. 3E, right panels).
UV exposure causes massive remodeling of H3K27ac landscape. (A) Bar plot showing total genome size in Mbp for the genome occupied by H3K27ac in untreated and UV-treated NIH3T3 cells calculated as the sum of the widths of all enriched H3K27ac peaks. The percentage displays the change of total genome size occupied by H3K27ac between UV-irradiated and untreated cells. (B) Same as in A, but for different genomic regions. (C) Box plot showing distribution of peak width for unique H3K27ac regions in untreated and UV-treated cells. The numbers of unique peaks for each condition are indicated on top. The box represents the 25–75th percentiles, and the median is indicated. The whiskers show 1.5 times the interquartile range (IQR) added to the 75th percentile (upper whisker) or subtracted from the 25th percentile (lower whisker). The notches represent median±1.57×IQR/(n0.5). **P<2.2×10−16 (Wilcoxon rank-sum test for non-normal distributions). (D) Distribution of H3K27ac peak enrichment in four genomic regions for untreated and UV-treated cells at unique sites (i.e. H3K27ac is at least 1.5-fold enriched above background either in the untreated or UV condition). Numbers above box plots indicate counts of unique peaks for this condition. **P<0.01 (Wilcoxon rank-sum tests for non-normal distribution; promoters, 0.7776; exons, 1.281×10−6; introns, <2.2×10−16; intergenic regions, <2.2×10−16). (E) Scatter plot showing dynamics of H3K27ac enrichment changes between untreated and UV conditions for all genomic regions together (left panel), for promoter regions (right, upper panel) and for intergenic regions (right, lower panel). Peaks with at least 1.5-fold enrichment change upon UV treatment are highlighted (gain in blue, loss in red). The numbers of peaks falling in each of these two categories are displayed on the plots. (F) Scatter plot showing relationship between H3K27ac enrichment changes at promoters (x-axis) and gene expression changes (y-axis). Quadrants indicated by dashed lines show changes in expression and H3K27ac enrichment greater than 1.5-fold, and numbers indicate the counts of genes falling in each category. (G) Boxplots summarizing the absolute values of expression changes for genes highlighted with the respective color in the four categories in F. Up/down indicates up-/ down-regulation of the gene expression, H3K27ac−/+ indicates H3K27ac enrichment loss/ gain 6 h following UV treatment. Significance was determined using a Wilcoxon rank-sum test for non-normal distributions with a cut-off of 0.01. **P<0.01 (blue versus green, 6.477×10−10; blue versus red, 4.718×10−12; green versus orange, 1.261×10−11; green versus red, <2.2×10−16; orange versus red, 1.059×10−13). (H) Genome browser plots of genes showing either concomitant increase or decrease in expression as well as H3K27ac enrichment at their promoters. Each track is shown from 0 to the indicated height. (I) Validation of RNA-seq and H3K27ac ChIP-seq results by independent RT-qPCR and ChIP-qPCR for genes selected from the red and blue categories of F (indicated by arrows). Mean±s.e.m. expression from RNA-seq is plotted as normalized read counts (n=3) on the left y-axis and mean±s.e.m. mRNA abundance determined by RT-qPCR is plotted as ΔCT, normalized to Ctcf on the right y-axis (n=3). ChIP-seq enrichments normalized to the input are plotted on the left y-axis (n=2), whereas ChIP-qPCR enrichment is plotted as immunoprecipitated DNA/input DNA (IP/input) on the right y-axis (n=4). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test for RT-qPCR, ChIP-qPCR and ChIP-seq, and as determined by DESeq for RNA-seq data). (J) Relative expression determined by RT-qPCR is plotted as fold change for NIH3T3 cells 30 min and 6 h after UV treatment relative to the untreated condition (ΔCT normalized to Ctcf, n=3, error bars represent s.e.m.). H3K27ac enrichments for the promoter of these genes were determined at the same stages by ChIP-qPCR. These results are normalized to input and then plotted as fold change to untreated condition (n=3, error bars represent s.e.m.). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test). (K) Same as in F, but H3K27ac enrichments at intergenic sites are plotted against the expression changes of the closest gene. (L) Boxplots summarize the absolute values of changes in expression associated with each of the four categories highlighted in K (similar to G). Significance was determined using a Wilcoxon rank-sum test for non-normal distributions with a cut-off of 0.01. *P<0.05; **P<0.01 (blue versus green, 1.4×10−4; blue versus red, 3.6×10−2; green versus orange, 2.464×10−6; green versus red, 1.942×10−9).
UV exposure causes massive remodeling of H3K27ac landscape. (A) Bar plot showing total genome size in Mbp for the genome occupied by H3K27ac in untreated and UV-treated NIH3T3 cells calculated as the sum of the widths of all enriched H3K27ac peaks. The percentage displays the change of total genome size occupied by H3K27ac between UV-irradiated and untreated cells. (B) Same as in A, but for different genomic regions. (C) Box plot showing distribution of peak width for unique H3K27ac regions in untreated and UV-treated cells. The numbers of unique peaks for each condition are indicated on top. The box represents the 25–75th percentiles, and the median is indicated. The whiskers show 1.5 times the interquartile range (IQR) added to the 75th percentile (upper whisker) or subtracted from the 25th percentile (lower whisker). The notches represent median±1.57×IQR/(n0.5). **P<2.2×10−16 (Wilcoxon rank-sum test for non-normal distributions). (D) Distribution of H3K27ac peak enrichment in four genomic regions for untreated and UV-treated cells at unique sites (i.e. H3K27ac is at least 1.5-fold enriched above background either in the untreated or UV condition). Numbers above box plots indicate counts of unique peaks for this condition. **P<0.01 (Wilcoxon rank-sum tests for non-normal distribution; promoters, 0.7776; exons, 1.281×10−6; introns, <2.2×10−16; intergenic regions, <2.2×10−16). (E) Scatter plot showing dynamics of H3K27ac enrichment changes between untreated and UV conditions for all genomic regions together (left panel), for promoter regions (right, upper panel) and for intergenic regions (right, lower panel). Peaks with at least 1.5-fold enrichment change upon UV treatment are highlighted (gain in blue, loss in red). The numbers of peaks falling in each of these two categories are displayed on the plots. (F) Scatter plot showing relationship between H3K27ac enrichment changes at promoters (x-axis) and gene expression changes (y-axis). Quadrants indicated by dashed lines show changes in expression and H3K27ac enrichment greater than 1.5-fold, and numbers indicate the counts of genes falling in each category. (G) Boxplots summarizing the absolute values of expression changes for genes highlighted with the respective color in the four categories in F. Up/down indicates up-/ down-regulation of the gene expression, H3K27ac−/+ indicates H3K27ac enrichment loss/ gain 6 h following UV treatment. Significance was determined using a Wilcoxon rank-sum test for non-normal distributions with a cut-off of 0.01. **P<0.01 (blue versus green, 6.477×10−10; blue versus red, 4.718×10−12; green versus orange, 1.261×10−11; green versus red, <2.2×10−16; orange versus red, 1.059×10−13). (H) Genome browser plots of genes showing either concomitant increase or decrease in expression as well as H3K27ac enrichment at their promoters. Each track is shown from 0 to the indicated height. (I) Validation of RNA-seq and H3K27ac ChIP-seq results by independent RT-qPCR and ChIP-qPCR for genes selected from the red and blue categories of F (indicated by arrows). Mean±s.e.m. expression from RNA-seq is plotted as normalized read counts (n=3) on the left y-axis and mean±s.e.m. mRNA abundance determined by RT-qPCR is plotted as ΔCT, normalized to Ctcf on the right y-axis (n=3). ChIP-seq enrichments normalized to the input are plotted on the left y-axis (n=2), whereas ChIP-qPCR enrichment is plotted as immunoprecipitated DNA/input DNA (IP/input) on the right y-axis (n=4). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test for RT-qPCR, ChIP-qPCR and ChIP-seq, and as determined by DESeq for RNA-seq data). (J) Relative expression determined by RT-qPCR is plotted as fold change for NIH3T3 cells 30 min and 6 h after UV treatment relative to the untreated condition (ΔCT normalized to Ctcf, n=3, error bars represent s.e.m.). H3K27ac enrichments for the promoter of these genes were determined at the same stages by ChIP-qPCR. These results are normalized to input and then plotted as fold change to untreated condition (n=3, error bars represent s.e.m.). *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001 (unpaired t-test). (K) Same as in F, but H3K27ac enrichments at intergenic sites are plotted against the expression changes of the closest gene. (L) Boxplots summarize the absolute values of changes in expression associated with each of the four categories highlighted in K (similar to G). Significance was determined using a Wilcoxon rank-sum test for non-normal distributions with a cut-off of 0.01. *P<0.05; **P<0.01 (blue versus green, 1.4×10−4; blue versus red, 3.6×10−2; green versus orange, 2.464×10−6; green versus red, 1.942×10−9).
Given that H3K27ac enrichment levels at promoters showed a high correlation to the expression status of a gene in untreated as well as in UV-irradiated cells (Fig. S3M), we next compared UV-induced alterations in H3K27ac levels at promoters to changes in gene expression. This revealed that most promoters losing H3K27ac also showed transcriptional down-regulation (n=390, red) (Fig. 3F–H). These genes were mainly enriched for processes such as regulation of transcription, chromatin organization or regulation of kinase activity (Fig. S3N; Table S1). Furthermore, a large number of genes whose promoters gained H3K27ac were also associated with increased transcription (n=286, blue) (Fig. 3F–H) and they had GO terms such as oxidation reduction, cell death, regulation of transcription and response to DNA damage stimulus (Fig. S3O; Table S1). For the genes that showed gain of H3K27ac but loss of expression and vice versa, the magnitudes of expression changes were less or similar to the other two classes of genes (Fig. 3G). We further validated several of these changes by independent RT-qPCRs and ChIP-qPCRs (Fig. 3F,I). Similar analysis at an earlier time point indicates that changes in H3K27ac enrichment either precede or correlate with the timing of expression changes (Fig. 3J).
Distal regions marked with H3K27ac are known to often represent active enhancers (Creyghton et al., 2010). To investigate whether UV-induced transcriptional changes also involve reprogramming of the enhancer landscape, we performed a comparative analysis of changes in distal H3K27ac sites with changes in expression for the nearest gene following UV treatment. A large fraction of the distal regions that showed loss of H3K27ac was associated with significant transcriptional down-regulation of the nearest gene (n=245, red) (Fig. 3K,L). This class was enriched for genes, for example, involved in regulation of transcription, chromatin organization and phosphorylation (Fig. S3P; Table S1). Similarly, in many distal regions, gain of H3K27ac was associated with an increase in expression of the nearest gene (n=180, blue) (Fig. 3K,L). This fraction contains many genes associated with the cellular response to stress, regulation of cell cycle, circulatory system process and regulation of cell death (Fig. S3Q; Table S1) and includes well-known UV-induced genes like Fosb and Egr1. For a minor fraction of distal regions, loss of H3K27ac was associated with an increase in expression (n=132, green) (Fig. 3K,L). Surprisingly, we also observed a fraction that showed gain of H3K27ac in intergenic regions but down-regulation in the expression of the nearest gene (n=224, yellow) (Fig. 3K,L). These contradicting associations of expression and H3K27ac might be due to a limitation in assigning enhancers to their target genes. This goes along with previous suggestions that enhancer-nearest gene pairing only holds true for ∼40% of genes (Andersson et al., 2014; Doyle et al., 2014; Samee and Sinha, 2014). However, as such contradictory findings were also observed when correlating H3K27ac enrichment changes at promoters to gene expression changes, it is likely that other chromatin features and regulatory mechanisms might be more important in determining the transcriptional state of these genes that we are unable to capture here or that the expression status of these genes is influenced by post-transcriptional mechanisms.
Overall, our observations suggest that UV exposure causes a genome-wide reorganization of the H3K27ac mark at regulatory elements such as promoters and enhancers, which underlie expression changes of crucial genes during the stress response.
Accessible distal enhancers are remodeled following UV treatment
We next investigated the dynamics of H3K27ac-enriched sites that also exhibited accessible chromatin (i.e. double-positive for FAIRE and H3K27ac), assuming that such regions represent sites of intense gene-regulatory activity and harbor a plethora of regulatory elements including transcription-factor-binding sites. The total count of these regions was drastically reduced following UV treatment (Fig. 4A; Fig. S4A). Overlap of untreated and UV double-positive sites revealed that 12,414 of such regions were lost following UV irradiation, whereas only 1912 were gained (Fig. 4B). This substantial loss of double-positive regulatory elements is in line with the loss of chromatin accessibility (Fig. 2A–E). We observed that these sites showed overall loss of accessibility while acquiring H3K27ac following UV treatment (Fig. S4B). The loss of FAIRE enrichment is observed at the uniquely occurring (i.e. enriched only in the untreated sample or the UV-treated sample) as well as the double-positive sites shared in both the untreated and UV-treated cells (‘common’ sites) upon UV treatment (Fig. S4C), whereas the overall gain of H3K27ac enrichment after UV arose from the double-positive sites common between the untreated and UV-treated cells (Fig. S4D). It is interesting to note that although the common sites lost FAIRE enrichment, they displayed a substantial gain in H3K27ac enrichment. Untreated-only double-positive sites mostly showed loss of H3K27ac enrichment following UV exposure, whereas UV-only double-positive sites showed a nearly exclusive gain of H3K27ac enrichment (Fig. S4E,F). Such loss and gain was even more prominent for peaks falling in double-positive intergenic regions (Fig. S4G,H). The untreated-only, double-positive peaks did not gain accessibility, whereas a majority of UV-only double-positive peaks gained chromatin openness (Fig. S4I–L).
Accessible distal enhancers are remodeled following UV exposure. (A) Bar plot showing the number of H3K27ac and FAIRE double-positive sites in untreated and UV-treated cells. The FAIRE peak was considered to be overlapping with the H3K27ac peak if there was an H3K27ac peak within a 500-bp distance from the FAIRE peak summit. (B) Venn diagram showing overlap of double-positive peaks in untreated and UV conditions. (C) Pie charts showing the genomic distribution of peaks that were found in both conditions (common peaks), uniquely in untreated cells (unique untreated peaks) or only in UV-treated cells (unique UV peaks). For comparison, a pie chart is shown displaying the partitioning of the NIH3T3 genome. (D) Bar plot showing enrichment of biological processes for genes near to common peaks. The bars length represent the number of genes (lower x-axis), whereas the P-value is shown as a line graph with respect to the upper x-axis. (E) Same as in D but for genes near to unique untreated peaks. (F). Same as D, but for genes near to unique UV peaks. (G) The table shows the percentage of double-positive regions found in each set of sites (unt., unique untreated peaks; UV, unique UV peaks; NR, not responding, double-positive peaks at least 1.5-fold enriched in untreated and UV samples and not changing more than 1.25-fold after UV exposure) containing the AP1, KLF4, TP53 or CTCF motif (upper panel). The conditions in which the motif was predicted to be enriched by Pscan with respect to the local and global background are highlighted in gray. The density plot shows the average distribution of the motif position relative to the center of FAIRE peaks at the double-positive sites (lower panel). (H) ChIP-qPCR for TP53 in untreated and NIH3T3 cells 6 h after UV treatment for non-target genes as well as potential target genes identified in the UV unique double-positive class. The results are plotted as immunoprecipitated DNA/input DNA (IP/input) (n=3; error bars represent s.e.m.). (I) ChIP-qPCR for CTCF in untreated and NIH3T3 cells 6 h after UV exposure for non-target genes as well as potential target genes identified in the common double-positive class (NR). The results are plotted as IP/input (n=3; error bars represent s.e.m.).
Accessible distal enhancers are remodeled following UV exposure. (A) Bar plot showing the number of H3K27ac and FAIRE double-positive sites in untreated and UV-treated cells. The FAIRE peak was considered to be overlapping with the H3K27ac peak if there was an H3K27ac peak within a 500-bp distance from the FAIRE peak summit. (B) Venn diagram showing overlap of double-positive peaks in untreated and UV conditions. (C) Pie charts showing the genomic distribution of peaks that were found in both conditions (common peaks), uniquely in untreated cells (unique untreated peaks) or only in UV-treated cells (unique UV peaks). For comparison, a pie chart is shown displaying the partitioning of the NIH3T3 genome. (D) Bar plot showing enrichment of biological processes for genes near to common peaks. The bars length represent the number of genes (lower x-axis), whereas the P-value is shown as a line graph with respect to the upper x-axis. (E) Same as in D but for genes near to unique untreated peaks. (F). Same as D, but for genes near to unique UV peaks. (G) The table shows the percentage of double-positive regions found in each set of sites (unt., unique untreated peaks; UV, unique UV peaks; NR, not responding, double-positive peaks at least 1.5-fold enriched in untreated and UV samples and not changing more than 1.25-fold after UV exposure) containing the AP1, KLF4, TP53 or CTCF motif (upper panel). The conditions in which the motif was predicted to be enriched by Pscan with respect to the local and global background are highlighted in gray. The density plot shows the average distribution of the motif position relative to the center of FAIRE peaks at the double-positive sites (lower panel). (H) ChIP-qPCR for TP53 in untreated and NIH3T3 cells 6 h after UV treatment for non-target genes as well as potential target genes identified in the UV unique double-positive class. The results are plotted as immunoprecipitated DNA/input DNA (IP/input) (n=3; error bars represent s.e.m.). (I) ChIP-qPCR for CTCF in untreated and NIH3T3 cells 6 h after UV exposure for non-target genes as well as potential target genes identified in the common double-positive class (NR). The results are plotted as IP/input (n=3; error bars represent s.e.m.).
Analysis of the genomic distribution of double-positive sites revealed, interestingly, that most of the common double-positive sites occurred at promoters, whereas the unique sites were mainly enriched at introns or intergenic regions (Fig. 4C). Inspection of the genes nearest to common, untreated-only and UV-only peaks revealed a substantial number of genes that uniquely gained or lost double-positive sites (Fig. S4M). Genes next to common peaks were mainly enriched in GO terms for housekeeping functions like RNA processing, protein localization, cell cycle regulation, chromatin organization, ribosome biogenesis, regulation of transcription and DNA repair (Fig. 4D; Table S1). Genes harboring double-positive sites only in the untreated condition include those involved in transport of metabolites and chromatin organization (Fig. 4E; Fig. S4N, Table S1). Genes close to double-positive UV-only peaks were enriched for GO term categories like vitamin metabolic processes, metal ion or protein transport, phosphorylation and negative regulation of transcription (Fig. 4F; Table S1). These findings are in line with the results of our transcriptome analysis indicating that cells in general reduce transcription and attempt to utilize available resources following UV exposure.
We next investigated whether untreated-only and UV-only double-positive regions were enriched for transcription-factor-binding motifs. Together with these two sets, we also used a control set of double-positive regions that are enriched in H3K27ac and accessibility in untreated and UV-irradiated cells. Motif analysis led to two levels of motif enrichment information: motifs enriched with respect to the local background (local enrichment, mainly reflecting motifs enriched near peak centers) and motifs enriched with respect to the global background (global enrichment, representing motifs found in surrounding regions). Interestingly, the motif enrichment analysis revealed sites for the same set of factors (FOSL2, FOS, JUNB; all AP1 components), which were found to be up-regulated in their expression upon UV treatment, upstream of the peak center in all three investigated sets of peaks (Fig. 4G; Table S2). However, there were small positional differences with respect to the peak center of the double-positive sites (Table S2). Regions downstream of the untreated-only peaks contained motifs for factors such as TFAP2A, KLF4, SP2 and EGR1 (Table S2). These factors are related to histone deacetylation (SP2) (Phan et al., 2004), repression of TP53 (KLF4) (Rowland et al., 2005), cell growth, apoptosis and DNA damage (SP1) and repression of apoptosis (EGR1) (Huang et al., 1998; de Belle et al., 1999; Zhang et al., 2001). The regions downstream of the UV-unique motifs showed a high enrichment for the TP53 motif (Table S2). TP53 has been shown to mediate UV-induced global chromatin changes (Rubbi and Milner, 2003). These UV-unique double-positive regions harboring TP53 motifs were further validated using ChIP-qPCR (Fig. 4H). Interestingly, these results showed that these sites might already be targeted by TP53 in untreated conditions. These findings are in line with a recent genome-wide analysis of TP53 binding, revealing TP53 to be found predominantly within enhancers and showing that these sites are already poised in the untreated conditions by binding of TP53, which then becomes functional upon phosphorylation in response to DNA damage signals (Younger et al., 2015). The control set of non-changing peaks showed high local enrichment for CTCF and such CTCF binding at target sites in the non-changing class was validated by ChIP-qPCR (Fig. 4I). Interestingly, these CTCF target sites did not show any changes in CTCF enrichment following UV, suggesting that binding of CTCF at these sites might protect them from UV-induced chromatin changes. Taken together, these results suggest that there is epigenetic remodeling of regulatory elements upon UV exposure, possibly involving the action of sequence-specific transcription factors, to either mediate or prevent changes following UV treatment.
UV exposure results in a dramatic reorganization of super-enhancers
Prompted by our observations that UV treatment results in genome-wide epigenetic reprogramming of regulatory elements, we next investigated whether the UV response also involves the reorganization of H3K27ac clusters called super-enhancers. Using our H3K27ac ChIP-seq data, we first identified such clusters of enhancers and then, using a publically available ChIP-seq dataset generated from ChIP perfomed on untreated NIH3T3 cells for histone 3 mono-methylated at lysine 4 (H3K4me1), confirmed their enrichment for H3K4me1 (Fig. S4O–Q) (Zhu et al., 2012). Interestingly, whereas the total number of super-enhancers hardly changed, UV treatment led to dramatic reorganization of these elements in which 360 genes lost and 337 genes gained super-enhancers (Fig. 5A,B). Furthermore, the gained super-enhancers were located closer to transcription start sites than those that were lost following UV treatment (Fig. 5C). Super-enhancers were found more frequently in introns and less frequently in intergenic regions after UV treatment (Fig. 5D,E). GO analysis of the nearest genes to untreated-only super-enhancers revealed that these genes were mainly related to cell morphogenesis, cell motion and cell proliferation (Fig. 5F; Table S1). Super-enhancers occurring uniquely in UV-treated cells were in vicinity of genes involved in cell cycle, apoptosis and macromolecule catabolic processes (Fig. 5G; Table S1). Furthermore, they were linked to kinase activity (Fig. 5G) as well as to various signaling pathways (Fig. 5H; Table S1). Further analyses revealed that unique genes near to untreated-unique super-enhancers significantly lost their expression upon UV irradiation, indicating that these super-enhancers contribute towards their higher expression in the untreated condition (Fig. 5I, left panel; Fig. S4R, left panel). Genes only occuring near UV-unique super-enhancers showed a slight trend towards gain in expression after UV irradiation (Fig. 5I, middle panel; Fig. S4R, middle panel). Moreover, the genes that retained super-enhancers, on average, did not show any changes in expression (Fig. 5I, right panel; Fig. S4R, right panel). These findings suggest that UV exposure results in a dramatic reorganization of the super-enhancer landscape, which possibly determines gene expression changes likely to be important for the UV-induced stress response.
UV exposure causes a dramatic reorganization of super-enhancers. (A) Bar plot showing total count of super-enhancers in untreated and UV conditions. (B) Venn diagram showing overlap of genes near to super-enhancers in untreated and UV conditions. (C) Box plot showing the absolute distance of super-enhancers to nearest transcription start site (TSS). The box represents the 25–75th percentiles, and the median is indicated. The whiskers show 1.5 times the interquartile range (IQR) added to the 75th percentile (upper whisker) or subtracted from the 25th percentile (lower whisker). The notches represent median±1.57×IQR/(n0.5). P-values are calculated using a Wilcoxon test. (D,E) Pie charts showing genomic distribution of super-enhancers detected in untreated NIH3T3 cells (D) or UV-irradiated cells (E). (F,G) Bar plot showing GO term enrichment for biological processes of genes near to super-enhancers unique to the untreated condition (F) or unique to the UV condition (G). Bar lengths represent the count of genes (primary x-axis), whereas the P-value is shown as a line graph with respect to the upper x-axis. (H). Same as in G, but GO term enrichment is calculated for pathways. (I) Boxplots showing expression of genes near to untreated unique super-enhancers (left), UV unique super-enhancers (middle), and common super-enhancers (right) in untreated (red) and UV conditions (blue). P-values are calculated using Wilcoxon signed-rank test.
UV exposure causes a dramatic reorganization of super-enhancers. (A) Bar plot showing total count of super-enhancers in untreated and UV conditions. (B) Venn diagram showing overlap of genes near to super-enhancers in untreated and UV conditions. (C) Box plot showing the absolute distance of super-enhancers to nearest transcription start site (TSS). The box represents the 25–75th percentiles, and the median is indicated. The whiskers show 1.5 times the interquartile range (IQR) added to the 75th percentile (upper whisker) or subtracted from the 25th percentile (lower whisker). The notches represent median±1.57×IQR/(n0.5). P-values are calculated using a Wilcoxon test. (D,E) Pie charts showing genomic distribution of super-enhancers detected in untreated NIH3T3 cells (D) or UV-irradiated cells (E). (F,G) Bar plot showing GO term enrichment for biological processes of genes near to super-enhancers unique to the untreated condition (F) or unique to the UV condition (G). Bar lengths represent the count of genes (primary x-axis), whereas the P-value is shown as a line graph with respect to the upper x-axis. (H). Same as in G, but GO term enrichment is calculated for pathways. (I) Boxplots showing expression of genes near to untreated unique super-enhancers (left), UV unique super-enhancers (middle), and common super-enhancers (right) in untreated (red) and UV conditions (blue). P-values are calculated using Wilcoxon signed-rank test.
DISCUSSION
It has been known for a long time that there is nucleosome rearrangement in response to UV (Smerdon and Lieberman, 1978) and several studies have followed to investigate how chromatin changes are involved in the UV response, mainly by applying microscopy approaches, biochemical assays or single-loci studies (Kruhlak et al., 2006; Polo, 2015; Yu et al., 2011, 2005). Recent evidence has suggested that chromatin plays a major role in the UV-induced cellular response including regulation of gene expression changes (Andrade-Lima et al., 2015; Dawes et al., 2014; Dinant et al., 2013; Duan and Smerdon, 2014; Hassan et al., 2014; Izhar et al., 2015; Li et al., 2014; Sesto et al., 2002; Zhang et al., 2014). However, a detailed genome-wide analysis to investigate the relationship of chromatin and gene expression changes was still missing. It is well established that distal regulatory elements play a crucial role in regulating gene expression; however, until now, no studies had described the changes in the distal regulatory landscape and its relationship to gene expression changes upon UV irradiation (Creyghton et al., 2010; Rada-Iglesias et al., 2011). Therefore, for the first time, in this study we comprehensively delineate the genome-wide interplay between chromatin accessibility, H3K27 acetylation and gene expression in mouse fibroblasts 6 h after UV irradiation using several high-throughput sequencing methods. In-depth computational analysis revealed a global chromatin compaction in response to UV, which was accompanied by a massive reorganization of the active histone mark H3K27ac at promoters as well as distal sites. Our data therefore highlight a hitherto unknown role of distal regulatory elements in modulating gene expression in response to UV irradiation. Furthermore, our results show that these chromatin changes often reflect expression changes occurring following UV treatment. Overall, our combinatorial analysis of these comprehensive datasets not only provides new insights into UV-induced chromatin remodeling regulating the expression response but also serves as a highly refined resource for the scientific community.
Our results revealed a preferential down-regulation of gene expression accompanied by an up-regulation of translation machinery components. These results might indicate a stress response mechanism in which cells employ alternative mechanisms to keep the relative protein levels constant by avoiding transcription of damaged DNA that might give rise to aberrant transcripts and non-functional proteins. It is also in line with recent findings, that transcription initiation of a majority of transcribed gene promoters is inhibited following UVB exposure (Gyenis et al., 2014). In the recent past, lincRNAs have been shown to play a role in gene regulation and in a plethora of other biological processes (Cheetham et al., 2013; Dinger et al., 2008; Xu et al., 2013). However, until now only a few lincRNAs have been implicated in the DNA damage response (Liu and Lu, 2012; Rashi-Elkeles et al., 2014; Younger et al., 2015). Our transcriptome analysis enabled us to reveal many lincRNAs that are differentially expressed upon UV irradiation, which might be crucial regulators of the stress response and be interesting candidates for further mechanistic analysis. These lincRNAs might act in a similar manner to lncRNAs upstream of the CCND1 promoter, which have been reported to respond to DNA damage signals and to regulate the transcription of CCND1 through recruitment of the histone acetyltransferase (HAT) inhibitor TLS (also known as FUS) to the CCND1 promoter (Wang et al., 2008a).
Previous studies have shown a global chromatin relaxation followed by either restoration or further compaction upon UV irradiation (Burgess et al., 2014; Rubbi and Milner, 2003; Soria et al., 2012). Our findings show a genome-wide loss of chromatin accessibility at 6 h after UV exposure, irrespective of the genomic location, which consolidates our previous observation indicating a potential cellular mechanism to protect genome and transcriptome integrity. The loss of accessibility is accompanied by a gain of H3, indicating that the progressively occurring chromatin compaction is achieved either by active chromatin remodeling or incorporation of new nucleosomes. Although previous studies have suggested that heterochromatic regions are less prone to DNA damage in response to ionizing irradiation, γ-rays or chemical agents (Falk et al., 2008; Seo et al., 2012; Takata et al., 2013; Wei Yu, Genome-wide analysis of DNA damage and repair, Technische Universität Darmstadt, PhD thesis, 2014), our and other published data indicate that UV-induced DNA damage might occur both in eu- and hetero-chromatic regions (Wei Yu, Genome-wide analysis of DNA damage and repair, Technische Universität Darmstadt, PhD thesis, 2014; Zavala et al., 2014).
Interestingly, loss of accessibility was less pronounced at promoters compared to other genomic regions. It is possible that some promoters are actively kept more accessible to maintain transcription or to be able to rapidly reactivate transcription following restoration of a normal cellular state. It is further conceivable that the compaction of promoters requires more energy and time given their occupancy by regulatory proteins as well as transcription machinery and their highly active chromatin status (Lenhard et al., 2012). Despite a genome-wide reduction in chromatin accessibility, a few promoters, including those for genes implicated in the DNA damage response, gain accessibility to allow transcription of the respective genes.
H3K27ac is a histone modification known to mark active promoters and enhancers, and found in our analysis to be drastically reorganized genome-wide following UV exposure (Creyghton et al., 2010; Rada-Iglesias et al., 2011; Shlyueva et al., 2014; Wang et al., 2008b). Furthermore, time-dependent acquisition of active chromatin status along with following expression changes argues that UV irradiation modulates chromatin in order to achieve required transcriptome changes. These observations are in line with previous studies that show changes in active chromatin status along with transcriptional modulation (Suganuma and Workman, 2011; Zhou et al., 2011). Strikingly, our data also hint that the epigenetic status prior to UV exposure influences expression as well as chromatin accessibility changes upon UV irradiation, indicating that a repertoire of epigenetic players cooperate to establish the cellular response to UV exposure.
In line with the global loss of chromatin accessibility, we also observed a dramatic loss of accessible (as determined by FAIRE) and H3K27ac-marked (double-positive) regulatory elements. Motif analysis at untreated-unique double-positive sites indicate that, for example, TFAP2A, KLF4, SP2 and EGR1 might be targeted to these sites to support normal cellular functions. Identification of TP53 at UV-unique double-positive sites is further in line with previous observations showing that TP53 acts as a chromatin accessibility regulating factor mediating UV-induced chromatin relaxation as well as a factor that directly binds to enhancers (Rubbi and Milner, 2003; Younger et al., 2015). Moreover, the presence of CTCF, which has previously been shown to protect cells from UV-induced apoptosis (Li and Lu, 2007), at the non-changing double-positive sites might confer protection against UV-induced alterations of the chromatin status.
Among regulatory regions, we also found that super-enhancers underwent a dramatic reorganization following UV treatment. They further seem to play an important role in transcriptional regulation of nearby genes, following the observation made previously that super-enhancers induce expression of proximal genes (Whyte et al., 2013). Many genes near to UV-induced super-enhancers have previously been shown to be important for the stress response, such as cell cycle, apoptosis and cell survival genes (e.g. Zfp36l1, Tgif, Ctgf and Cyr61) (Thakurela et al., 2015). The super-enhancers might be gained next to them to mainly prevent their down-regulation upon UV exposure or even lead to an up-regulation of their expression. Moreover, UV-induced super-enhancers might be involved in the regulation of many signaling pathways of which some are shown to be activated in response to UV (Engelberg et al., 1994; Lu et al., 2003; Stokes et al., 2007). These findings further suggest that in response to stress, cells rearrange their regulatory landscape to allow desired gene expression programs. This expands our knowledge of the gene regulatory networks induced by UV irradiation.
Taken together, our study serves as a comprehensive resource of how chromatin remodeling of promoters and distal regulatory regions relate to expression changes upon UV exposure and provides new insights into the epigenetic responses to UV damage. Further work should involve investigating the functional role of new target genes, including lincRNAs, in the UV-induced DNA damage response. In addition, given our findings showing that the genomic regions undergoing epigenetic reprogramming might serve as target sites for sequence-specific transcription factors, future work should aim to unravel their functional impact at these sites in response to UV irradiation.
MATERIALS AND METHODS
Cell culture and UV treatment
NIH3T3 cells (NIH/3T3; ATCC® CRL-1658™, ATCC) were cultured at 37°C under 7% CO2 and 88% relative humidity. The culture medium contained Dulbecco's modified Eagle's medium (DMEM) supplemented with 10% fetal bovine serum, 2 mM L-glutamine and 1× non-essential amino acids. At 2 days before the experiment, cells were seeded at a density of 3000 cells/cm2. Then, the culture medium was removed and the cells were washed twice with Dulbecco's phosphate-buffered saline (D-PBS). Any residual solution was carefully removed, and the uncovered cells were irradiated with 80 J/m2 UVC (254 nm) in a CL-1000 Ultraviolet Crosslinker (UVP). Fresh culture medium containing DMSO (1:2000) was added to the cells, which were cultured again for 6 h or as indicated.
ChIP assay and FAIRE assays
ChIP and FAIRE assay was performed as previously described with small adaptations (Sahu et al., 2015).
ATAC
50,000 untreated and UV-treated NIH3T3 cells were applied to the ATAC assay according to Buenrostro et al., 2015. Transposed DNA was amplified for 11 cycles and purified using AMPure XP beads according to the manufacturer's protocol. The DNA was eluted in 15 μl elution buffer (EB, 10 mM Tris-Cl, pH 8.5) and 0.01 μl was used per qPCR measurement.
RNA isolation, cDNA synthesis and quantitative RT-qPCR
Total RNA was prepared using Trizol and reverse transcribed with the First Strand cDNA Synthesis Kit (Thermo Scientific). Transcripts were quantified by PCR using SYBR Green PCR MasterMix on a ViiA7 PCR machine (Life Technologies). The sequences of all primers used in this study are provided in Table S3.
Histone isolation and western blot analysis
Histones were isolated according to Abcam's protocol (Abcam Inc., Cambridge, MA). 10 μg of proteins were run on 15% polyacrylamide gels, transferred onto a PVDF membrane and probed with the respective antibodies.
Dot blot
Human H3.3 peptide containing H3K27ac and the H3S28p modification (JPT, SP-His_0679; aa1.42) were spotted on a 0.45 μM AmershamTM ProtranTM nitrocellulose membrane, blocked in 5% BSA in TBST (0.1% Tween) and then probed with respective antibodies.
Reagents and antibodies
Please refer to Table S4.
RNA-seq
RNA was isolated in biological triplicates with a Purelink RNA Mini Kit and ribosomal RNA was removed using a Ribo-Zero rRNA removal kit. 50-bp reads and single-end sequencing of the RNA-seq libraries, prepared with a TruSeq RNA Library Prep kit, were performed on an Illumina Hi-Seq 2000 platform. Fastq files generated from the sequencing were processed using TopHat (version 2.0.9) with default parameters for alignment to mouse genome mm9 (available from UCSC) resulting in 14×106–22×106 reads per sample (Trapnell et al., 2010). FastQC was used for quality control and HTSeq (version 0.5.4p1) to calculate number of mapped reads on protein-coding genes (n=19,069) or lincRNAs (n=15,694). The normalized expression and differentially expressed genes between the two conditions were identified using DESeq with a FDR of 0.1 (Anders and Huber, 2010). Values showing no tag counts in either condition were excluded from the analysis (resulting in 18,784 protein-coding genes and 8857 lincRNAs).
ChIP-seq and FAIRE-seq
ChIP- and FAIRE-seq experiments were performed in biological duplicates. Fastq files generated by the sequencing platform were processed using Bowtie (version 0.12.7) with default parameters for alignment to mouse genome mm9 (available from UCSC) (Langmead, 2010). A master bam file was generated by merging the bam files from the replicates in both the untreated and UV conditions. MACS2 was used to call the peaks from the merged bam file without input (Zhang et al., 2008). Peak enrichments for any peak given by MACS2 or for specific genomic regions (promoters, exons, introns or intergenic regions) were then calculated using the QuasR package (Gaidatzis et al., 2015). Peaks that were enriched at least 1.5-fold above input were considered for further analysis. Genomic regions (promoters, introns, exons and intergenic regions) were defined using information from UCSC by applying a hierarchical approach. Promoters were defined first at −800 to +200 bp around the transcription start site (TSS), then exon and intron information was extracted from the UCSC gene annotation file and all remaining regions were considered as intergenic. To build the correlation plots, an overlap of 20% was required for two peaks to be at the same site.
Public datasets from the NCBI GEO database used were: SRR1014989, SRR1015026, SRR1015028, SRR1015029, SRR1015032, SRR1187052, SRR1187056, SRR1187061, SRR118706-4-5, SRR1187064, SRR350001.
Gene ontology analysis
DAVID was used to perform functional annotation clustering using biological process and cellular compartments (Huang et al., 2009). From each DAVID cluster, the sub-categories showing the highest number of genes and corresponding P-value were chosen as the representing values. Pathway analysis was performed using ToppGene (Chen et al., 2009). Gene Set Enrichment analysis on differentially expressed genes was performed using GSEA package from GenePattern (Subramanian et al., 2005).
Motif analysis
Motif analysis was performed using the Pscan-ChIP tool applying mouse (mm9) genome assembly while using a mixed background set against JASPAR (Zambelli et al., 2013). To identify the exact genomic locations of the motifs predicted by Pscan, we applied findMotifs.pl of Homer package (Heinz et al., 2010).
Super-enhancer identification
Super-enhancers were identified using the Homer package (Heinz et al., 2010), which applies a similar approach as described in the first paper reporting super-enhancers (Whyte et al., 2013). In brief, peak calling was performed using default options and then peaks within 6 kb were stitched together. A score for super-enhancers was calculated using the total normalized read count in the ChIP sample compared to the input sample. These regions were then sorted and normalized based on the highest super-enhancer score.
Accession numbers
All the next-generation sequencing datasets used in this study have been submitted to GEO and will be publically available under accession number GSE66286.
Acknowledgements
We thank members of the Tiwari laboratory for cooperation and critical feedback during this project. Support by the Core Facilities of the Institute of Molecular Biology (IMB), Mainz, is gratefully acknowledged, especially the genomics and bioinformatics core facilities.
Footnotes
Author contributions
S.S. performed experiments and wrote the manuscript. D.F. analyzed data and wrote the manuscript. S.T. analyzed data. S.K.S. performed experiments. A.G. helped with the progress of the experiments during revision. V.K.T. designed the study, analyzed data and wrote the manuscript.
Funding
Research in the laboratory of V.K.T. is supported by Wilhelm Sander Stiftung [grant number 2012.009.2]; the EpiGeneSys RISE1 program; Marie Curie CIG [grant number 322210]; and Deutsche Forschungsgemeinschaft (DFG) [grant number TI 799/1-1].
References
Competing interests
The authors declare no competing or financial interests.