Chromium disrupts chromatin organization and CTCF access to its cognate sites in promoters of differentially expressed genes

ABSTRACT

Hexavalent chromium compounds are well-established respiratory carcinogens used in industrial processes. While inhalation exposure constitutes an occupational risk affecting mostly chromium workers, environmental exposure from drinking water is a widespread gastrointestinal cancer risk, affecting millions of people throughout the world. Cr(VI) is genotoxic, forming protein-Cr-DNA adducts and silencing tumor suppressor genes, but its mechanism of action at the molecular level is poorly understood. Our prior work using FAIRE showed that Cr(VI) disrupted the binding of transcription factors CTCF and AP-1 to their cognate chromatin sites. Here, we used two complementary approaches to test the hypothesis that chromium perturbs chromatin organization and dynamics. DANPOS2 analyses of MNase-seq data identified several chromatin alterations induced by Cr(VI) affecting nucleosome architecture, including occupancy changes at specific genome locations; position shifts of 10 nucleotides or more; and changes in position amplitude or fuzziness. ATAC-seq analysis revealed that Cr(VI) disrupted the accessibility of chromatin regions enriched for CTCF and AP-1 binding motifs, with a significant co-occurrence of binding sites for both factors in the same region. Cr(VI)-enriched CTCF sites were confirmed by ChIP-seq and found to correlate with evolutionarily conserved sites occupied by CTCF in vivo, as determined by comparison with ENCODE-validated CTCF datasets from mouse liver. In addition, more than 30% of the Cr(VI)-enriched CTCF sites were located in promoters of genes differentially expressed from chromium treatment. Our results support the conclusion that Cr(VI) exposure promotes broad changes in chromatin accessibility and suggest that the subsequent effects on transcription regulation may result from disruption of CTCF binding and nucleosome spacing, implicating transcription regulatory mechanisms as primary Cr(VI) targets.

KEYWORDS:

• chromatin organization
• nucleosome architecture
• ATAC-seq
• CTCF
• hexavalent chromium
• transcription

Abbreviations

5’ UTR	=	5’ Untranslated Region
AP-1	=	Activator Protein-1
ATAC	=	Assay for Transposase-Accessible Chromatin
ATCC	=	American Type Culture Collection
BORIS	=	Brother of the Regulator of Imprinted Sites
ChIP	=	Chromatin Immunoprecipitation
Cr(VI)	=	Hexavalent chromium, represented by K2CrO4
CTCF	=	CCCTC-Binding Factor
DANPOS2	=	Dynamic Analysis of Nucleosome and Protein Occupancy by Sequencing, version 2
DSB	=	Double-Strand Break
ENCODE	=	The Encyclopedia of DNA Elements Consortium
FAIRE	=	Formaldehyde-Assisted Isolation of Regulatory Elements
FDR	=	False Discovery Rate
MNase	=	Micrococcal Nuclease
NFR	=	Nucleosome-Free Regions
NUC	=	Nucleosomal Insertions
TAD	=	Topologically Associated Domain
TSS	=	Transcription Start Site

Introduction

Hexavalent chromium compounds are well-known respiratory carcinogens, widely used in industrial processes, including among others, electroplating, stainless steel production, and welding. Numerous epidemiological and targeted studies have shown a high incidence of lung cancer among workers chronically exposed to chromium by inhalation or dermal contact [1,2]. Inhalation and dermal contact exposure, however, affect only a small fraction of the population, mainly some half-a-million occupationally exposed individuals. In contrast, anthropogenic and natural means spread chromium contamination beyond the industrial sites of intended use, causing the exposure of the general population to low levels of this metal, primarily through food and water consumption [3]. As a consequence, oral exposure by ingestion is widespread, affecting millions of people throughout the globe, causing gastrointestinal effects that include indigestion, stomach and abdominal pain, ulcer formation, gastritis, gastrointestinal cancer, and increasing the risk of many other types of human cancers [3–6]. Mouse studies in our laboratory, among others, have shown that chronic ingestion of low Cr(VI) levels in drinking water lead to enterocyte hyperplasia, a possible compensatory phenotype resulting from chronic villous wounding [7–9]. A large number of animal studies support the conclusion that the molecular properties of hexavalent chromium result in carcinogenic outcomes regardless of the exposure route.

Despite this well-characterized phenotype, the molecular mechanisms of action that potentiate tumorigenesis from chromium exposure remain unclear. Chromium is genotoxic but appears to be non-mutagenic [10]. Toxicity from chromium exposure is associated with chromium (VI) compounds, even though chromium (III) is the most prevalent form in the environment and in biological tissues [11]. Cr(VI) enters cells via the sulfate anion transporter system [12] and is reduced through 1- or 2-electron reduction steps to its intermediate oxidation states, Cr(V), Cr(IV), and ultimately Cr(III), by ascorbic acid, glutathione, and cysteine [13–15]. As the product of Cr(VI) reduction, Cr(III) causes radical-mediated DNA damage [16], double strand DNA breaks [17,18], increasing γH2A.X foci in euchromatic regions, loss of tumor suppressor gene expression, and altered transcriptional responses [19–23], disrupting signaling pathways that might be involved in transcriptional inhibition [24]. Cr(VI) reduction allows it to bind to several reactive amino acids, including cysteine, tyrosine and histidine, causing the formation of Cr-DNA adducts and DNA-Cr-protein crosslinks [25–32]. Chromium (VI) has also been shown to persistently activate the MAP kinases ERK1 and ERK2 by a redox-sensitive mechanism [26,27] and to induce the activation of NF-κB in cultured Jurkat cells by generation of free radicals from the reduction of Cr(VI) [28,29]. These observations suggest that the primary Cr(VI) targets involved in the progression of tumorigenesis might be the regulatory mechanisms responsible for guiding proper transcriptional programming.

