Perturbation of the labile Zn²⁺ pool leads to changes in chromatin accessibility

To determine which TFs may be activated in Zn²⁺ replete or Zn²⁺ deficient conditions, we performed ATAC-seq on MCF10A cells subjected to perturbations of Zn²⁺. We first used a genetically encoded FRET sensor (NLS-ZapCV2) (Fiedler et al, 2017) that is specific for measuring labile (freely exchangeable) Zn²⁺ in the nucleus. Addition of 30 μM ZnCl₂ to the extracellular media caused a relatively slow but measurable rise in labile Zn²⁺ (Fig 1A–C). To limit any secondary effects that may occur from downstream activation of TFs, we opted to collect cells for ATAC-seq after 30 min of ZnCl₂ treatment. In situ calibration of the sensor combined with the binding parameters (K_D’ = 5.3 nM, n = 0.29 [Sanford et al, 2019]) revealed that at 30 min there was a change in labile Zn²⁺ from 150 pM to 75 nM. Addition of 50 μM of the Zn²⁺ chelator tris(2-pyridylmethyl)amine (TPA) resulted in a rapid decrease in nuclear Zn²⁺ to less than 1 pM (the minimum concentration of labile Zn²⁺ that can be quantified using NLS-ZapCV2) (Fig 1D–F). To keep treatments consistent, we also treated cells for 30 min with 50 μM TPA.

We performed ATAC-seq to assess which regions in the genome are differentially accessible in Zn²⁺ deficient and Zn²⁺ replete states. Treatment with TPA or ZnCl₂ induced broad changes in global chromatin accessibility (Fig 2A). TPA treatment impacted accessibility more than ZnCl₂ treatment (813 versus 517 differentially accessible regions with TPA versus ZnCl₂) (Fig 2B, SdataF2). Fig 2C shows genomic tracks for the top hit for TPA and ZnCl₂ treated cells. In both treatments, there was an overall trend towards reduced accessibility, with 593 regions and 334 regions having negative log₂FoldChanges, corresponding to 73% and 64% of the differentially accessible regions for TPA and ZnCl₂, respectively. Interestingly, whereas treatment with TPA led to broad changes in accessibility, suggestive of non-specific global decreases in accessibility, elevated Zn²⁺ led to a smaller number of highly significant changes in accessibility in select regions, suggesting that perhaps elevation of Zn²⁺ leads to changes in specific regions of chromatin (Fig 2A). One of the top hits showing increased accessibility with ZnCl₂ treatment was a region (chr16:56623012–56628480; log₂FC = 1.14, P_adj = 4.47 × 10⁻⁷) directly aligned with the metallothionein isoform MT1E, a gene regulated by the MTF1 in response to increased levels of Zn²⁺ (Fig 3). In addition, the region overlapping a second metallothionein isoform, MT2A, also showed increased accessibility though to a lesser extent (chr16:56605751–56611835; log₂FoldChange = 0.611, P_adj = 0.018). These results show that perturbations of the labile Zn²⁺ pool can alter the landscape of accessible chromatin in a short (30 min) time frame, with Zn²⁺ depletion leading to global decreases in chromatin accessibility and Zn²⁺ increases leading to more specific increases in accessibility. Furthermore, we observe changes in accessibility in genomic regions corresponding to key zinc regulatory genes (MT1E and MT2A), as would be expected upon perturbation of metal homeostasis.

TFEA reveals differential enrichment of TF motifs with open chromatin

To determine which TFs may be activated or repressed in each treatment, we performed TFEA (Rubin et al, 2021). Briefly, regions of open chromatin were annotated as peaks using HMMRATAC (Tarbell & Liu, 2019) and peaks from two technical replicates were combined using the muMerge algorithm (Rubin et al, 2021). Then, each of 1279 TF motifs from HOCOMOCO (Kulakovskiy et al, 2018) were subjected to motif scanning using Find Individual Motif Occurrences (FIMO) to identify potential TF-binding motifs within the detected peaks (Grant et al, 2011), and differential expression analysis using DESeq2 was performed to generate a ranked list of peaks that were differentially accessible between control and treatment conditions. This was used to generate an enrichment score (E-score) to determine which TF motifs were differentially enriched in regions of newly open chromatin between the control and treatment conditions. For this analysis we separated promoter regions, defined as 1,000 bp upstream of transcription start sites from distal regulatory regions, defined as regions more than 1,000 bp from a transcription start site.

