Decoding the Epigenetics and Chromatin Loop Dynamics of Androgen Receptor-Mediated Transcription

Abstract Androgen receptor (AR)-mediated transcription plays a critical role in normal prostate development and prostate cancer growth. AR drives gene expression by binding to thousands of cis-regulatory elements (CRE) that loop to hundreds of target promoters. With multiple CREs interacting with a single promoter, it remains unclear how individual AR bound CREs contribute to gene expression. To characterize the involvement of these CREs, we investigated the AR-driven epigenetic and chromosomal chromatin looping changes. We collected a kinetic multi-omic dataset comprised of steady-state mRNA, chromatin accessibility, transcription factor binding, histone modifications, chromatin looping, and nascent RNA. Using an integrated regulatory network, we found that AR binding induces sequential changes in the epigenetic features at CREs, independent of gene expression. Further, we showed that binding of AR does not result in a substantial rewiring of chromatin loops, but instead increases the contact frequency of pre-existing loops to target promoters. Our results show that gene expression strongly correlates to the changes in contact frequency. We then proposed and experimentally validated an unbalanced multi-enhancer model where the impact on gene expression of AR-bound enhancers is heterogeneous, and is proportional to their contact frequency with target gene promoters. Overall, these findings provide new insight into AR-mediated gene expression upon acute androgen simulation and develop a mechanistic framework to investigate nuclear receptor mediated perturbations.