To develop an understanding of the molecular mechanisms of action that potentiate the health effects of chromium exposure, we previously carried out targeted and genome-wide studies of the transcriptional regulatory changes mediated by Cr(VI) in the Hepa-1c1c7 murine hepatoma cell line. We found that Cr(VI) crosslinked the chromatin remodeling complexes of HDAC1-DNMT1 to the proximal promoter of several inducible genes, inhibiting their expression by maintaining a closed chromatin state and preventing RNA polymerase II recruitment [30–32], suggesting that Cr(VI) disrupts chromatin architecture, a key component in the regulation of transcription. To examine the impact of Cr(VI) on global chromatin accessibility, we used FAIRE-seq after treatment with long-chronic or short-acute chromium concentrations. At the low concentration, accessibility correlated with increased gene expression mainly from AP-1 promoters, while at the high concentration, accessibility correlated with repression, with approximately 80% of the identified peaks representing downregulated genes enriched for domains containing predominantly CTCF motifs [33]. The differential accessibility of AP-1 and CTCF following Cr(VI) treatment suggested that transcriptional perturbations might be potentiated by chromium disrupting key factors associated with genomic integrity and hence affecting epigenomic regulation through chromatin re-organization and nucleosome dynamics.

To test this hypothesis, we utilized two complementary approaches, MNase-seq and ATAC-seq, to measure changes in dynamic nucleosome positioning and native chromatin accessibility, respectively, in response to Cr(VI) treatment. We found that chromium induced changes in nucleosome occupancy and fuzziness (amplitude) and that differentially accessible domains exhibited a gradual increase dependent on chromium concentration. In agreement with our previous findings, AP-1 and CTCF motifs were highly enriched in every condition tested and frequently found within promoters of differentially expressed genes. Our findings suggest that Cr(VI) exposure disrupts the normal spacing of accessible nucleosomes and inter-nucleosomal chromatin links surrounding CTCF motifs and highlight the need to further investigate direct and indirect interactions that may exist between chromium ions and CTCF.

Results

Cr(VI) treatment alters the occupancy of dynamic nucleosomes in promoter regions

MNase-seq analysis, using the Dynamic Analysis of Nucleosome and Protein Occupancy by Sequencing, version 2 (DANPOS2) Toolkit [34], sorts dynamic nucleosomes into three categories based on the strength of their occupancy at a given site, positional shifts of 10 or more nucleotides, and fuzziness, a measure of signal amplitude. Using a false discovery rate (FDR) threshold of 0.01, dynamic nucleosome profiles in Hepa-1c1c7 cells treated with 2 µM potassium chromate for 72 hours were called representing changes in occupancy, position, and fuzziness in Cr-treated cells compared to control. Differential occupancy of nucleosome position was the most frequent call, comprising 91,192 sites, followed by fuzziness, comprising 33,088 sites, and 2,438 position shifts, with relatively few nucleosome positions fitting two criteria and none all three (Figure 1A). We used HOMER [35] to identify the genomic domains corresponding to dynamic nucleosome changes sorted by three criteria: total pool (whole genome), intergenic regions [>1.5 Kb from transcription start site (TSS)], and promoter regions (≤1.5 Kb from TSS). Annotation of these regions identified CTCF as an enriched target in promoter regions that exhibited differential occupancy following Cr(VI) treatment (Figure 1B). The impact of exposure on dynamic nucleosome positions surrounding CTCF motifs in promoter regions was visualized by using mean signal plots for nucleosome positions extended ±1.5 Kb in either side of the labeled motifs after accounting for directionality in CTCF motif orientations, a feature tied to functionality and chromatin looping [36] (Figure 1C). CTCF sites in promoter regions that exhibited changes in occupancy were typically associated with nucleosomes immediately adjacent to the motif, a feature that was not seen at the level of the whole genome (Figure 1B). Given that most CTCF sites occur outside of promoter regions and not all are likely to be affected by chromium, our findings are not unexpected, since the signal for whole genomic sites affected by Cr(VI) may be dampened or lost when averaged with sites that are otherwise unaltered.