Perturbation of Zn²⁺ led to many changes in the enrichment of TF motifs in regions of open chromatin. We found 685 motifs were differentially enriched with a P_adj < 1 × 10⁻⁷, corresponding to 507 unique TFs, 434 of which were zinc-finger TFs. The pattern of perturbation and the types of TFs that were most significantly affected were notably different at promoters compared with distal regulatory regions. Fig 4 shows the enrichment score in ZnCl₂ (elevation of labile Zn²⁺ to 75 nM) compared with the enrichment score in TPA (depletion of labile Zn²⁺ to ∼1 pM) for individual TFs in either promoter regions (Fig 4A, SdataF4.1) or intergenic regulatory regions (Fig 4B, SdataF4.2). Strikingly, in distal regulatory regions, a large number of TFs show significant motif enrichment when chromatin becomes more open in response to elevated Zn²⁺ (211 TFs, 187 of which are zinc-finger TFs). 83 of these TFs show motif depletion in regions of open chromatin in low zinc, suggesting reciprocal activation of these TFs in response to Zn²⁺ perturbation. In contrast to the intergenic region, at promoters the changes are dominated by TFs whose motifs are enriched in regions of open chromatin upon Zn²⁺ depletion. 268 of the 418 TFs with motifs that are differentially enriched in promoter regions with P_adj < 1 × 10⁻⁷ show enrichment in TPA. 158 of these TFs are also depleted in high Zn²⁺, suggesting reciprocal activation upon zinc perturbation. Notably, zinc-finger TFs (ZF TFs) are under-represented in this group (101 of 268 or 38%, compared with the 48% of the human TFs that are ZF TFs), suggesting that differential enrichment of these TF motifs may result from global changes in chromatin organization in low Zn²⁺. Fig 4C summarizes the differential enrichment results at promoters and intergenic regions, showing that a given Zn²⁺ perturbation tends to have an opposite effect in each region (e.g., elevation of Zn²⁺ causes enrichment of motifs with open chromatin in intergenic regions and depletion of such motifs at promoters).

Elevation of zinc has the biggest effect on chromatin around TF motifs at distal regulatory regions whereas depletion of zinc has the biggest effect on chromatin around TF motifs at promoters. There were a significant number of motifs (and unique TFs) that were differentially enriched in both promoter regions and intergenic regions, however there was a notable shift in the relative significance of these motifs. Fig 4D shows the relative significance (P_adj ranked from high to low) for TF motifs that were differentially enriched at both promoters and intergenic regions in high Zn²⁺. The plot reveals major changes in the rank order of TF motifs that were enriched in high Zn²⁺. Fig 4E shows a similar comparison for overlapping TF motifs that were depleted in low Zn²⁺. These results reveal that perturbation of Zn²⁺ affects TF motif enrichment differently at promoters versus distal regulatory regions.

We noticed that a large number of the differentially enriched motifs corresponded to Krüppel-type zinc finger (ZNF) TFs. ZNF proteins contain a median of 12 zinc fingers, suggesting they bind relatively long and specific sequences of DNA. Still, they are largely uncharacterized. A comprehensive catalog of ZNF proteins found that one third of all ZNF TFs contain a Krüppel associated box (KRAB) domain which recruits histone deacetylase proteins giving rise to transcriptional repression (Huntley et al, 2006). Given that we found 86% of the differentially enriched TFs in intergenic regions were ZF TFs, we examined what fraction of these also contained KRAB domains. We found that 97 TFs (40%) contain KRAB domains, suggesting this class of TFs is slightly overrepresented in the repertoire of TFs that bind in distal regulatory regions predicted to be differentially accessible by zinc.