INTRODUCTION
The androgen receptor (AR) is a ligand-dependent transcription factor that plays a critical role in regulating gene expression in the prostate 1 .In its inactive form, the AR resides in the cytoplasm where it is stabilized by heat-shock chaperone proteins 2 .After binding androgens, such as testosterone or dihydrotestosterone (DHT), the AR undergoes an allosteric modi cation and translocates into the nucleus [2][3][4] .Once there, the AR binds to speci c cis-regulatory elements (CREs) on DNA through an interplay of chromatin accessibility, pioneer factors such as FOXA1, and sequence motifs [5][6][7][8][9] .The majority of these AR-bound CREs are proposed to function as enhancers as they are located distal from gene promoters 10 and are brought into close physical proximity by chromatin loops 11,12 .The enhancer activity at these CREs is associated with various epigenetic features, including chromatin accessibility, transcription factor binding, and post-translational histone modi cations, such as H3K27ac and H3K4me3 [13][14][15][16][17][18][19] .CREs are proposed to impact transcription through chromatin contacts with the target promoter that cause the recruitment of co-regulatory proteins and transcriptional machinery [20][21][22][23] .This typically involves multiple AR-bound enhancers and CREs that interact with the target promoter [24][25][26][27][28][29][30] .The contribution of individual CREs has been controversial with some studies suggesting that CREs work additively to increase gene expression [31][32][33] , while others propose that many enhancers are redundant and only contribute at speci c developmental stages [34][35][36][37] .The complex kinetic interplay of epigenetic modi cations, co-regulatory proteins and chromatin loops across multiple CREs following transcription factor binding is poorly understood.
To explore how AR binding impacts epigenetic modi cations and chromatin looping at regulatory elements in response to an acute perturbation, we generated an extensive multi-omics experimental dataset following androgen stimulation that was integrated into a graph-based framework.We demonstrated that AR binding sequentially induces an increase in FOXA1 and H3K27ac signals, that was followed by an increase in chromatin accessibility.We further show that AR does not induce new chromatin loops, but instead increases the contact frequency between gene promoters and selective ARbound enhancers.From these results, we generated and validated a multi-enhancer model, where a small subset of pre-established dominant CREs with increased chromatin contact frequency exhibits an elevated dynamic response to androgen stimulation, which signi cantly contribute to gene expression.These results provide a foundation for understanding how enhancers respond to an acute perturbation.
Assay for transposase-accessible chromatin with sequencing (ATAC-seq) LNCaP cells were isolated and subjected to modi ed ATAC-seq as previously described 38 .Brie y, 50,000 nuclei were pelleted and resuspended in ice-cold 50µl of lysis buffer (10 mM Tris-HCl, pH 7.5, 10mM NaCl, 3mM MgCl2, 0.1% NP-40, 0.1% Tween20, and 0.01% Digitonin).The subsequent centrifugation was performed at 500 g for 10 min at 4°C.The nuclei pellets were resuspended in 50µl of transposition buffer (25µl of 2× TD buffer, 22.5µl of distilled water, 2.5 µl of Illumina Tn5 transposase) and incubated at 37°C for 30 min with shaking at 1000 rpm for fragmentation.Transposed DNA was puri ed with the MinElute PCR Puri cation kit (Qiagen).DNA was puri ed using Qiagen MinElute (#28004), and the library was ampli ed up to the cycle number determined by 1/3rd maximal qPCR uorescence.ATAC-seq libraries were sequenced with 150-BP PE high-throughput sequencing on an Illumina HiSeq 2500 Sequencing platform (Novogene).
HiC combined with capture ChIP-seq (HiChIP) HiChIP was performed as previously described 39 .Trypsinized LNCaP cells (10 million) were xed with 1% formaldehyde at room temperature for 10 min and quenched.The sample was lysed in HiChIP lysis buffer and digested with MboI (NEB) for 4 h.After 1 h of biotin incorporation with biotin dATP, the sample was ligated with T4 DNA ligase for 4 h with rotation.Chromatin was sonicated using Covaris E220 (conditions: 140 PIP, 5% DF, 200 CB) to 300-800 bp in ChIP lysis buffer (1% NP-40, 0.5% sodium deoxycholate, 0.1% SDS and protease inhibitor in PBS) and centrifuged at 13,000 rpm.for 10 minutes at 4°C.Preclearing 30µl of Dynabeads protein A/G for 1 h at 4°C was followed by incubation with antibodies (H3K27ac, Diagenode, 3C15410196; H3K4me3, Diagenode 3C15410003).Reverse-crossed IP sample were pulled down with streptavidin C1 beads (Life Technologies), treated with Transposase (Illumina) and ampli ed with reasonable cycle numbers based on the qPCR with a 5-cycle pre-ampli ed library.The library was sequenced with 150-BP PE reads on the Illumina HiSeq 2500 Sequencing platform (Novogene).
Small Capped Nascent RNA Sequencing (Start-seq) For Start-seq, LNCaP cells were grown and collected as described as above.Cell pellets were washed with ice-cold 1x PBS. 1 million cells were treated with 1.5mL of Nun Buffer (0.3M NaCl, 1M Urea, 1% NP-40, 20mM HEPES pH 7.5, 7.5mM MgCl2, 0.2M EDTAm, protease inhibitors and 20U/mL SUPERase-IN) for 30 min on ice with frequent vortexing.Chromatin was pelleted by centrifugation at 12500 g for 30 min for 4°C.After 3 times-washing with 1mL ice-cold chromatin washing buffer (50mM Tris-HCL pH7.5 and 40U/mL SUPERase-In) and additional 0.5mL of Nun buffer.After centrifugation for 5 min at 500 g, 0.5m TRIzol was added to the remaining pellet.Libraries were prepared according to the TruSeq Small RNA Kit (Illumina).To normalize samples, 15 synthetic capped RNAs were spiked into the Trizolpreparation at a speci c quantity per 10^7 cells, as in 40 .The library was sequenced with 150-BP PE reads on the Illumina HiSeq 2500 Sequencing platform (Novogene).

CRISPRi
Guide RNAs (gRNAs) were designed for each enhancer and promoter region using CRISPR-SURF 41 .Cas-OFFinder was used to eliminate off-target gRNAs 42 .LNCaP cells stably expressing dCas9-KRAB (Addgene #89567) were seeded in a 6-well plate at a density of 200K cells per well.For transfection, a total of 500-1500ng DNA was used, divided according to the number of available gRNAs.Transfection was performed using Mirus TransIT-X2.After transfection, the media was replaced and 2ng/µl of puromycin was added for selection.Following 72 h of selection, the medium was changed to charcoal-stripped serum for androgen deprivation.After 48 h, the cells were treated with 10nM DHT for 4 h and then trypsinized.RNA extraction and cDNA preparation was performed using the LunaScript® RT SuperMix Kit.The androgeninduced expression was quanti ed using qRT-PCR and each experiment was conducted in triplicate.The sequences of gRNAs and qRT-PCR primers can be found in Supp Table 1.

ChIP-seq and ATAC-seq analyses
ChIP-seq and ATAC-seq were processed through the ChiLin pipeline 48 .Brie y, Illumina Casava1.7 software used for base calling and raw sequence quality and GC content were checked using FastQC (v0.10.1).The Burrows-Wheeler Aligner (BWA 49 , v0.7.10) was used to align the reads to the human genome hg19.Then, MACS2 50 (v2.1.0.20140616) was used for peak calling with an FDR q-value threshold of 0.01.Bed les and Bigwig les were generated using bedGraphToBigWig 51 and the union of narrow and broad peaks from ChIP-seq were used as anchors to call loops.The following quality metrics were assessed for each sample: (i) percentage of uniquely mapped reads, (ii) PCR bottleneck coe cient to identify potential over-ampli cation by PCR, (iii) FRiP (fraction of non-mitochondrial reads in peak regions), (iv) peak number, (v) number of peaks with 10-fold and 20-fold enrichment over the background, (vi) fragment size, (vii) the percentage of the merged peaks with promoter, enhancer, intron, or intergenic, and (viii) peak overlap with DNaseI hypersensitivity sites.For datasets with replicates, the replicate consistency was checked by two metrics: 1. Pearson correlation of reads across the genome using UCSC software wigCorrelate after normalizing signal to reads per million and 2. percentage of overlapping peaks in the replicates.

HiChIP Loop calling
After trimming adaptors from the HiChIP datasets using Trim Galore (v0.5.0) (https://github.com/FelixKrueger/TrimGalore),we used HiC-Pro (v3.1.0) 52as previously described in Giambartolomei and Seo et al. 39 .This aligned the reads to the hg19 human genome, assigned reads to MboI restriction fragments, and removed duplicate reads.We used the following options: <MIN_MAPQ = 20, BOWTIE2_GLOBAL_OPTIONS = -very-sensitive-end-to-end-reorder, BOWTIE2_LOCAL_OPTIONS =very-sensitive-end-to-end-reorder, GENOME_FRAGMENT = MboI_resfrag_hg19.bed,LIGATION_SITE = GATCGATC, LIGATION_SITE = "GATCGATC," BIN_SIZE = "5000.">All other default settings were used.To build the contact maps, the HiC-Pro pipeline selects only uniquely mapped valid read pairs involving two different restriction fragments.We applied FitHiChIP (v10.0) 53for bias-corrected peak and DNA loop calls.FitHiChIP models the genomic distance effect with a spline t, normalizes for coverage differences with regression, and computes statistical signi cance estimates for each pair of loci.We used the FitHiChIP loop signi cance model to determine whether interactions are signi cantly stronger than the random background interaction frequency.As anchors to call loops in the HiChIP analyses, we used 842367 regions for H3K27ac and 136939 regions for H3K4me3, which resulted from merging ChIP-seq narrow and broad peaks comprised 2 replicates for each broad and narrow peak and each of the ve-time points (0m, 30m, 4h, 16h, 72h) for H3K27ac, and 1 replicate for each broad and narrow peak and each of the same time points for H3K4me3.We used a 5 kb resolution and considered only interactions between 5 kb and 3 Mb.We used the peak to all for the foreground, meaning at least one anchor needed to be in the peak rather than both.The corresponding FitHiChIP options speci cation is < IntType = 3 > For the global background estimation of expected counts (and contact probabilities for each genomic distance), FitHiChIP can use either peak-to-peak (stringent) or peak-to-all (loose) loops for learning the background and spline tting.We speci ed the suggested option to merge interactions close to each other to represent a single interaction when their originating bins are closer.The corresponding FitHiChIP options speci cations are < UseP2PBackgrnd = 0 > and < MergeInt = 1> (FitHiChIP (L + M)).We used the default FitHiChIP q value < 0.01 to identify signi cant loops.We used hicpro2higlass (v3.1.0)to convert allValidPairs to .coolles after having speci ed the hg19 chromosome sizes, using the following command: <hicpro2higlass.sh-i sample.allValidPairs-r 5000 -c chrsizes -n > The reads from HiChIP data were merged from every time point for H3K27ac and H3K4me3 separately, and the reference loop sets were called with the same parameters above (resulting loop number, n = 296,326; n = 278,491 respectively).Next, the contact frequency values of each time point at loops are captured from the .coolles using the Python package, Cooler (v0.9.3).The count table was then TMM normalized among time points to reduce the between-sample variation.

Annotation of Cis-Regulatory Elements / De ning Background Genes
CREs were de ned as +/-2.5 kb around the summit of the accessible region from ATAC-seq peaks at any given time point.Promoters were de ned with a multi-step process.First, we identi ed the highest expression transcript (Gencode v19) 45 isoforms using Salmon (v0.14.1) 44 , see RNA-seq analysis.Next, the start locations -according to strand information-of the highest expressing transcripts were collected and the summit was extended +/2.5 kb.If they overlap with a de ned CRE, we assigned the overlapping CRE as the active promoters.From these annotated active promoters, the genes of HALLMARKS_ANDROGEN_RESPONSE from the Molecular Signature Database (MSigDB) 54 were de ned as AR-regulated genes.The transcriptome was divided into four quartiles and AR-regulated genes were compared with similarly expressing genes (Supp Fig. 1A, B).Overlapping CREs to an AR ChIP-seq peaks at any time point were de ned as potential AR-bound enhancers.The median number AR-bound enhancers (E + AR; ~2) of AR-regulated genes were less than AR-free enhancers (E-AR; ~14) (Supp Fig. 1E).Similarly, a CRE was considered FOXA1-bound, if it intersects with a FOXA1 peak at any time point (Supp Fig. 1C).To identify down-regulated genes, we calculated the log2foldChange induction between 16h and 0m, considering those with a value below − 1 as downregulated genes (Supp Fig. 3A, B).

Constructing the Graph Network
When we called signi cant loops from each time point separately, there were few called loops overlapping between time points potentially due to both the loop calling methodologies and experimental noise present in HiChIP data 55 .As we could observe the matching loops in the contact matrices of all time points (Supp Fig. 5), we instead generated a reference loop set, normalized the count matrix, and compared contact frequencies similar to published work 20 .A custom R (v4.1.1)script (https://github.com/lacklab/AR_transcription)was used to extract interaction pairs of annotated CRE regions (BED) according to the reference loop set (BEDPE; both H3K27ac and H3K4me3) using (v1.34) 56.The graph structure is built in custom Python (v3.9.7) script (link) using the NetworkX (v3.1) package 57 .Brie y, any pair of CREs are included as nodes in the network with TMM normalized HiChIP contact frequency as weighted edges, for both H3K27ac and H3K4me3 at every time point.

Epigenetic changes of CREs and kinetic changes of enhancers of AR-regulated genes
The average signal (AR, FOXA1, H3K27ac, H3K4me3, ATAC) at the promoters and CREs were collected using a custom python package (https://github.com/breambio/bluegill).To reduce the batch effect across time, TMM normalization was applied individually for each epigenetic factor 58 .To visualize the temporal change upon androgen stimulation, line plots were drawn (Fig. 2).For the selected six ARregulated genes (Fig. 3F), both the average signal (AR, FOXA1, H3K27ac, H3K4me3, ATAC) over the rstdegree interacting enhancers of each gene promoter and contact frequency (CF) between every gene promoter and corresponding enhancer were collected for every time point.The standard deviation (SD) of each feature along the time domain was calculated and min-max normalized.

GenomicInteractions
Promoter-centric analysis was done with the query loop sets, AR-regulated genes' promoter loops (P + AR), and random highly expressed genes' promoter loops (P-AR; n = 100 ; seed = 7).This was compared to reference loop sets (n = 100) that was randomized (1000 iterations; seed = 7) from highly expressed genes' promoter loops (Fig. 3B).Fold change was calculated for each iteration between the average contact frequency of loops in the query and reference.The same approach was used to calculate the contact frequency fold change from the enhancer viewpoint.The loop sets were selected according to the E-P pairs.While selecting the query and reference loops, the number of randomly selected promoters was xed (n = 100).The loops between AR-regulated gene promoters to AR-bound (E + AR) and AR-free enhancers (E-AR), and between AR-independent gene promoters to AR-free enhancers (E-AR) were compared to randomized loops from an AR-free enhancer to highly expressed gene promoters (1000 iterations; seed = 7).
The transcripts (Gencode v19) were extracted GTF le.The maximum signal within a range of +/-500 bp of TSS was gathered from the forward Start-seq track for plus-stranded transcripts, and the reverse for minus-stranded transcripts.Next, the log2foldchange (LFC), compared to 0m, was calculated for every time point (30m, 4h, 16h, 72h).If a transcript was found to be LFC > 1 at any time point, it was considered as differential up regulation.Next, the nascent expression levels at all time points for those transcripts (the union of differential up-regulation) were z-normalized to capture the highest expression time point, which is assigned to time-based expression groups (Fig. 4A).Similar to enhancer viewpoint analysis, the rst-degree AR-bound enhancer contact frequency to the gene promoters was compared to randomly selected contact frequencies (1000 iterations) between the highly expressed gene promoters (n = 100) and their rst-degree AR-free enhancers.

Contact frequency, expression correlation
To explore the contact frequency vs. expression, we performed three rst-degree summarization models (average, maximum, and sum).Brie y, we identi ed all rst-degree contacts of every gene promoter, then we summarized the contact frequencies with the corresponding function of the models for every gene.To reduce the signal noise, we rst sorted genes according to expression, binned them into equal-sized sets (k = 25), plotted bins as a scatter plot with average log expression on the x-axis, and averaged summarized contact frequency on the y-axis (Fig. 5B, Supp Fig. 7).

Random Forest Regressor
The binned (k = 25) contact frequency and expression data (see above) were used to train a random forest regressor from the Sklearn (v1.3.1) 59with "random_state = 0".The permutation importance of each feature is also calculated within the Sklearn 60 .

Gini Index Analysis
The chromatin contact frequency between gene promoter and rst-degree interacting elements was identi ed and the Gini index was calculated for a 16h time-point (Fig. 5C).To calculate the Gini index a custom python function was utilized (https://github.com/lacklab/AR_transcription).The mean absolute difference is the average absolute difference between all pairs of items in the population.The relative mean absolute difference is obtained by dividing the mean absolute difference by the population's average ( ) to account for differences in scale (e.q 1), where is contact frequency of a loop i, and there are n loops of a promoter.

Statistical analysis
The probability of nding k loops from each CREs was determined by hypergeometric distribution (Supp Fig. 4).When comparing the target set with the background set, a non-parametric Mann-Whitney-U test was applied.Correlation coe cients and their p-values were calculated with Spearman's test.All statistical tests were performed using the SciPy (v1.11.1)Python package 62 .

Generation of androgen-stimulation kinetic dataset and construction of regulatory networks
To characterize the temporal impact of AR binding on epigenetic features and chromatin looping, we generated an extensive kinetic multi-omics experimental dataset following androgen stimulation.We treated LNCaP cells with dihydrotestosterone and collected cells at 5 different time points (0m, 30m, 4h, 16h, 72h).At each time point, multiple features were characterized including gene expression (RNA-seq), chromatin accessibility (ATAC-seq), transcription factor binding and post-translational histone modi cation (AR, FOXA1, H3K27ac and H3K4me3 ChIP-seq), chromatin looping (HiChIP) and capped nascent RNA (Start-seq) (Fig. 1A).From these datasets, CREs (n = 78,522) were de ned from accessible sites (ATAC-seq) 15 .Based on known gene annotations and AR ChIP-seq, CREs were annotated as either promoters (n = 13,575), AR-free CREs (n = 58,265), or AR-bound CREs (n = 6,682).The interaction between these CREs was de ned from consensus H3K27ac (enhancer-centric; n = 296,326) and H3K4me3 (promoter-centric; n = 278,491; see methods) HiChIP loop sets.Supporting that these histone marks are associated with different functional CREs, only 23% of loops were found in both H3K27ac and H3K4me3 HiChIP (Fig. 1B).As expected, H3K27ac loops were predominantly between enhancer-enhancer (E-E) pairs, accounting for 43.2% of the loops, followed by enhancer-promoter (E-P) pairs, which constituted around 40.4% (Fig. 1C).In contrast, H3K4me3 loops were primarily between E-P pairs, making up 58.8% of the total loops.This is consistent with the known associations of these histone marks to promoter and enhancer CREs 63-65 .To allow more quantitative analyses of these large-scale chromatin interactions we transformed the resulting HiChIP looping data into a graphical network with each node representing a CRE and the edges being the chromatin loops between these two elements (Fig. 1D).With this regulatory network, we then overlaid the various multi-omics datasets to provide a framework that can interrogate the impact of AR binding.We investigated androgen-driven AR-mediated gene transcription based on the previously characterized hallmark androgen-responsive genes (Supp Fig. 1A, B).As expected, these were induced compared to similarly expressed random AR-independent genes (Fig. 1E).This comprehensive multi-omics dataset and graphical regulatory network provided a framework to quantitatively investigate the temporal impact of AR binding on epigenetic features, chromatin looping, and gene expression.

Androgen stimulation leads to stepwise epigenetic changes
With this structured regulatory network, we then explored how AR binding affects each epigenetic feature and its relationship with gene transcription.The AR-regulated genes' promoters (P + AR) and their looped AR-bound enhancers (E + AR) were compared to randomly matched from highly expressed ARindependent genes' promoters (P-AR) and their looped AR-free enhancer (E-AR) (Fig. 2A).We observed that E + AR displayed a strong AR and FOXA1 signal following androgen stimulation, that reached its peak at 4h and then subsequently decreased at 16h and 72h.Emphasizing FOXA1's pioneering activity, E + AR displayed elevated FOXA1 signal at the initial time point (0m).As expected, the AR and FOXA1 signal did not signi cantly change at AR-independent promoters or AR-free enhancer.Interestingly, those FOXA1-bound CREs that were not co-occupied by AR had no change in the relative FOXA1 signal (Supp Fig. 1C, D), suggesting that AR potentially in uences FOXA1 occupancy 66 .As expected, we observed a higher H3K4me3 ChIP-seq signal at promoters compared to enhancers.This signal was largely unaffected by AR binding, but there were selective genes, including KLK3, which exhibited an increasing H3K4me3 mark at its promoter following androgen treatment (Supp Fig. 2).We observed an elevated H3K27ac signal at promoters compared to enhancers (Fig. 2A).Further, the H3K27ac signal increased speci cally at E + AR, while it remained unchanged at E-AR.This change at E + AR was also observed for chromatin accessibility (ATAC-seq), though the maximum signal (16h) was found to occur after AR and FOXA1 peak occupancy (4h).We also investigated the epigenetic changes of AR down-regulated gene promoters and their AR-bound enhancers (Supp Fig. 3A, B).Interestingly, we observed that the AR-bound CREs (E + dAR), that were looped to the promoter of down-regulated genes (P + dAR) displayed a very similar increase in AR, FOXA1, H3K27ac, and chromatin accessibility (Fig. 2B).There was no statistically signi cant change at any time point in the AR-bound enhancers of either up or downregulated genes (p > 0.1).Overall, these kinetic datasets show that there is a sequential process that occurs following androgen stimulation where AR, FOXA1, and H3K27ac signals selectively increase at looped enhancers before inducing subsequent changes in chromatin accessibility.

AR-bound enhancers increase contact frequency to ARupregulated gene promoters
We rst investigated how chromatin looping changed following AR activation.Initially, we compared the number of loops formed following androgen stimulation.We found that both promoters of AR-regulated genes (P + AR) and their looped AR-bound enhancers (E + AR) did not exhibit any signi cant changes (p > 0.05) in the number of loops compared to background P-AR/E-AR (Fig. 3A).Supporting this, there was no signi cant difference in the distribution of the number of loops on gene promoters during androgen treatment (Supp Fig. 4A).These results demonstrate that AR binding does not cause a substantial rewiring of chromatin looping.
While androgen treatment did not signi cantly change the number of loops, we did observe an increase in contact frequency at AR-bound enhancers looped to AR-regulated genes (Supp Fig. 5).To quantify these changes, we calculated the fold change in contact frequency compared to bootstrapped (b = 1000) random AR-independent genes from both promoter and enhancer viewpoints (Fig. 3B).We observed that AR activation increased the contact frequency of loops at AR-regulated promoters (P + AR) over time, with a peak at 16h (H3K27ac: p < 0.001; H3K4me3: p < 0.001) (Fig. 3C).In contrast, the change in contact frequency of loops to promoters of AR independent genes (P-AR) remained stable.From an enhancer viewpoint (Fig. 3D), AR-bound enhancer CREs (E + AR) showed a signi cant increase in both H3K27ac and H3K4me3 contact frequency (H3K27ac: p < 0.001; H3K4me3: p < 0.001).In contrast, AR-free CREs (E-AR; grey) that were connected to the same AR-regulated gene promoters did not exhibit any signi cant change.Further, no change in contact frequency was observed in the CREs (E-ARi; black) that interact with AR-independent gene promoters.While chromatin contact frequency increased between AR-bound enhancer CREs (E + AR) and upregulated hallmark gene promoters (P + AR), we did not observe any change in loops between AR-bound CREs (E + dAR) and downregulated gene promoters (Fig. 3E).This is particularly striking as, AR-bound CREs (E + AR, E + dAR) exhibited a similar pattern in their epigenetic pro les, regardless of whether they looped to an upregulated or downregulated gene (Fig. 2A, Fig. 2B).
Given that a similar trend is observed at all AR-bound CREs, this suggests that epigenetic modi cations must be combined with chromatin looping to contribute to gene expression.
To interrogate the kinetic changes in epigenetics and contact frequency we focused on several upregulated androgen-responsive genes (KLK2, KLK3, NKX3-1, UAP1, ABCC4, and DHCR24) (Fig. 3F).As expected, AR-free enhancers (E-AR; grey) had minimal changes in epigenetics and contact frequency.Interestingly, we observed that AR-bound enhancers (E + AR; orange) had broad heterogeneity in response to treatment.The same enhancer typically had both the highest change in epigenetic features and contact frequency.Our ndings demonstrate that AR-bound enhancers increase contact frequency with AR-regulated gene promoters in response to androgen treatment.However, this response is heterogeneous and there is signi cant variability among AR-bound enhancers.

Association of nascent transcription to epigenetic changes and contact frequency
To determine if the change in contact frequency precedes or occurs simultaneously with AR-mediated gene transcription, we characterized the kinetic rate of androgen-induced gene expression.We identi ed up-regulated androgen-induced genes from nascent RNA (Start-seq) and categorized them based on the maximal expression (30m, 4h, 16h or 72h) (Fig. 4A).We chose to group these genes based on nascent RNA, as RNA-seq provides only steady-state quanti cation of mRNA transcripts 67 .Between these groups, no signi cant difference was observed in the maximal signal at AR-bound enhancers for AR, FOXA1 and H3K27ac ChIP-seq or chromatin accessibility (Fig. 4B), suggesting that the epigenetic features do not dictate the timing of maximal nascent RNA production.However, we observed a temporal relationship between nascent transcription and chromatin contact frequency.Speci cally, the maximal nascent transcription (4h and 16h groups) occurred simultaneously with an increase in H3K27ac and came before the peak of H3K4me3 (16h) maximal contact frequency at AR-bound enhancers (Fig. 4C).These observations indicate that the change in contact frequency does not precede RNA transcription.Overall, these ndings suggest that although AR and FOXA1 rapidly bind at enhancers upon androgen treatment, the maximal transcription occurs either simultaneously with or before a maximal alteration in chromatin contact frequency.

Multi-enhancers in uence transcription proportional to contact frequency
Having observed marked heterogeneity in chromatin loop contact frequency at AR-bound enhancers, we began to explore the impact of multi-enhancer contacts on gene expression.We expand the scale of our analysis across all detectable genes (n = 13,575), and characterized how loop frequency correlates to gene expression.Most genes had multiple CREs interacting with target promoters with a median frequency of ~ 15 interactions per gene (Supp.Figure 6A).With this large dataset, we evaluated three models correlating H3K27ac and H3K4me3 contact frequency to gene expression where: all interactions contribute equally to gene expression (average), only a single strong interaction affects gene expression (maximum) and all interactions contribute to gene expression in an additive manner (sum) (Fig. 5A).While all models strongly correlated H3K27ac and H3K4me3 contact frequency to gene expression (p < 0.0001; Spearman's test), we found that the sum model (H3K27ac; r = 0.896, H3K4me3; r = 0.904) and maximum model (H3K27ac; r = 0.839, H3K4me3; r = 0.876) correlated signi cantly better than the average model (H3K27ac; r = 0.643, H3K4me3; r = 0.679) across all time points (Fig. 5B and Supp Fig. 7).To quantify the importance of each feature, we built a random forest regressor to predict the expression from these operations (Accuracy; r 2 = 0.805) and found that the additive model (sum) best-predicted expression (Supp Fig. 6B) 68 .However, as we observed similar correlation with the maximum model, which scored only a single chromatin loop, this suggested that the contact frequency of chromatin loops to a promoter were markedly unbalanced and that there was a "dominant" interaction.To characterize this, we quanti ed the inequality of chromatin loop contact frequency to a target promoter and found that both H3K27ac and H3K4me3 were strongly unbalanced (Gini > ~ 0.5; Supp Fig. 6C).These ndings suggest that enhancers interacting with gene promoters do not have a uniform distribution in contact frequency, and that there are "dominant" loops that strongly correlate with expression (Fig. 5D).
To better understand how these potential dominant loops are dynamically affected by acute perturbation, we characterized the CRE-promoter interactions of androgen-regulated genes (n = 88; Supp Table 2).Dominant loops were identi ed for every gene promoter by rst scaling the contact frequency of interacting CREs in the range (0, 1), and selecting the highest ones with a threshold (0.8).The dominant loops were not solely based on proximity to the gene promoter as they were the closest CRE in only 27% of AR-regulated genes (n = 24) at any time point.Suggesting that the dominant loops are "primed" before AR binds, the same dominant loops were commonly found at every time point (Supp Fig. 8).Furthermore, we found that the dominant AR-bound CREs (D+) had signi cantly (p < 0.0001) higher dynamic change in contact frequency than non-dominant AR-bound CREs (D-) that interacts with AR-regulated gene promoters (Fig. Interestingly, the dominant loops were highly gene speci c.When we characterized two AR-regulated genes, KLK3 and KLK2, that share many looped CREs (Fig. 5E), we observed genespeci c disparities in contact frequency changes from an AR enhancer (ARBS3) to either the KLK2 or KLK3 promoters.While ARBS3 displayed the strongest and most dynamic contact frequency with the KLK2 promoter, it did not exhibit the strongest change for the KLK3 promoter.Instead, the loop with wellknown upstream AR-enhancer 69 (ARBS2) dominates others.Interestingly, both dominant loop enhancers were not the CRE with the highest change in AR peak height, suggesting that there are additional factors that contribute to changes in contact frequency.Overall, these results demonstrate that AR binding does not affect chromatin looping equally and that those CREs with the most dynamic contact frequency potentially have a greater impact on gene transcription.

CRISPR-based perturbations con rm the existence of "dominant" chromatin loops
To validate these correlative models, we tested the effects of dominant chromatin loops on androgen-induced gene expression.Utilizing CRISPRi, a derivative of the CRISPR/Cas9 system that suppresses enhancer activity without altering DNA sequences, we inhibited all AR-bound enhancers (n = 20) that interact with (Supp Table 1, Supp 9) the previously described AR-regulated genes (KLK2, KLK3, NKX3-1, UAP1, ABCC4, and DHCR24).Across all tested genes we found that inhibiting those CREs, that had a dominant chromatin loop, signi cantly impacted androgen-induced transcription (Fig. 6A).Importantly, as we observed highest dynamic change in contact frequency at ARBS3 and ARBS2 for KLK2 and KLK3, respectively (Fig. 5E), the highest gene expression inhibition is observed at these ARBSs.
Supporting the kinetics of the dominance model, we observed an inverse correlation between the gene expression inhibition and the dynamic change in contact frequency of loops during androgen stimulation (Fig. 6B).This correlation was greater than other genomic features including AR, FOXA1, H3K4me3, H3K27ac and chromatin accessibility (Supp Fig. 10).These results demonstrate that not all CREs contribute equally to androgen stimulation and that their contact frequency to promoters correlates with their impact on gene expression.Overall, these results underscore the importance of spatial genome organization in transcriptional regulation and validate our proposed dominance model.

DISCUSSION
Transcription factors bind to speci c DNA sites and regulate gene expression through the recruitment of co-regulatory proteins that activate transcriptional machinery 70,71 .Yet despite extensive research, many questions remain about how this complex process occurs, particularly as there are multiple CREs that interact with each promoter.Using the ligand-activated AR as a model system, we characterized how transcription factor activation changes the epigenetics and chromatin looping.Similar to published studies, our work showed that the AR stabilizes/recruits FOXA1 and increases both H3K27ac and chromatin accessibility (Fig. 2A) 2,5,7,72,73 .However, by characterizing a kinetic dataset, we showed that the change in chromatin accessibility occurs after both FOXA1/AR binding and H3K27ac posttranslational modi cations, suggesting that this is mediated by the recruitment of additional chromatinremodeling proteins.Those enhancers that are not bound by AR do not exhibit these changes.However, these epigenetic changes do not solely drive gene expression, as we observed a consistent pattern in all AR-bound CREs, regardless of either the direction of expression (upregulated and downregulated), or timing of maximal nascent transcription (Fig. 2; Fig. 4B).These observations suggest a sequential mechanism, independent of gene expression, in which AR activation recruits speci c co-regulatory proteins that alter histone modi cation and chromatin accessibility.
This work also characterized how AR binding affects genome organization.This is a controversial eld, with earlier studies proposing that steroid hormone receptors signi cantly reorganize chromatin looping when activated 74 .However, recent work has suggested that gene expression may occur through alreadyexisting interactions 23 .Our research found that androgen treatment does not cause rewiring of chromatin contact loops but instead increases the contact frequencies of previously established loops (Fig. 3C, D; Supp 5).Interestingly, we observed that maximal chromatin looping happens either during or after the maximal nascent gene expression (Fig. 4C).This suggests that increasing chromatin looping does not precede, but instead likely occurs simultaneously with gene transcription.Supporting this result, recent work observed that higher nascent RNA production is associated with a higher frequency of chromatin contacts 22 .Changes in the chromatin contact frequency have also been shown to be associated with differential gene expression 20,23 .Although several studies observed temporal changes in chromatin looping that occur before maximal RNA expression 21,23 , this distinction is likely due to the experimental methodology used, as RNA-seq primarily quanti es mature mRNA, whereas Start-seq captures only nascent mRNA.Highlighting the consistent pattern of epigenetic features of AR-bound enhancers (Fig. 2; Fig. 4B), can infer chromatin looping is an additional mechanism which regulates gene expression independently of AR binding.The importance of chromatin loop contact frequency is highlighted when we look at AR-downregulated genes, which show similar changes in epigenetic alterations and chromatin accessibility but not contact frequency (Fig. 2B; Fig. 3E).This suggests that AR binding recruits additional factors that potentially increase the contact frequencies of pre-established loops between enhancers and their target gene promoters to regulate the gene expression.
Numerous studies demonstrate that multiple enhancers contribute to the expression of a single gene .
However, there is no consensus about how each CRE contributes to gene expression.Our work suggests that connected enhancers contribute unevenly, and there exists dominant loops that have the largest impact on gene expression.These ndings align with the assumptions made in the activity-by-contact (ABC) model 76 , which presumes that the enhancer's impact on gene expression is correlated with the strength of the contact between them.Both our validation (Fig. 6), and a recent large-scale CRISPRi study from Barshard et al. 22 , demonstrate that those enhancers with higher contact frequency to target promoters are more likely to be functionally important which suggest a "dominance" model (Fig. 5D).
However, further improvements of this model can help us to understand the effect of individual CREs, thus allowing us to evaluate how multiple transcription factors impact gene expression.
Given the complex nature of this system, we had to make several assumptions.First, we chose to de ne individual regulatory units based on chromatin accessibility, since de ning enhancers and promoters by histone modi cations is prone to false positives 77,78 .Within these CREs, we focused on hallmark androgen response genes (Fig. 1E), as these have been extensively demonstrated to be regulated by AR.
Next, we used HiChIP, a protein-centric HiC method, rather than conventional HiC to provide enhanced resolution and speci city in capturing protein-bound chromatin 79 .By characterizing both promoter-centric (H3K4me3) and enhancer-centric (H3K27ac) chromatin loops, we believe this provided us with different genomic perspectives of chromatin loops and reduced potential biases (Fig. 1B, C) 21,80,81 .It is important to note that while H3K27ac is strongly impacted by AR binding we observed that H3K4me3 was relatively static (Fig. 2; Fig. 4B), suggesting that the change in contact frequency (Fig. 3C-E; Fig. 4C) following androgen treatment was not due to capture e ciency.Finally, as gene expression regulation circuits are challenging to integrate due to the intricate network structures formed by diverse CREs (Fig. 1D), we utilized regulatory networks wherein the nodes and edges represent the CREs and chromatin loops, respectively.Several studies have employed graph-based approaches to connect genomic regions and address various high-dimensional biological inquiries [82][83][84] .The versatility of this approach underscores the quantitative advantages of utilizing graphs to characterize chromatin loops.Despite these limitations, this study represents one of the largest experimental datasets (Fig. 1A) to characterize AR-mediated transcription.
Overall, this paper provides new insight into several important aspects of AR-mediated gene expression.
We show that AR binding triggers a temporal cascade which increases FOXA1 and H3K27ac that affects chromatin accessibility.Further, we demonstrate that AR does not introduce novel chromatin loops, but instead increases the contact frequency between AR-bound enhancers and their target promoters.
However, the effect of each enhancer on gene expression is markedly heterogeneous and proportional to promoter contact frequency.These ndings suggest that while AR binding to DNA induces a step-wise epigenetic alteration, the impact of bound enhancers is strongly dependent on the contact frequency of the established chromatin loops with the target promoter.Activation of the androgen receptor leads to a delayed increase in histone modi cations and chromatin accessibility (A) TMM normalized ChIP-seq or ATAC-seq signal were compared across all time points at different regulatory elements including: Promoters of AR up-regulated genes (P+AR: red), promoters of AR-independent genes (P-AR: black), AR-bound enhancers (E+AR: yellow) of AR-regulated genes, AR-free enhancers (E-AR: grey) of AR-independent genes.(B) A similar analysis was done for the promoters of AR down-regulated genes (P+dAR: blue), and their AR-bound enhancers (E+dAR: yellow).

Figure 3 AR
Figure 3

Figure 5 Multi
Figure 5