Figure 1. Occupancy of dynamic nucleosomes in promoter regions affects CTCF motif accessibility. Hepa-1c1c7 cells were treated with either H₂O or 2 µM K₂CrO₄ for 72 hours before collection and monosome preparation by MNase digestion. (A) Represents dynamic nucleosomes that were called with the DaNPOS2 toolkit [34] with a threshold of q = 0.01. (B) Dynamic nucleosomes were analyzed for motif enrichment using HOMER [35] and separated based on genomic annotation at two different thresholds, q = 0.1 and q = 0.01. Enriched motifs are provided as the -log₁₀ P value of the known motif, with its rank compared to other motifs listed in brackets. (C) CTCF motifs identified in promoters based on changes in nucleosome occupancy were extracted and motif orientation was corrected using previous annotations prior to generating average signal plots of nucleosome positions with deepTools2 Suite [57]. (D) The bidirectional promoter region between Trove2 and Uchl5 is provided to illustrate the complementarity between different sequencing techniques. Signal tracks in blue (control) and red (Cr-treated) represent either the nucleosome position generated using MNase-seq (upper panel) or ATAC-seq insert signal (lower panels). ATAC-seq tracks are represented using three separate criteria; “composite” is the average signal recorded for each bp with no insert-size filtering, while the middle track represents short inserts (less than 100 bp in length), and the lowest nucleosome-spanning inserts (greater than 180 bp in length).

To confirm the findings of nucleosome positioning analysis, we used ATAC-seq, a sensitive measurement of chromatin accessibility that uses transposase insertion to identify changes in native chromatin accessibility [37]. Treatment with low chromium concentrations, ranging from 0.1 to 1.0 µM showed a good complementarity between nucleosome positioning and ATAC-seq-defined chromatin accessibility. An example of this correspondence is shown in Figure 1D for the bidirectional promoter domain of the Trove2 and Uchl5 genes in chromosome 1. Within a narrow range of 227 bp, a CTCF-site located between the two TSS, was depleted of nucleosomes and differentially accessible following exposure to Cr(VI), strongly indicative of altered CTCF binding at this location.

ATAC-seq identifies broad, chromium-dependent changes in chromatin accessibility

Chromium-dependent differentially accessible regions were called using an FDR of 0.0001, which identified 14,401, 11,130, and 15,894 peaks with an average length of 417 ± 249 bp at 0.1, 0.5, and 1.0 µM chromium, respectively. Our previous work using FAIRE-seq had shown distinct accessibility profiles for regulatory elements that were dose-dependent [33]; therefore, we compared peaks for intersecting segments to determine whether Cr(VI) affected chromatin accessibility consistently at multiple low doses. After intersection analysis, the total number of segments for the peaks at 0.1 µM was 34,602, while the peaks identified at 0.5 µM and 1.0 µM Cr consisted of 22,186 and 39,868 segments, respectively, with an average of 2.3 ± 0.3 segments per peak. The calculations for multiple intersections result in the segmentation of close, unique peaks, with most peaks consisting of several partially overlapping segments between treatments, while retaining a substantial number of unique, differentially accessible regions, as shown in the Venn diagram in Figure 2A. Of note, a core set of 7,720 segments common to all three Cr treatments were consistently identified when compared to control. This trend indicates that specific regions may be susceptible to changes in accessibility following low doses of chromium with a certain degree of stochasticity, as intersecting peaks frequently resulted in overhanging segments, a potential indication of slight variations in the localization of the structure.

Figure 2. ATAC-seq identifies changes in chromatin accessibility associated with Cr(VI) exposure. Hepa-1c1c7 cells were treated with 0.1, 0.5, and 1.0 µM of K₂CrO₄ for 72 hours followed by nuclei isolation and addition of the Tn5-transposition reaction. (A) Differentially accessible peaks were called against an untreated control using MACS2 [52] with a threshold of q = 0.0001 and overlaps were calculated using the multiIntersectBed command in the Bedtools Suite [53]. The total number of input peaks for each treatment were 14,401, 11,130, and 15,894 for Cr 0.1, 0.5, and 1.0 µM prior to measuring the intersections. (B) Normalized log₂ signal tracks comparing treated and control were generated using deepTools2 [57] with a bin size of 50 bp. Mean signal values were measured across normalized peak distances set to 1,000 bp for their respective treatments. Regions were considered to be significantly opened or closed following treatment if the mean was R ≥ 0.5 or R ≤ −0.5, respectively. (C-D) Each set of differentially accessible peaks was analyzed using HOMER [35] for annotation statistics and identification of significantly enriched motifs, listed as the –log₁₀ P value of the known motif for ease of description [35]. The rank of enrichment for each motif logo is provided in brackets.

To estimate the degree of change following treatment, signal tracks representing the log₂ ratio between treatment and control were used to measure the degree of accessibility for each set of peaks relevant to their respective treatment. Regions were considered to be open if the mean of the log₂ signal was greater than 0.5 and closed if smaller than -0.5. Interestingly, the lowest Cr(VI) concentration (0.1 µM) identified the largest number of opened regions, 10,369 out of 14,401 total regions, followed by 2,733 in the 0.5 µM concentration, and 6,073 in the 1 µM samples (Figure 2B). Closed regions appeared in a minimal number of cases, with a trend towards moderate changes as the concentration increased. This pattern suggests that increasing the Cr(VI) treatment affects the degree of chromatin permissibility, either through direct action of Cr-mediated interference or through a cellular response to exposure.