As noted above, many of the TF motifs that emerge from TFEA correspond to TFs that either directly bind Zn²⁺ or are known to be correlated with metal homeostasis. That elevation of Zn²⁺ led to enrichment of motifs for zinc-finger TFs in regions of open chromatin in intergenic regions suggests that the functionality of zinc-finger TFs may be sensitive to increases in the labile Zn²⁺ pool in the nM range. On the other hand, functional annotation of TF motifs that are enriched in open chromatin in TPA did not reveal any meaningful trends, suggesting that perhaps these conditions lead to non-specific global changes in chromatin accessibility. For the motifs that are enriched in elevated Zn²⁺ and depleted in Zn²⁺ deficiency, many correspond to poorly annotated Zn²⁺ finger proteins, such as ZNFs and ZSCANs. However, there were some hits that are well characterized and that were particularly intriguing. For example, the CCCTC-binding factor (CTCF) is an 11-zinc-finger protein involved in chromatin organization. Two CTCF motifs were globally enriched in elevated Zn²⁺ (E-score = 0.44 and 0.38, P_adj = 1 × 10⁻¹⁹⁷ and 1 × 10⁻¹⁷⁶) and depleted in TPA (E-score = −0.29 and −0.23, P_adj = 1 × 10⁻⁷¹ and 1 × 10⁻⁶⁴) in intergenic regions. Our laboratory recently used single molecule microscopy to show that CTCF senses low Zn²⁺ conditions and becomes significantly more mobile when Zn²⁺ is depleted, suggesting decreased association with chromatin (Damon et al, 2022). We also observed that a motif of the zinc-finger DNA methyltransferase DNMT1 was strongly enriched in elevated Zn²⁺ (E-score = 0.52, P_adj = 1 × 10⁻⁵⁶) and depleted in TPA (E-score = −0.36, P_adj = 1 × 10⁻²⁹) in intergenic regions. DNMT1 is not a canonical TF but methylates DNA during DNA synthesis, and our laboratory has previously shown that DNA synthesis is impaired in Zn²⁺ deficiency (Lo et al, 2020).

Correlation of TFEA enrichment with TF occupancy for the candidate TF p53

Whereas TFEA identifies TF motifs that become enriched in newly opened chromatin upon a given treatment, it does not directly reveal whether TF occupancy on DNA has changed. To probe how TFEA enrichments correlate with TF occupancy, we selected a candidate TF and used chromatin immunoprecipitation to pull down regions of DNA bound to the TF, followed by quantitative real-time PCR (ChIP-qPCR) to quantify how much of the TF is bound to a given target gene. Our two main goals were to determine whether differential enrichment identified by TFEA could be used to predict changes in TF occupancy and to evaluate a bioinformatic plus ChIP-qPCR workflow that could be used to validate TFEA predictions. As a TF target, we selected p53 (encoded by the TP53 gene) because of its differential motif enrichment (Fig 5A) and because there are multiple ChIP/ChIP-seq datasets indicating robust vetting of the antibody for ChIP.

Because ATAC-seq identified more than 50,000 peaks, we sought to narrow the number of possible p53 binding sites against which to perform ChIP-qPCR. Hence, we used a previously published p53 ChIP-seq dataset (GSM3378513) (Karsli Uzunbas et al, 2019) where MCF10A cells were treated with 5 μM Nutlin-3A for 6 h to activate p53. We reasoned that this ChIP-seq dataset would allow us to identify bona fide regions of the genome to which p53 binds. Furthermore, we intersected the ChIP-seq dataset with a PRO-seq dataset (GSE227931) where MCF10A cells were treated with Nutlin-3A (10 μM, 3 h) to identify regions of active transcription marked by the presence of bidirectional transcripts (Azofeifa & Dowell, 2017). This reduced the number of potential p53 binding sites from approximately 13,500 (as identified in the ChIP-seq dataset) to 2,164 (ChIP-seq intersected with PRO-seq). We then took this list and queried whether any of the validated p53 binding sites were differentially accessible in response to Zn²⁺ in our ATAC-seq dataset (Fig 5A). Differential accessibility analysis using DESeq2 revealed that with TPA treatment, 62 sites were differentially accessible (P_adj ≤ 0.1), with 51 of these being less accessible and 11 of these becoming more accessible (Fig 5B, Table S1, SdataF5). Conversely, treatment with ZnCl₂ resulted in six regions being differentially accessible, with five of these being less accessible and one being more accessible (Fig 5B, Table S2, SdataF5). No sites with significant changes in accessibility were shared between both treatments.