To determine whether differential accessibility was associated with specific domains, we measured annotation statistics and motif enrichment for the global peak sets using HOMER. Regions associated with transcription, specifically those located near the promoter and 5UTR, were highly enriched compared to intergenic noncoding regions (Figure 2C). Motif enrichment analyses identified general similarities comparable to the profiles found in our previous FAIRE-seq results. At an FDR of 0.01, the subset of highly enriched motifs for each treatment consisted of AP-1 at ranks 1 and 2, followed closely by CTCF at ranks 4, 5, and 4 for the increasing concentrations; at the more stringent FDR of 0.0001, CTCF ranked at 6, 19, and 5 (Figure 2D).

Proximal genomic regions of differentially accessible CTCF sites exhibit increased AP-1 motif density

The AP-1 transcription factor is involved in the activation and regulation of multiple genes in signaling pathways related to stress responses, cell proliferation, and apoptosis. Cells treated with Cr(VI) have been shown to activate AP-1 and contribute to cytotoxicity through specific pathways of MAP kinase signaling [26,38]. Consistent with our previously established findings, the AP-1 motif and its related motif, Jun-AP1, were the most enriched sites observed in ATAC-seq identified peaks (Figure 2D). To determine whether AP-1 and CTCF sites consistently co-occurred in the enriched peaks, differentially accessible CTCF sites were extracted from each peak set and used to assess motif densities of Jun-AP-1/AP-1 within ±500 bp of an identified CTCF site. As an internal control we tested the co-occurrence with BACH2, the next ranked motif, and the association with an independent determination of potential co-occurrence of the same sites in the CTCF ChIP-seq ENCODE dataset. Interestingly, domains proximal to differentially accessible CTCF sites were frequent contributors to Jun-AP-1/AP-1 motif enrichment compared to ENCODE-identified sites, while BACH2 displayed minimal changes in co-occurrence (Figure 3). These data are consistent with possible functional associations between CTCF and AP-1.

Figure 3. Proximal regions surrounding differentially accessible CTCF sites exhibit increased AP-1 motif density. CTCF sites were extracted from each set of differentially accessible peaks and subsequently used to measure the density of AP-1, Jun-AP1 and BACH2 motifs relative to each site ± 500 bp using HOMER [35]. The top panel represents the average frequency score for each motif in 20-bp bins. The panel below shows a smoothed model for the same data, generated using the smooth feature and generalized additive model method in ggplot2 [65] with no specified bin size. Jun-AP1 and AP-1 exhibit very similar trends among all conditions, thus only the Jun-AP1 results are shown. The number of CTCF sites represented in each graph from left to right is 1,354; 862; 1,525; and 21,901 respectively. ENCODE motif locations were derived from the adult CTCF ChIP-seq data, ENCSR000CBU.

Chromium causes differential accessibility to CTCF sites in gene promoters associated with changes in gene expression

CTCF has been shown to play a pivotal role in several critical cellular processes including transcription regulation, epigenetic signal insulation, and nuclear organization of topologically associated domains (TADs), often in conjunction with the cohesin complex [39,40]. Following this lead, we chose to examine whether sites with differentially accessible CTCF motifs in promoter regions displayed changes in expression. Chromium treatment elicits broad changes in chromatin accessibility associated with concentration-dependent changes in gene expression [23]. Annotation of CTCF sites identified by ATAC-seq indicated that between 809 and 1,428 of these sites were mapped as being in a promoter domain (±1.5 Kb from TSS), an intergenic domain (>1.5 Kb from TSS in a noncoding region), or an intragenic domain (>1.5 Kb from TSS in a coding region) (Figure 4A). After Cr(VI) treatment, 25–30% of the sites mapped to promoter domains, suggesting a possible flux in the transcriptional status of the associated genes. To investigate these potential changes in gene expression, we compiled a set of 445 unique genes (Supplemental Table 3) with differentially accessible CTCF sites within their promoter regions and used it to query our previously published RNA-seq data of Hepa-1c1c7 cells exposed to either a low concentration, chronic Cr(VI) treatment (0.5 µM) [23] or an acute, high concentration (25 µM) [33]. We did not find a significant correlation between CTCF accessibility to gene promoters and preference for gene repression or induction at either chromium concentration; rather, approximately 50% of the genes affected at either concentration showed increased levels of expression and the remaining 50% decreased. Several genes involved in tumor suppression or DNA repair processes, including Rev1, Rad1, and Cenpc, exhibited notable concentration-dependent changes in expression, as well as genes in the cohesin complex, such as Smc3 and Stag1, which exhibited strong downregulation at both concentrations compared to non-treated controls (Figure 4B).