From this list of differentially accessible sites, we selected six p53 binding sites against which to perform ChIP-qPCR. As shown in Fig 6, the six target sites (ERGIC1, NFIB, SFN, EGR1, PLD5, LRIG3-DT, Fig 6A–F, respectively) were selected because they showed differential chromatin accessibility via ATAC-seq for at least one of the Zn²⁺ perturbations and were associated with active transcription in MCF10A cells in response to Nutlin-3A based on PRO-seq bidirectional tracks (Fig 6). We also included the well-established p53 target gene CDKN1A (p21) as a positive control for immunoprecipitation and a negative control primer set that would not be expected to show enrichment with and without the p53 antibody.

We performed ChIP-qPCR against the target regions using a p53 specific antibody. We incubated chromatin with Protein A/G beads to assess the overall background signal and calculated the %IP of samples incubated for each p53 target (Fig 7A) and for the positive and negative qPCR controls (Fig 7B). As expected, there was a significant increase in the %IP for the +p53 antibody samples compared with the −p53 antibody samples. The signal to noise ratio (SNR) was calculated by dividing the %IP of the +p53 antibody sample by the %IP of the beads-only (−p53) sample for each target (Fig 7C). Of our gene targets, CDKN1A (p21) had the strongest overall signal and SNR, as expected given that it is a high affinity p53 target. The negative control (Negative Control Primer Set 1) showed no difference between the plus and minus antibody condition (Fig 7B), further validating the ChIP conditions. Most gene targets had comparable signal and SNR across conditions (average SNR > 5), with ZnCl₂ samples often displaying the largest SNR value except for NFIB and PLD5, which had larger SNR values in the control (no zinc treatment) samples. EGR1 showed the lowest overall SNR across all conditions, with the TPA and control conditions displaying an average value lower than five (TPA: 3.52; Control: 3.24; ZnCl₂: 7.91).

The fold enrichment of p53 binding in our TPA and ZnCl₂ treatments relative to the control was calculated by dividing the %IP of the +p53 treated (TPA or ZnCl₂) ChIP samples by the +p53 of untreated (Control) ChIP samples for each target. The selected gene targets showed changes in p53 binding as a consequence of Zn²⁺ perturbation, where the fold enrichment values (indicative of increased or decreased binding) largely agreed with their respective changes in accessibility (for accessibility changes where P_adj < 0.1). Fig 7D shows the fold enrichment ratio for each of the six targets along with the ATAC-seq prediction of increased or decreased accessibility (>1 or < 1, respectively). For five of the six targets, the mean enrichment ratio generally agreed with the accessibility predicted by ATAC-seq. ERGIC1 was the one target gene where differential p53 binding was opposite the accessibility changes by ATAC-seq. In TPA, accessibility was increased (Fig 6A), whereas binding by p53 decreased (0.49 versus 1 in TPA versus control). Overall, our results indicate that Zn²⁺-induced changes in the enrichment of TF motifs in regions of open chromatin can be associated with changes in TF binding to target genes.

Molecular cloning

PiggyBac-NLS-ZapCV2 was generated by digesting PiggyBac-NES-ZapCV2 (Lo et al, 2020) with EcoRI and SalI. NLS-ZapCV2 was amplified from pcDNA(3.1)-NLS-ZapCV2 with EcoRI (5′ end) and SalI (3′ end) restriction sites using the following primers: ggaattGAATTCGCTTGGTACCGAGCATGCC and gaagcgGTCGAGCCACTGTGCTGGATATCTGCAGAA TTC. The NLS-ZapCV2 insert and PiggyBac-NES-ZapCV2 (Lo et al, 2020) were subsequently digested with EcoRI and SalI, and the NLS-ZapCV2 insert was ligated into the PiggyBac backbone.

Cell culture