Figure 4. Differentially accessible CTCF motifs in promoters are associated with changes in expression. (A) The table describes the annotation statistics and CTCF enrichment analysis from HOMER [35]. Peaks were separated into “Promoter”, “Intergenic”, or “Intragenic” and compared to the total number of differentially accessible regions per treatment to estimate the percentage of peaks in each category. Additionally, the sum of peaks per category that contained one or more CTCF motifs are provided. (B) Peaks positive for the CTCF motif within ±1.5 Kb of a transcriptional start site (TSS) were selected and the closest gene ID was used to generate a list of 923 potential genes of interest. These were filtered to remove duplicate genes called by each treatment resulting in a total of 445 unique genes (see Supplemental Table 3). The log₂ fold change values from previously published RNA-Seq datasets using doses of 0.5 µM or 25 µM K₂CrO₄ were used to assess changes in transcription for the filtered set of genes [23,33]; 400 of the 445 genes were annotated with expression values and each condition was independently analyzed back to controls in its respective experiment.

CTCF binding motifs with disrupted accessibility due to Cr(VI) treatment are evolutionarily conserved sites

Numerous studies have shown that changes in CTCF binding patterns may result in severe adverse consequences to chromatin structure. Comparisons between tumor-derived cell lines and their normal counterparts have shown that tumorigenesis can reshape local chromatin organization and alter gene expression through changes in CTCF binding patterns [41]. To assess whether the differentially accessible CTCF sites induced by chromium had a physiological relevance we used two separate approaches. First, we measured their degree of conservation by calculating the per-base phastCons score for all unique CTCF motifs ± 40 bp (an 80 bp window) identified by ATAC-seq (Figure 5A). PhastCons provides a relative measure of DNA-element conservation derived from comparisons between 60 vertebrate species (follow link to UCSC Browser in Supplemental Table 2). The average score for our 2,044 mappable elements was approximately 0.6; using the same method, the 27,168 unique CTCF motifs derived from three published ENCODE mouse liver datasets returned an average conservation score of roughly 0.5, indicating that the majority of the affected CTCF-elements identified by ATAC-seq were moderately to highly conserved (Figure 5A). Second, we used the ENCODE ChIP-seq datasets from E14.5, PND0 and 8-week adult livers to compare the specific locations of CTCF motifs within peaks from our in vitro sites to those found in vivo. A total of 19,581 commonly occupied CTCF motifs were identified, suggesting the existence of a core set of binding sites important in the development and maintenance of the liver (Figure 5B). To examine whether Cr treatment affects these functional sites, we used CTCF ChIP-seq analysis after a short, acute treatment with 25 µM potassium chromate for 90 minutes and compared the peaks called between control and treatment to the set of 19,581 common motifs in ENCODE. We found that 734 sites were common between each condition with an additional 2,369 motifs being called in either treatment. Of note, a much smaller set of motifs, between 520 and 424, were unique to our ChIP-seq samples (Figure 5C).

Figure 5. CTCF sites with altered accessibility are conserved, functional sites in the mouse liver. (A) The locations of CTCF motifs (20 bp sites) within each set of peaks was determined using HOMER [35], then pooled and subsequently filtered for unique values to generate a comprehensive list of CTCF sites affected by Cr(VI) treatment. DeepTools2⁵⁷ was used to compute the conservation scores in the mm10 phastCons bigwig file for each motif and its surrounding region ±40 bp. Regions with no value were not included. The summary plot represents the mean phastCons score per bp. (B) Conservative, IDR-thresholded peak files from three mouse liver CTCF ChIP studies in ENCODE were obtained and scanned for the same 20 bp CTCF motif (Find accession numbers in Supplemental Table 1) [61,62]. Locations were extracted and measured for overlap as previously described. (C) de novo motif locations called in a separate ChIP-Seq experiment using 25 µM Cr(VI) were compared for overlapping sites in the core, common set of ENCODE-derived motifs (19,581).

Chromatin surrounding CTCF motifs is increasingly accessible after Cr(VI) treatment

The use of Tn5 Transposase in ATAC-seq generates inserts of varied lengths capable of informing about the nucleosome-level characteristics of a targeted region. Depending on the nucleosome positioning and degree of chromatin accessibility, paired transposition events may occur within the length of DNA between nucleosomes as well as traverse nucleosomes of open chromatin, providing insight into their positions and occupancy.

To understand how increasing concentrations of Cr(VI) affected the extent of chromatin accessibility surrounding CTCF sites, reads from ATAC-seq were divided into “nucleosome-free” (NFR) and “nucleosome-spanning” insertions (NUC). NFRs included 7-12 million insertions of less than 100 bp in length, and NUCs included 12-16 million insertions greater than 180 bp, accounting for 146 bp of nucleosomal DNA plus 20 bp at either side of the nucleosome. To account for differing numbers, each set of reads was normalized to 1X sequencing depth prior to measuring the signal across differentially accessible CTCF sites in the 1 µM chromium concentration for control and treated conditions (Figure 6A). The NFR signal surrounding each set of sites was relatively unchanged and was consistent with the enrichment of shortened fragments in CTCF sites as previously observed by others [37]. In contrast, the signal of NUC fragments was enriched in the 1 µM Cr(VI) concentration and exhibited stronger average values across its differential accessible CTCF motifs when compared to the control. To ensure accuracy of our interpretation, the mean number of transposition sites per bp surrounding differentially accessible CTCF motifs was calculated using two ranges, ±100 and 500 bp, respectively (Figure 6B). In each set of criteria, the base pairs directly involved in the motif exhibited a strong reduction in the number of transposition events, suggesting that CTCF was bound. Following Cr(VI) treatment, the proximal regions immediately surrounding CTCF exhibited an increased number of mean transpositions in the NFR insertions despite minimal changes occurring at the motif itself. Analysis of the distal regions (±500 bp) revealed no discernible differences, suggesting that insertion events are localized to the region immediately surrounding CTCF. These findings strongly support the conclusion that low dose (1 µM) of Cr-exposure alters chromatin structure through disruption of the correct spacing between nucleosomes.

Figure 6. Cr(VI) disrupts normal chromatin accessibility profiles. Bigwig files representing the enrichment of nucleosome free regions (NFR) and nucleosome-traversing (NUC) insertions were generated from bam files using deepTools2⁵⁷. Fragments less than 100 bp were considered as nucleosome free regions (NFR, ∼7-12 M reads) and those greater than 180 bp were considered as nucleosome-traversing (NUC, ∼12-16 M reads), based on previously established cutoffs [37]. Following allocation, signal tracks were normalized to 1X sequencing depth. (A) The signal surrounding 1,524 differentially accessible CTCF motifs in the 1 µM Cr(VI) treatment compared to control was calculated using a bin size of 10 bp for the motif ± 500 bp, then plotted as a heatmap to examine trends. (B) Atactk [60] was used to calculate the mean transpositions per bp surrounding CTCF motifs identified in the 1 µM Cr(VI) treatment. Control and treated bam files were queried using CTCF motifs allocated based on their orientation and subsequently corrected to provide a single, forward-facing orientation. Each bam file was tested using two separate bin criteria at single base-pair resolution (run simultaneously), noted as NFR (1-100 bp) and NUC (180-1,000 bp). Positive and negative strand cut-point values were summed and subsequently used to calculate a mean cut point value based on their initial means.

Discussion

In this study, we show that Cr(VI) disrupts chromatin organization and perturbs the accessibility of CTCF to its cognate sites. Disruption is not unidirectional: in some chromatin regions sites that were opened in control cells become inaccessible after chromium exposure and in other regions the reverse is true. Results from our earlier work using FAIRE suggested that the effect of Cr(VI) on chromatin was closely associated with changes in the transcriptome [33] and identified broad changes in chromatin accessibility that were dependent on dose and length of exposure. The work reported here confirms and extends those earlier conclusions and shows that changes in chromatin structure resulting from Cr(VI) treatment disrupt nucleosome occupancy and frequently alter CTCF's accessibility to its cognate binding sites in regulatory elements.

Many of the changes in chromatin accessibility induced by Cr(VI) were common to all three metal concentrations tested, suggesting that these were likely to be specific, deterministic effects of chromium. However, the majority of the changes observed were uniquely found in cells exposed to one or another chromium concentration in large part due to small degrees of shifts between overlapping peaks that result in overhanging segments, supporting the contention that these structural changes induced by Cr(VI) are determined with a degree of stochasticity. Notwithstanding, our results indicate that chromatin changes induced by Cr(VI) largely occur in functional regions associated with transcription, and show strong enrichment for promoters, suggesting that Cr-induced differential chromatin accessibility impacts the transcriptome. For the purpose of this article, we have followed our previously published work using FAIRE [33] and focused on Cr(VI)-mediated effects on CTCF-site accessibility. Interestingly, ATAC-seq analysis of differentially accessible regions clearly detected the enrichment of AP-1 and CTCF motifs, as previously found using FAIRE-seq [33]. Given recent findings that characterize a multivalent set of functions for CTCF as a key regulator of transcription and chromatin organization [42], it is likely that CTCF sites are Cr(VI) targets that disrupt important epigenetic regulatory mechanisms required for maintaining proper gene expression profiles. Additionally, recent studies have proposed a role for AP-1 in the formation of activation hubs and dynamic chromatin loop formation associated with transcriptional regulation [43]. Our findings from ATAC-seq analysis show that Cr(VI) treatment impairs accessibility to CTCF sites with a significantly increased density of co-occurring AP-1 sites, but not BACH2, another highly enriched motif, or the global control with ENCODE data. Both AP-1 and CTCF have been associated with multiple chromatin remodeling systems [44], consistent with our hypothesis that Cr(VI) treatment elicits broad structural changes in chromatin. In our study, the binding status of these CTCF-proximal AP-1 sites is unknown and it remains unclear whether or not a functional association between AP-1 and CTCF exists in chromium-mediated differentially accessible regions. Interestingly, it has been suggested that AP-1 recruitment can reorganize local chromatin resulting in the eviction of CTCF [45,46]; it is therefore plausible that chromatin remodeling events driven by one may impact binding of the other, and warrants further investigation to examine whether Cr(VI) targets this process in rerouting the transcriptional status of regulated genes.