Wildtype MCF10A cells (#CRL-10317; ATCC) were cultured at 5% CO₂ in DMEM/F12 (#11320033; Thermo Fisher Scientific) supplemented with 5% Horse Serum (#16050122; Thermo Fisher Scientific), 1% penicillin/streptomycin (#15070063; Thermo Fisher Scientific), 20 ng/ml EGF (#PHG0313; Thermo Fisher Scientific), 0.5 mg/ml hydrocortisone (#H0888; Sigma-Aldrich), 100 ng/ml Cholera toxin (#C8052; Sigma-Aldrich, and 10 μg/ml insulin (#12585014; Thermo Fisher Scientific). To generate MCF10A cells stably expressing PiggyBac-NLS-ZapCV2, ∼1.5 million wildtype MCF10A cells were incubated with 1 μg PiggyBac-NLS-ZapCV2 and 200 ng Super PiggyBac Transposase (#PB210PA-1; Systems Biosciences) and electroporated using a Neon Electroporation System (Thermo Fisher Scientific). Cells were plated into antibiotic-free media (the same composition as above, minus antibiotic) for 24 h and then underwent puromycin selection using 0.5 μg/ml puromycin for 7 d before experiments.

Live cell fluorescence microscopy

FRET sensor calibrations were performed on a Nikon widefield microscope equipped with a Ti-E perfect focus system, an iXon3 EMCCD camera (Andor), mercury arc lamp, Lambda 10-3 filter changer (Sutter Instruments), and a 20x Plan Apo air objective (Nikon). Filter settings used were: CFP (434/16 excitation, 458 dichroic, 470/24 emission), YFP FRET (434/16 excitation, 458 dichroic, 535/20 emission), and YFP (495/10 excitation, 515 dichroic, 535/20 emission). Camera exposure times of 300 ms, EM multiplier of 300, and 1 MHz readout speeds were consistent across all channels.

Quantification of nuclear Zn²⁺ levels

MCF10A cells stably expressing PiggyBac-NLS-ZapCV2 were equilibrated at 37°C, 0% CO₂ in phosphate-free HEPES-buffered Hanks balanced salt solution (HHBSS) containing 1.26 mM CaCl₂, 1.1 mM MgCl₂, 5.4 mM KCl, 20 mM HEPES, and 3 g/liter glucose for 30 min before imaging. Baseline resting FRET ratios (R_rest) were collected for 10 min, at which point ZnCl₂ was added to and final media concentration of 30 μM for approximately 50 min. Cells were then washed 3x in phosphate-free HHBSS to remove excess ZnCl₂ and 50 μM of tris(2-pyridylmethyl)amine (TPA) diluted in phosphate-free HHBSS was added to chelate Zn²⁺ and determine the sensor’s minimum FRET ratio (R_min). Cells were washed 3x in phosphate-free HHBSS (with no CaCl₂ or MgCl₂) to remove excess TPA, and a solution containing 119 nM buffered Zn²⁺ (Tsien & Pozzan, 1989), 1.5 μM pyrithione, and 0.001% wt/vol saponin was then added to determine the maximum FRET ratio of the sensor (R_max). Images were then processed in MATLAB R2017b using previously published analysis code (Lo et al, 2020; Anson et al, 2021). Briefly, nuclei were computationally detected using the fluorescence signal in the YFP FRET channel and the background-corrected FRET ratio was calculated by dividing the background-corrected acceptor (YFP FRET) signal by the background-corrected donor (CFP) signal. For each cell, average labile Zn²⁺ concentrations were calculated using the equation [Zn2+]=Kd(R−Rmin)(Rmax−R)1/n with a K_d = 5.3 nM and Hill coefficient (n) = 0.29.

ATAC-seq sample preparation and sequencing

∼50,000 MCF10A cells were pre-incubated in phosphate-free HHBSS for 30 min at 37°C and 0% CO₂, followed by treatment with either 50 μM TPA, 30 μM ZnCl₂, or phosphate-free HHBSS for 30 min. Nuclei isolation and transposition followed the OMNI ATAC-seq protocol as previously reported (Corces et al, 2017) using the Illumina Tagment DNA TDE1 Enzyme & Buffer Kit (#20034197; Illumina). Libraries were generated using Illumina Nextera CD Indexes (#20018708; Illumina) and sequenced in paired-end mode (2 × 75 bp) on an Illumina Nextseq 500 at the BioFrontiers Next Generation Sequencing Facility (CU Boulder). ATAC-seq data were collected in two biological replicates for three different conditions: TPA treatment, ZnCl₂ treatment, and no treatment (control). For differential accessibility analysis, the TPA and ZnCl₂ treatment conditions were compared to a common control dataset.

ATAC-seq data processing

Quality control of demultiplexed fastq files was performed using fastQC (v0.11.8) (Babraham Bioinformatics, 2022). Adapters were trimmed using BBDuk (BBDuk Guide, 2022) using the following flags: ktrim = r, qtrim = 10, k = 23, mink = 11, hdist = 1, maq = 10, minlen = 20, tpe, tbo. Resulting trimmed fastq files were mapped in a paired-end fashion with HISAT2 (Kim et al, 2019) (v2.1.0) to hg38. Resulting BAM files were indexed with samtools (v1.8) (Li et al, 2009). Regions of open chromatin were annotated using HMMRATAC (Tarbell & Liu, 2019); for downstream analyses (e.g., differential accessibility, TFEA, etc.), the entire open region, rather than the annotated summit, was used as an input. For assessing which regions of open chromatin were differentially accessible, peaks from all treatments were merged using muMerge (Rubin et al, 2021). Peaks that overlapped with blacklisted genomic regions (Amemiya et al, 2019) were removed using bedtools (v2.28.0). Mapped read counts were then counted using featureCounts in R (v4.0.5) and differential accessibility was determined using DESeq2 (Love et al, 2014).

For assessment of differential TF motif enrichment we used TFEA (Rubin et al, 2021). Within the TFEA pipeline, regions of interest (ROIs) were defined by merging HMMRATAC-called peaks from two biological replicates using muMerge. Read coverages within each ROI were then calculated with HTSeq (v. 0.13.5) and ROIs were ranked according to the significance of differential accessibility with DESeq2 (v. 1.26.0, R v. 3.6.1). For analysis of only promoter or nongenic regions, annotation files were generated from Refseq annotations. Promoters were defined as 1,000 base regions directly upstream from gene transcription start sites, and nongenic annotations were created using bedtools intersect (v. 2.28.0) and included any genomic region outside of a gene or promoter. Ranked differential region lists were then filtered with bedtools to overlap promoter or nongenic annotation files prior to TF motif scanning. TF motif scanning was then performed within the ±1,500 bp window surrounding each ROI center using FIMO (meme v. 5.0.3) FIMO (Grant et al, 2011) at a P-value cutoff of 10⁻⁵ and a curated set of TF position weight matrices (Lambert et al, 2018). TFEA then computed motif enrichment scores for each motif for each treatment relative to the control.

For a single region to be “statically significant” there has to be a very large change in that region. In TFEA, like Gene Set Enrichment Analysis (Subramanian et al, 2005), the assay depends on a ranked list of regions/genes. Therefore, it is not based on the set of significant regions, but on the direction all regions are shifted in differential expression.

PRO-seq and library preparation methods

MCF10a cells were grown up in DMEM/F12 media with supplements as described (Levandowski et al, 2021) to a confluency of 70% before treatment with Nutlin-3a (10 μM Nutlin-3a in DMSO, #S8059; Selleck) or vehicle (0.1% DMSO) for 3 h before the nuclei were harvested. 10 × 10⁶ nuclei were prepared for nuclear run-on and library preparation as described (Kwak et al, 2013) with modifications (Fant et al, 2020). Samples were sequenced on the Illumina NextSeq 500 platform (single-end × 75 bp). These data are publicly available on NCBI GEO at GSE227931.

Selection of p53 binding sites for ChIP-qPCR

To find putative p53 binding sites to probe using ChIP-qPCR, we downloaded the BED file containing ChIP-seq peaks from GSM3378513 (Karsli Uzunbas et al, 2019) (p53 ChIP-seq in MCF10A cells treated with Nutlin-3a) from the Cistrome Data Browser (ID #105311; CistromeDB). CistromeDB assesses ChIP-seq datasets on six quality control metrics (sequence quality, mapping quality, library complexity, ChIP enrichment, fraction of reads in peaks, and overlap with DNase hypersensitive sites), and this dataset passed all QC metrics; therefore, we used the BED file available from the CistromeDB website. To increase confidence in p53 ChiP-seq peaks that correlate with regions of productive transcription, we used bedtools intersect with the flag “−F 0.5” to only retain ChIP peaks that overlapped at least 50% with regions of productive transcription. This reduced the number of p53 ChIP peaks from 13,453 to 2,163. We then used the featureCounts module in R to count the number of ATAC-seq reads mapped to these peaks and performed differential accessibility analysis using DESeq2. The list of significant (P_adj ≤ 0.1) hits was then manually curated to select for regions that were proximal to or within genes undergoing bidirectional transcription using the Nutlin-3A PRO-seq dataset. Of these, we selected six candidates whose target genes were of interest and lent themselves to successful qPCR primer validation via amplification efficiency testing.

The six candidates were selected based on the following criteria: (1) they exhibited differential accessibility in at least one of the Zn²⁺ perturbations, (2) they overlapped a p53 binding site identified by ChIP-seq, (3) they correlated with regions of active transcription upon activation of p53 by nutlin 3a, (4) primers could be designed within the introns of the gene or at most ∼1,000 base pairs up/downstream of the gene. Within the p53 targets in this list is an intronic region of endoplasmic reticulum-Golgi intermediate compartment 1 (ERGIC1), which encodes a cycling membrane protein whose potential function is to transport cargo from the ER to the Golgi (Breuza et al, 2004). In TPA treated cells, the accessibility at this genomic region was slightly, but significantly increased, whereas treatment with ZnCl₂ did not significantly alter its accessibility. In addition, an intronic region of the Nuclear Factor I B (NFIB), a TF that regulates differentiation and proliferation of the central nervous system (Chen et al, 2017), showed significantly reduced accessibility with TPA. We also targeted a p53 binding site upstream of stratifin (SFN), which, when activated by p53 because of detected DNA damage, leads to cell-cycle arrest at the G2/M transition (Taylor & Stark, 2001; Arakawa et al, 2022). This region saw a significant increase in accessibility in TPA. A region that saw a large increase in accessibility under ZnCl₂ conditions was a p53 binding site downstream of the early growth response 1 (EGR1) gene, which encodes a TF containing three C2H2-type zinc fingers that can regulate p53 activity upon ionizing radiation exposure (Inoue & Fry, 2018). Cells treated with ZnCl₂ also displayed a slight but significant reduction in the accessibility at the intronic regions of PLD5 and LRIG3-DT. Phospholipase D family member 5 (PLD5) has been shown to be regulated by a tumor suppressor microRNA that is itself a client of p53 (Liu et al, 2021), whereas the leucine-rich repeats and immunoglobulin-like domains divergent transcript (LRIG3-DT) encodes an integral plasma membrane protein whose expression has been shown to positively correlate with the expression of p53 in cervical intraepithelial neoplasia (Lindström & Hellberg, 2014).

As a positive control for immunoprecipitation, we used previously published (Menendez et al, 2013) primers targeting the p53 response element within the promoter of CDKN1A (p21), a well-established p53 target (Macleod et al, 1995; He et al, 2005). For CDKN1A, the p53 binding site was not included in our DESeq2 analysis, as the primary region of active transcription within the p53 PRO-seq dataset was downstream of the p53 binding site we chose. Visually, however, the region appeared to be slightly more accessible in TPA. As a negative control for immunoprecipitation, the Human Negative Control Primer Set 1 (#71001; ActiveMotif) amplifies a 78 base pair fragment from a gene desert on human chromosome 12 and would not be expected to show enrichment when incubated with the p53 antibody. This primer set is recommended as a negative control for almost all transcription factors.

ChIP-qPCR

Wildtype MCF10A cells (2 × 15 cm dishes with approximately 12 million cells/dish) were treated as noted above for ATAC-seq for 30 min, followed by washing with room temperature PBS, crosslinking using 1% methanol-free formaldehyde (#28908; Thermo Fisher Scientific) for 10 min and quenching for 8 min using 125 mM glycine. Cells were washed using cold PBS and nuclei were isolated using a hypotonic buffer (10 mM Tris–HCl pH 7.5, 50 mM NaCl, 2 mM EDTA, 1% NP-40, 1 mM DTT, and 1x protease inhibitors [#78441; Thermo Fisher Scientific]). Nuclei were then lysed and subjected to chromatin fragmentation using a BioRuptor (Diagenode) sonicator (7 x 10 min pulses, high setting) resulting in a mean chromatin length of ∼200 bp, followed by a spin at 21,130g for 15 min at 4°C to pellet insoluble material. Chromatin was aliquoted, flash frozen in liquid nitrogen, and stored at −70°C.

To immunoprecipitate p53-bound chromatin, 220 μl of chromatin (from approximately 7 × 10⁶ MCF10A cells) were pre-cleared per condition by incubating with 30 μl of Protein A/G beads (#SC-2003; Santa Cruz. #B2323; Lot) for 2 h at 4°C on a rotator. Then, 1/10^th of the pre-cleared chromatin from each condition was saved as “Input” and flash frozen, and half of the remaining chromatin from each condition was rotated overnight at 4°C with 4.5 μg of anti-p53 antibody (#554293; BD Biosciences, #9046678; Lot). The other half of the remaining chromatin from each condition was also rotated overnight with no antibody. Next, chromatin samples (+/− antibody for TPA, Control, and ZnCl₂) were incubated for 2 h with Protein A/G beads (15 μl of beads per 100 μl of chromatin; blocked overnight with 0.5 mg/ml bovine serum albumin and 0.4 mg/ml yeast RNA). The beads were then washed (resuspension followed by centrifugation) twice with 400 μl of each a low salt buffer (20 mM Tris–HCl pH 8, 150 mM NaCl, 2 mM EDTA, 1% Triton X 100, 0.1% SDS), a high salt buffer (20 mM Tris–HCl pH 8, 500 mM NaCl, 2 mM EDTA, 1% Triton X-100, 0.1% SDS), a high salt buffer containing LiCl (20 mM Tris–HCl pH 8, 250 mM LiCl, 1. mM EDTA, 1% sodium deoxycholate, 1% NP-40), and finally with TE buffer (10 mM Tris–HCl pH 8, 1 mM EDTA). Beads were then resuspended in 200 μl of elution buffer (50 mM Tris–HCl pH 8, 10 mM EDTA, 1% SDS) and incubated for 1 h at 37°C with shaking. Frozen inputs were thawed and diluted to a total volume of 200 μl of elution buffer and incubated for 1 h at 37°C with shaking. To reverse crosslinks, the resulting supernatant from chromatin samples and inputs were incubated overnight in 200 mM NaCl at 65°C. Next, samples were diluted 1:1 with Milli Q water to lower the salt concentration and proteins were digested using 160 μg of Proteinase K at 55°C for 1 h. DNA was purified via phenol:chloroform:isoamyl alcohol (#15593031; Thermo Fisher Scientific) extraction and ethanol precipitation. DNA pellets of chromatin samples (+/− antibody for TPA, Control, and ZnCl₂) and of inputs were then resuspended in 27 or 100 μl of diethylpyrocarbonate (DEPC)-treated water, respectively.

For qPCR, 10-fold serial dilutions (1:10, 1:100, 1:1,000) of each input from each condition were prepared to allow for an absolute quantification analysis of the data. The qPCR reactions were prepared as follows: 600 nM forward primer, 600 nM reverse primer, 3 μl of DNA (sample, diluted input, or a no-treatment control [DEPC-water]), and 1X SYBR Select qPCR master mix (#4472918; Thermo Fisher Scientific), with a final volume of 21 μl. qPCR was then performed in technical duplicate (10 μl/reaction) on a Bio-Rad C1000 Touch thermal cycler using the recommended SYBR Select protocol and an annealing temperature of 57°C. C_T values between technical replicates were averaged and the %IP for each chromatin sample and primer set was calculated by generating a standard curve of the log quantities (based on dilution factor) of input chromatin versus their C_T values, performing a regression analysis to calculate the quantity of DNA in the sample, and finally accounting for all dilutions made throughout the ChIP protocol given the equation: %IP=(DNA quantity)*(Volume of ChIP eluate)(Volume of input eluate)*(Volume of chromatin used for ChIP)(Volume of chromatin used as input). All qPCR primers used in this study are listed in Table S3. The entirety of the ChIP-qPCR protocol was run on six biological replicates.

The signal:noise was calculated by dividing the %IP in the +Ab sample by the %IP in the −Ab sample for each target. The enrichment ratio was calculated by dividing the %IP in the +Ab sample for a given target by the %IP in the +Ab sample for the respective control sample. An inter-quartile statistical test was run on the enrichment ratio data to determine statistical outliers.