To test the potential connections between CTCF, Cr(VI), and transcription we searched for genes with differentially accessible CTCF sites in their promoter among the genes in our previously published RNA-seq data [23,33]. We found that these genes frequently exhibited changes in expression. Several of these genes are involved in processes associated with chromatin maintenance and DNA repair and genes required for efficient loading of the cohesin complex onto chromatin, sister chromatid cohesion, DNA replication and repair, and genomic compartmentalization [47]. Given the strong association between CTCF and cohesin, it is plausible that broad changes in transcription associated with Cr(VI) exposure are, in part, due to disruption of the mechanisms responsible for chromatin structure, organization, and regulation. Indeed, depletion of the STAG1 subunit of cohesin in mice leads to the development of tumors and aneuploidy [48], also a hallmark of prolonged Cr(VI) exposure [19].

Changes in CTCF binding profiles between normal and tumorigenic hepatocyte cell lines are well known [41] and functional CTCF/cohesin sites may be targets for mutation [49]. This consideration raises the concern that the Cr-dependent sites identified in tumorigenic Hepa-1c1c7 cells might not be relevant to physiological functions in normal cells. To address this concern, we measured the degree of overlap between Hepa-1c1c7 control and Cr-treated CTCF ChIP-seq samples, following two separate approaches. First, we calculated the average conservation score of differentially accessible motifs, using the assumption that highly conserved sequences are likely to play crucial roles in function. Second, we used ENCODE-validated mouse liver CTCF ChIP-seq samples to derive a set of core constitutively bound motifs present for different developmental time points. In each case, we found that the CTCF motifs affected by Cr(VI) were moderately to highly conserved and that the majority of these sites were also present in the ENCODE datasets. These findings point to the possibility that conserved, physiologically relevant CTCF sites may be deterministic Cr(VI) targets.

Currently, there is little insight into the mechanisms of action of environmental toxicants at the level of three-dimensional chromatin organization. The work described here highlights the need of further analysis using targeted approaches to improve our understanding of how changes in genomic accessibility surrounding CTCF binding due to Cr(VI) exposure affect the regulatory functions of the transcriptome, and whether these actions contribute to toxicity. Our findings provide further evidence that environmental exposures in general, and to Cr(VI) in particular, impact the regulation of the transcriptome and suggest that disruption of chromatin organization may be an important factor in understanding environmental effects on health and disease.

Materials and methods

Cell culture

Hepa-1c1c7 (mouse hepatoma, ATCC) were cultured in alpha-minimum essential medium (Gibco, 12000-022) supplemented with 5% fetal bovine serum (Sigma, F2442) and 1% penicillin-streptomycin (Gibco, 15140-122) at 37°C under 5% CO₂ conditions. Cells were passaged, treated, and collected between 80–90% confluency.

Reagents

Potassium chromate (Fisher Chemical, P220-500), hereinafter simply referred to as chromium, was freshly prepared in a volume of double-distilled H₂O followed by filtration through a 0.22 µm filter prior to use in cell culture. Anti-CTCF antibodies were used for ChIP- analyses (Cell Signaling, 2899S). A list of buffers used for each experiment can be found in Supplemental Table 1.

ATAC-seq

Hepa-1c1c7 cells were treated with 0.1, 0.5, or 1.0 µM potassium chromate for 48 hours before preparation for sequencing using a protocol adapted from the original by the Greenleaf lab [37]. In brief, cells were harvested immediately after chromium treatment, pelleted, washed with 1X PBS, and lysed in cell lysis buffer. Nuclei were re-suspended in transposition buffer containing Tn5 Transposase (Illumina, FC-121-1030) and incubated at 37°C for 30 minutes before eluting transposed DNA using Qiagen™ MinElute Kit (Qiagen, 28004). Following elution of transposed DNA fragments, the Nextera library was generated using an Illumina Nextera DNA Library Prep Kit (San Diego) after 9 PCR cycles for library amplification and indexing. Following amplification, samples were sequenced using an Illumina HiSeq platform to generate ∼75 M pair-end 2 × 100 bp reads. Adapter sequences were trimmed from FASTQs and the paired-end reads were subsequently aligned to the mm10 genome using BOWTIE2 [50]. PICARD Tools [51] was used to remove duplicates and estimate the library size. Reads were filtered for a minimum alignment quality score greater than 30, paired-end matching, and those mapping to the mitochondria, unmapped contigs, and the Y chromosome were removed.

Nucleosome Positioning / MNase-Seq

Hepa-1c1c7 cells were treated with 2 µM potassium chromate for 72 hours, replacing the medium daily. When 80% confluent, cells were washed twice with ice-cold 1X PBS and lysed in homogenization buffer to release nuclei, which were re-suspended in washing buffer and treated with 60 units of MNase for 6 minutes. Nucleosomes were separated by electrophoresis in 1% agarose gels and DNA was purified from mononucleosomes, deproteinized with Proteinase K (Roche, 03-115-879-001), and prepared for sequencing using the same PrepX DNA Library kit and Apollo 324 NGS automatic library prep system (Takara Bio USA, 400075): see detailed protocol for ChIP-seq below. The ligated library was enriched by 6 cycles of PCR. Samples were sequenced using an Illumina HiSeq platform to generate ∼360 M pair-end 2 × 100 bp reads.

ChIP-seq

ChIP-seq was performed by methods previously described with slight modifications [31]. Cells were fixed in 1% formaldehyde for 10 min at 37°C, quenched with 0.125 M glycine for 10 min at room temperature and subsequently collected in PBS containing protease inhibitors (Roche, 05-056-489-001) and allowed to lyse on ice for 10 min. The nuclei were pelleted, resuspended in nuclear lysis buffer and sonicated using a Diagenode Bioruptor (Diagenode, UCD-200) for 6 sets of 6 minutes, 30 seconds on/30 second off cycles on the High setting. An aliquot was taken from each sample and pooled into two replicates to represent total input fractions. Sonicated chromatin was diluted with IP dilution buffer, divided into equal aliquots and incubated overnight at 4°C with the appropriate antibody. Antigen-antibody complexes were collected with MagnaChIP A/G beads (EMD Millipore, 16-663) and washed using 1X dialysis buffer and IP wash buffer before elution. The eluents were treated for 2 hours with 10 µg/ml RNase A (Sigma, R5503) and NaCl at a final concentration of 0.3 M to reverse crosslinks, followed by 30 µg/ml Proteinase K at 56°C before DNA enrichment using the ChIP DNA Clean and Concentrator Kit (Zymo, D5205). Purified DNA samples were quality checked and subsequently sequenced. The quality of the DNA fragments was analyzed by Bioanalyzer DNA High Sensitivity chip (Agilent, 5067–4626), which showed the expected size distribution (∼200 bp peak size). DNA concentrations were measured by Qubit Fluorometric Quantitation using dsDNA HS Assay Kit (Thermofisher, Q32851). A ChIP-seq script was selected for sequencing library preparation using PrepX DNA Library kit and Apollo 324 NGS automatic library prep system (Takara Bio USA, 400075). Overhangs in samples of double-strand DNA fragments (∼5 ng) were first blunt-ended, and adenylated at the 3' ends for TA ligation to the sequencing adaptors. The ligated library was enriched by 8 cycles of PCR using index-specific primers, followed by automated AMPure XP beads (Beckman Coulter, A63882) purification using the Apollo 324 system. To check the quality and yield of the library, one microliter of purified PCR library was analyzed by Bioanalyzer with DNA High Sensitivity chip. To accurately quantify the library concentration for cluster generation, the library was 1:104 diluted with dilution buffer (10 mM Tris-HCl, pH 8.0 with 0.05% Tween 20), and qPCR analyzed with NEBNext Library Quant Kit (New England Biolabs, E7630L) using ABI 9700HT real-time PCR system (Lifetech, 4329001). Individually indexed libraries were proportionally pooled for clustering in the cBot system (Illumina, SY-301-2002). Libraries at the final concentration of 16 pM were clustered onto a flow cell using Illumina TruSeq SR Cluster kit v3 (Illumina, GD-401-3001), and sequenced for 50 cycles using TruSeq SBS Kit (Illumina, FC-401-3002) on Illumina HiSeq system (SY-405-1001). As recommended by Illumina for studies targeting transcription factors, we generated ∼10 M reads per sample.

Peak calling, differential accessibility, and conservation analysis

Differentially accessible peaks were called on processed bam files compared to non-treated controls using MACS2 [52] and the outputs were filtered using ENCODE's consensus exclusion blacklist (see Supplemental Table 2). Nucleosome Positioning results were prepared and processed using the DANPOS2 Suite [34]. Overlap analysis and general functions related to bed files were performed using the BEDTools Suite [53]. Annotation of peak files and extraction of differentially accessible motif locations was performed using HOMER [35]. Additionally, CTCF motifs used to generate signal plots of surrounding nucleosomes and transpositions were re-oriented to face the same direction based on annotations provided by HOMER [35]. The mm10 phastCons [54,55] bigwig file was downloaded from the UCSC Genome Browser [56] and used to calculate the conservation scores of differentially accessible CTCF motifs identified in ATAC-seq. Data preparation and normalization for signal tracks relating to ATAC-seq was performed and visualized using the deepTools2 Suite [57] and the Integrative Genomics Viewer [58,59]. Analysis of transposition events was performed using Atactk [60].

Data accessibility

A list of the datasets used can be found with their ENCODE [61,62] and/or GEO accession ID [23,33,63,64] in Supplemental Table 2 along with current links to programs and additional files. Genome-wide ATAC-seq, MNase-seq, and ChIP-seq data have been submitted to the GEO database with accession number GSE104566 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104566).

Disclosure of potential conflicts of interest

The authors declare no conflict of interest.

Acknowledgements

We thank Dr. Ying Xia, Dr. Chia-I Ko, Dr. Hongxia Zhang, Yunxia Fan, and Matthew de Gannes for a critical reading of the manuscript. This research was supported by NIEHS grants R01 ES010807, and by the NIEHS Center for Environmental Genetics grant P30 ES06096. A.V.H. is supported by the NIEHS Training Grant T32 ES007250.

Additional information

Funding

This research was supported by NIEHS [grant number R01 ES010807], and by the NIEHS Center for Environmental Genetics [grant number P30 ES06096]. A.V.H. is supported by the NIEHS Training [grant number T32 ES007250].