Enhanced transcriptional heterogeneity mediated by NF-κB super-enhancers

The transcription factor NF-κB, which plays an important role in cell fate determination, is involved in the activation of super-enhancers (SEs). However, the biological functions of the NF-κB SEs in gene control are not fully elucidated. We investigated the characteristics of NF-κB-mediated SE activity using fluorescence imaging of RelA, single-cell transcriptome and chromatin accessibility analyses in anti-IgM-stimulated B cells. The formation of cell stimulation-induced nuclear RelA foci was abolished in the presence of hexanediol, suggesting an underlying process of liquid-liquid phase separation. The gained SEs induced a switch-like expression and enhanced cell-to-cell variability in transcriptional response. These properties were correlated with the number of gained cis-regulatory interactions, while switch-like gene induction was associated with the number of NF-κB binding sites in SE. Our study suggests that NF-κB SEs have an important role in the transcriptional regulation of B cells possibly through liquid condensate formation consisting of macromolecular interactions.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 expression heterogeneity through the modulation of transcriptional burst due to histone acetylation deposition and the presence of the initiation Ser5p RNA polymerase II [21,22]. In embryonic stem cells, it was reported that elements causing noise in gene expression include promoter type, conflicting chromatin states, lack of active histone modification, and SE existence [23]. In contrast, a stochastic model of SE-mediated transcriptional regulation exhibit a low variation pattern in transcription noise than a typical enhancer (TE) does [24]. Therefore, it remains unclear whether SE can truly induce an enhanced transcriptional heterogeneity and, if it is the case, what is the underlying mechanism resulting in heterogeneity.
In this study, we performed single-cell fluorescence imaging to confirm the visualization of NF-κB (RelA)-mediated SE formation and assessed the biochemical properties and dynamics of RelA foci upon anti-IgM stimulation in chicken DT40 B cells. We further utilized scATACseq (single-cell assay for transposase accessible chromatin with sequencing) to investigate the changes in SE activity upon cell activation through chromatin accessibility, and predicted the cis-regulatory interactions through co-accessibility analysis [25]. We observed that the formation of NF-κB nuclear aggregates upon nuclear translocation in anti-IgM-stimulated cells exhibits switch-like dynamics and was sensitive to a BET bromodomain inhibitor JQ1 and an LLPS inhibitor, 1,6-hexanediol. Single-cell RNA and ATAC-sequencing analyses further revealed that fold-change expression of the target genes and the enhanced transcriptional heterogeneity are correlated with the number of gained cis-regulatory interactions after cell stimulation. Our study also confirmed that BCR-mediated switch-like gene expression is associated with the number of putative NF-κB binding sites in a SE region. The overall results suggest that these distinct quantitative transcriptional responses originate from two mechanisms: open chromatin capacity, which allows direct binding of NF-κB to DNA, and the higher-order biomolecular assembly leading to enhanced DNA contacts. These findings also suggest that the biological functions of NF-κB SE emerge from a combination of chromatin status and binding properties of NF-κB, in particular, possibly mediated by liquid condensate formation, to facilitate divergent cellular responses.

NF-κB proteins form SE-like nuclear foci
Initially, we observed the NF-κB foci formation in DT40 cells upon nuclear translocation of GFP-tagged RelA protein, a component of the NF-κB heterodimer, under a fluorescence microscope. We observed that nuclear NF-κB maximally formed foci 20-30 min after anti-IgM stimulation (Figs 1A and S1A and S1B). As expected, an anti-IgM dose-dependent effect was observed in the numbers of RelA foci per cell, presenting a bimodal distribution across doses (S1C Fig), which was observed in a previous study [7]. Fitting the median number of foci of the dose-response curve to the Hill function resulted in a Hill coefficient of N = 4.33, suggesting that the formation of foci proceeded cooperatively (Fig 1B). Since cooperative binding of transcription factors is a marker of SE-mediated gene regulation [17,26], our results suggest that NF-κB may also activate the target genes cooperatively through the formation of these foci.
To further investigate the properties of the RelA foci, we focused on the relationship between NF-κB and BRD4, which is known to bind to active SEs [27]. We generated DT40 cells co-expressing RelA-GFP and mKate2-BRD4S (see Materials and methods). These cells revealed co-localization of BRD4S and NF-κB upon anti-IgM stimulation in a time-dependent manner (Fig 1C and 1D). We further inhibited BRD4 activity by adding JQ1, a BET bromodomain inhibitor, 60 min before anti-IgM stimulation. The number of RelA foci increased to a maximum value at 20 min, where the difference between the control and JQ1 treated cells was

PLOS GENETICS
Single-cell transcriptome of NF-κB super-enhancers negligible, then significantly decreased afterwards (Fig 1E and 1F). At the same time, JQ1 inhibition disrupted BRD4S puncta formation (S1D Fig). This indicates that RelA foci formation itself is BRD4-independent, but foci maintenance is BRD4-dependent. This is consistent with a previous result indicating that BRD4 maintained active NF-κB through RelA binding [16]. Conversely, inhibition of IκB kinase (IKK) using IKK-16 prevented the formation of RelA foci altogether (Fig 1E and 1F), confirming that RelA foci formation is NF-κB pathway-dependent.
We further investigated the presence of LLPS-like properties in RelA foci. Earlier experimental approaches examined the presence of IDRs in RelA [28]. We used PONDR VLXT [29] to analyze the IDRs and observed that RelA had high disorder scores in 250-500 amino acid residues (54.1% of all residues), which was more or less comparable to BRD4 (85%) (S1E Fig). To examine whether the RelA foci exhibited LLPS-like properties, we treated cells with 1,6-hexanediol (1,6-HD), a compound known to promote the dissolution of liquid-like condensates [30]. Fluorescence imaging analysis (Figs 1G and S1F) presented that treatment with 1,6-HD after 20 min of anti-IgM stimulation dramatically reduced the number of foci and that washing the medium recovered the foci. Although our assessment of the foci was not definitive as previously reported for other SE regulated condensates in an isolated system [17,20], these results suggest that the RelA foci have similar properties to SE and may be responsible for anti-IgM-dependent SE formation and target gene expression.

Enhancement of heterogeneity in gene expression upon anti-IgM stimulation
To examine the above possibility, we next analyzed the transcriptional responses of DT40 single-cell populations (n = 453) after one hour exposure to various anti-IgM doses (0, 0.01, 0.1, 1, and 10 μg/mL; S1 Table) using RamDA-seq, a method for single-cell RNA sequencing (scRNA-seq) [31]. Using Seurat, we obtained two distinct cell clusters (Fig 2A) [32]. This suggests the presence of two discrete cell populations over the different doses of anti-IgM, such as those obtained in the imaging analysis (S1C Fig). We further performed a pseudo-time analysis to confirm that two cell populations were discrete (S2A and S2B Fig). We classified these clusters as activated (red cluster) and inactivated (blue cluster), as determined by their mean expression levels of representative NF-κB target genes, such as NFKBIA, CD83, and TNFAIP3 [33] (Fig 2B).
Prediction of cell activation through logistic regression analysis on the number of foci per cell revealed similar quantitative profiles in cell activation predicted from RNA sequence data, in which most cells were activated at 1 and 10 μg/mL and inactivated at 0 and 0.01 μg/mL doses of anti-IgM (Fig 2C and S2 Table). This suggests that the cell states determined from the scRNA-seq and imaging analysis were closely related.
Since the cells exhibited digital states, we investigated the change of gene expression variability across different concentrations of anti-IgM stimulation to confirm whether gene expression correlates with phenotypical differences. Therefore, we calculated the Fano factor of 1,337 differentially expressed genes (DEGs) that were upregulated in activated cells compared with inactivated cells from the scRNA-seq analysis at each dose point. The Fano factor (O 2 /μ), which is typically used to calculate changes in transcriptional bursting [34,35], provides a measure of deviation from the Poisson distribution where O 2 /μ = 1. Further, we clustered genes by heterogeneity dynamics by scaling Fano factor changes per gene to identify patterns of gene expression heterogeneity across concentrations (Fig 2D). We observed that representative NF-κB negative modulators, such as NFKBIA [2] and TNFAIP3 [36], present in cluster 2, have high heterogeneity at 0.1 μg/mL and thus, correlate with cellular heterogeneity which is supposed to be highest at 0.1 μg/mL ( Fig 2C). Interestingly, cluster 4 presents an increasing

PLOS GENETICS
Single-cell transcriptome of NF-κB super-enhancers heterogeneity across anti-IgM concentrations, having the highest heterogeneity at 10 μg/mL, where most cells are phenotypically similar (Fig 2C). This cluster was observed to contain B cells activation-related genes, such as CD83 and BIRC3. The pseudo-time analysis also revealed that CD83 has a higher gene expression heterogeneity upon cell activation than NFKBIA (S2C Fig). After further inspection, we observed that cluster 4 presents a significantly high Fano factor across anti-IgM concentrations than other clusters (S2D Fig). More importantly, genes in cluster 4 were enriched for immune cell activation-related functions (S2E Fig). From these findings, we selected NFKBIA and CD83 as representative genes for further analysis.
To investigate the relationship between the RelA foci numbers and target gene mRNA in the same cell upon anti-IgM stimulation, we quantified their cellular mRNA using smRNA--FISH (single-molecule RNA-FISH) (Fig 2E). There was a positive correlation between the RelA foci number and mRNA expression of both genes at all dose points (R = 0.49 for NFKBIA, R = 0.52 for CD83) (S3A and S3B Fig), and thus, NF-κB foci formation is regarded as an indication of the expression of these target genes. We also confirmed that the RelA foci numbers per cell across doses of the CD83 and NFKBIA smRNA-FISH samples were similar (S3C Fig), suggesting that technical differences between the samples were minimal. It should be noted that the correlation coefficient in this study represents only a rough estimate of the correlation as there was a difficulty in accurately counting the number of RelA foci due to the degradation during smRNA-FISH sample preparation and differences in the timing of gene expression and RelA translocation. Moreover, to investigate the relationship between the RelA foci and transcriptional activity, we performed intronic smRNA-FISH to visualize nascent transcripts of CD83 and NFKBIA and observed colocalization between CD83 and NFKBIA nascent transcripts and RelA foci ( Fig 2F). However, using this method, we only observed colocalization in a subset of cells. Therefore, we believe that more rigorous methods, such as realtime visualization of transcription, might be required to validate the effect of RelA foci on the target gene loci and their respective mRNA expression.
For these two genes, we visualized the gene expression across different doses of anti-IgM stimuli through both smRNA-FISH and RNA-seq ( Fig 2G). For CD83, we see a lower basal gene expression at 0 μg/mL anti-IgM than NFKBIA and an increasingly wide gene expression distribution across anti-IgM concentration (Figs 2G and S3C). In contrast, NFKBIA exhibited a bimodal gene expression at 0.1 μg/mL, while there was a uniformly high gene expression across cells at 1 and 10 μg/mL. To quantify these differences, we further calculated the Fano factor change between cells at doses of 10 and 0 μg/mL, where cells were predominantly activated and inactivated, respectively ( Fig 2H). We observed a significantly higher fold-change increase in the Fano factor associated with CD83 (RNA-seq: 2.7, RNA-FISH: 10.1) compared to NFKBIA (RNA-seq: 1.1, RNA-FISH: 1.7). In addition, the bimodality of NFKBIA expression was well reproduced using smRNA-FISH especially at an intermediate anti-IgM concentration (0.1 and 0.01 μg/mL), while CD83 expression presented a single broad peak across different concentrations (S4A Fig). This evidence suggests that the transcriptional regulation between these two genes upon cell activation is potentially controlled by different modes.
One may be concerned that the measurement of the Fano factor is affected by changes in gene expression levels. However, the coefficient variation (CV), which is O2/μ, inversely correlates with gene expression magnitude [37]. This is especially evident in genes such as CD83, where the initial expression is close to zero, and expression upon activation is long tailed (Fig  2H), the change in gene expression distribution is inevitably dependent upon mean expression (S4A Fig). We showed that the change in distribution of gene expression is better recapitulated using Fano factor than CV (S4A and S4B Fig). This is especially evident in genes, such as CD83, where the initial expression is close to zero, and expression upon activation is longtailed ( Fig 2H).

Gained SEs are responsible for the strong induction of B cell developmentrelated genes
To clarify the underlying transcriptional mechanism, we explored epigenetic regulation of NF-κB mediated SEs and their effects on transcription. Typically, SE is identified using rankordering of super-enhancers (ROSE) algorithm [13] using transcription factors, H3K27Ac ChIP-seq, or DNase-seq data [13,14,26]. In our research, we utilized ATAC-seq data to investigate a large stretch of enhancers similar to SEs since both DNase-seq and ATAC-seq fundamentally identify open chromatin [38] and the fold-change in H3K27Ac and ATAC signals are correlated in BCR-stimulated B cell SEs [8].
Stitching of ATAC peaks was performed for peaks within 5 kb of each other, shorter than the default setting (12.5 kb) since ATAC peaks are often broader than TF ChIP-seq peaks [39]. The ROSE algorithm was further implemented [13,40] to rank the enhancers of the cells with and without anti-IgM stimulation ( Fig 3A). Changes in chromatin accessibility between unstimulated and stimulated cells were acquired by obtaining merged peaks and calculating the fold-changes of the signals between both conditions (see Materials and methods). We observed that anti-IgM stimulation evoked a significant change in SE but not in TE ( Fig 3B). Changes in SE signals can be observed in representative genes CD83 and NFKBIA (Fig 3A and 3C), and both the genes presented one of the highest gene expression fold changes among SE-controlled genes ( Fig 3D). The ATAC signal fold-change more robustly correlates with gene expression fold-change in SE than TE (R = 0.29 for SE, R = 0.084 for TE) (Fig 3E). We confirmed that the SEs detected using ATAC-seq were comparable to the previous report using H3K27Ac ChIPseq of mouse primary B cells [8], where SEs have a longer genomic length than TEs (S5A Fig ChIP-seq data [8]. Therefore, we concluded that the SEs identified from ATAC-seq are comparable to the SEs identified from H3K27Ac ChIP-seq data in this experimental setting. Interestingly, NFKBIA was classified as a SE target in control ( Fig 3A) and presented fewer accessibility differences upon stimulation (Fig 3C and 3E), indicating that NFKBIA expression might be less influenced by SEs and chromatin status than CD83.
We further classified the enhancers as gained, lost, or unchanged by the quantiles for TEs and SEs ( Fig 3F). In SE, the mean RNA fold-change is greater in higher dosage in gained SE and lost SE, while in TE, the difference is negligible (Fig 3G). This result confirms that large chromatin structure changes are associated with dynamic gene regulation in SEs.
To investigate the role of SE in the whole transcriptome, we took the intersection of DEGs and SE and TE ( Fig 3H). Among all DEGs, 16.7% and 48.4% of DEGs are controlled by SE and TE, respectively, and about half of SE and TE are involved in the regulation of DEGs. Expectedly, gain and loss of SEs presented a higher upregulation and downregulation of associated DEGs, respectively, compared to those of TEs ( Fig 3I). Further, examination of biological functions of the gained SE-associated genes that have enhanced expression exhibited the enrichment of functions related to immune cell activation and regulation (S5E Fig), confirming that SE-activated genes are responsible for B cell activation.

SEs are also defined by threshold gene expression tied to NF-κB binding
NF-κB and PU.1 coexistence have been suggested to be one of the defining features of SEs in B cell activation [8]. In this study, we demonstrated that there are higher NF-κB (median motif occurrences: 16 for SE, 0 for TE; median motif occurrences per 10 4 bp: 3.4 for SE, 0.0 for TE) Interestingly, between the enhancer categorization to gained, unchanged, and lost, we observed that there is a difference in NF-κB motif density (4.7 for gained, 3.2 for unchanged, 2.7 for lost), but not in PU.1 (10.7 for gained, 10.8 for unchanged, 11.3 for lost) at SE (Fig 4A and 4B). Therefore, this result updates the previous study [8] where co-localization of NF-κB and PU.1 are relevant to SEs, albeit NF-κB has a stronger effect in determining changes in chromatin accessibility and SE formation after B cell activation.
SEs have exhibited switch-like gene expression defined by high dose-response Hill coefficient (N) [8,24]. We investigated this property by obtaining the Hill coefficient of gene expression across the anti-IgM concentration of our data. We performed optimization to the Hill function for the genes (0.2 < fold-change) (S7 Fig) and further categorized the genes based on the Hill coefficient (see Methods) ( Fig 4C). We further analyzed the motifs at the enhancers of these genes categorized by the Hill coefficient ( Fig 4D). The threshold gene expression in B cell has been reported to be tied to the number of NF-κB binding through RelA ChIP-seq and motif analysis [8]. Interestingly, we observed that motif density of NF-κB significantly increased across Hill coefficient in gained SE (6.3 for high, 4.6 for mid, 4.1 for low), albeit not TE. In contrast, PU.1 motif density at enhancers did not present significant changes across different Hill coefficients (12.7 for high, 12.3 for mid, and 11.2 for low). CD83 and NFKBIA con-

Gene expression heterogeneity is potentially induced by gained cisregulatory interactions
Further, we focused on the TEs and SEs in DEGs between activated and inactivated cells to examine the transcriptional heterogeneity. We observed that SE-associated genes have a higher Fano factor at all anti-IgM concentration ranges ( Fig 5A) and that this trend is proportional to the stimulation dose (Fig 5B), as compared with TE. Further, there is a significant change in the Fano factor increase across different enhancer categorizations, where gained SE has a significantly higher Fano factor (S9A Fig). Fano factor change between 10 μg/mL and 0 μg/mL correlates with the increase in ATAC signal (R = 0.49 for SE, R = 0.19 for TE) (Fig 5C). Therefore, an increase in chromatin accessibility at enhancers is followed by an increase in heterogeneity in gene expression. However, it must be noted that the heterogeneity measured by the Fano factor, is also mean-dependent. In contrast, the classification of genes into high and low Fano factors does not agree with Hill coefficient classification (S9B Fig).
Due to the complex arrangement of high-order chromatin structures involving SE-mediated gene expression, it remains difficult to generalize the contribution of chromatin accessibility of SEs to transcriptional regulation [41]. Therefore, we utilized single-cell ATAC-seq data for co-accessibility analysis using Cicero [25] to elucidate the putative chromatin contact between enhancers and genes. We predicted accessible genomic regions that are potentially in close physical proximity by determining their co-accessibility scores and changes upon 10 μg/ mL anti-IgM stimulation. Co-accessibility calculations were performed separately for the stimulated and unstimulated samples. The peaks used for co-accessibility analysis were obtained by merging shortly stitched peaks (108 bp) from both samples. Regions with 0.1 � co-accessibility scores in the stimulated samples and 0.05 � differences between the stimulated and unstimulated samples were retained to detect significant increases in genomic interactions. To determine whether co-accessibility correlates with gene expression amplitude and heterogeneity, we annotated each gene with co-accessible pairs above the threshold and observed a positive correlation between the number of regions with co-accessible pairs and their RNA fold changes (Fig 5D) or Fano factor ratios ( Fig 5E). Here, we use the Fano factor change confirming its agreement with CV within DEGs (R = 0.71) (S9C Fig). RNA fold-changes (R = 0.4 for gained SEs, R = 0.15 for gained TEs) and Fano factor (R = 0.46 for gained SE, R = 0.13 for gained TE) positively correlated with the number of regions with gained co-accessible pairs for SEs. These results suggest that the increased number of cis-regulatory genomic interaction is a common mechanism for enhanced gene expression and transcriptional heterogeneity.
To closely examine CD83 and NFKBIA, we visualized the co-accessible region pairs within thresholds between peaks residing ± 1 kb of each gene annotated transcriptional start site (TSS) and other positions in the same chromosome (Figs 5F and S10A). Surprisingly, no coaccessible regions within annotated SE paired with the peaks ± 1 kb of the NFKBIA-annotated

PLOS GENETICS
Single-cell transcriptome of NF-κB super-enhancers

PLOS GENETICS
Single-cell transcriptome of NF-κB super-enhancers TSS were observed above the 0.05 score threshold, while several regions associated with CD83 were observed (Fig 5F). Stronger and more long-range interactions were gained for CD83 compared to NFKBIA (S10A Fig). These results suggest the increase in the formation of chromatin interactions upon cell activation for CD83 than NFKBIA, contributing to more enhanced expression heterogeneity.
Finally, we performed time-course RT-qPCR analysis on these two genes in the presence of an IKK inhibitor and JQ1 to investigate their dependence on the NF-κB pathway or requirement of BRD4 for their SE formation. CD83 expression pattern revealed a sustained pattern, while NFKBIA expression revealed a transient pattern on both qPCR and microarray data [42] (Figs 5G and S10B), suggesting that the differences in transcriptional regulation between these genes are potentially linked with their biologically different functions (NFKBIA as a feedback regulator of signaling pathway vs. CD83 as a B cell differentiation marker). Further, we confirmed that both genes were NF-κB pathway-regulated since the IKK-16 treatment suppressed the expression of both genes (Fig 5H, left). In contrast, JQ1 treatment significantly suppressed CD83 expression proportionally with higher JQ1 concentration, while NFKBIA expression slightly increased (Fig 5H, right), suggesting that CD83 gene expression was more highly controlled by BRD4, but NFKBIA was not. While we speculate that this difference is caused secondary to the initial accessibility in the enhancer region of NFKBIA in the absence of stimuli ( Fig 3A) and that the chromatin accessibility change was comparably small upon stimulation (Fig 3C), there might be other unknown mechanisms causing the increase in gene expression, which have been previously reported [43].

Discussion
In this paper, we aimed to reveal the SE activation and transcriptional regulatory mechanism controlled by NF-κB transcription factor using single-cell gene expression and chromatin accessibility analyses, supplemented by fluorescence imaging to characterize the biochemical properties of NF-κB foci. We used DT40 B cells, which were derived from chicken bursal lymphoma. Previous reports have revealed that the BCR-NF-κB activation pathway and its dynamics in mouse primary B cells are reproducible in DT40 cells [7,44]. We have also demonstrated a consensus in SE classification for target genes such as CD83 and NFKBIA, which are expressed at a slightly different time after anti-IgM stimulation (S10B Fig), in these analyses regardless of the difference in analytical approaches (H3K27ac ChIP-seq vs. ATAC-seq) (S5C Fig).
Initially, we characterized the properties of RelA foci through fluorescence imaging analysis and observed certain properties, such as sensitivity to LLPS perturbation, BRD4 dependence, and positive cooperativity. These analyses thus revealed that NF-κB nuclear aggregates exhibited properties similar to those of well-known SE condensates [19,20,24,45]. We conducted these analyses on DT40 chicken B cells, which are cultured at a higher temperature condition than murine or human cells. Therefore, this culture condition may have contributed to the more favorable formation of the RelA foci, especially concerning the LLPS-like properties.
SE-mediated gene regulation has been reported to cause transcriptional heterogeneity or noise [23]. However, another report presented that transcriptional noise is reduced in SEmediated gene regulation [24]. Through experimental and computational approaches, we revealed that NF-κB-mediated SE formation strongly modulated transcriptional heterogeneity. We observed that genes with an increased heterogeneity upon increasing stimulation dose are enriched with cell-type-specific immune regulatory genes (S2E Fig), supporting a previous report stating that heterogeneity in gene expression is tied to biological functions and may be used by cells as a bet-hedging or a response distribution mechanism [46,47], where cells exhibit heterogeneity to enable a response to changing environment and also allowing dose-dependent fractional activation respectively. This was observed in CD83, a B cell activation marker, demonstrating the involvement of heterogeneity in B cell development.
As a mechanism to produce gene expression heterogeneity in phenotypically identical cells, we observed that co-accessibility, which is concordant with genomic contacts, is correlated to the Fano factor change, indicating that gene expression heterogeneity possibly stems from cisregulatory genome interactions. NF-κB activation has been reported to be tied to the increase in heterogeneity in some genes and is attributed to the accumulation of Ser5p RNA polymerase II [21], which has also been reported to accumulate at enhancer regions [48]. This accumulation of RNA polymerase II is proposed to assist in gene expression activation through enhancer-promoter contact [49]. Our results support these conclusions since co-accessibility or putative cis-regulatory interactions correlate to Fano factor changes. However, we also observed a weak correlation, which may be caused by the inequality of different enhancer elements, such as presenting additive properties reported in gene activation in Drosophila development [50]. The enrichment of the TATA motif has also been proposed to generate transcriptional heterogeneity [23]. However, we observed a higher occurrence of TATAbox in genes associated with lost SE (S11 Fig), which might have caused gene expression heterogeneity in unstimulated cells. This heterogeneity might be owing to the difference in Pol II loading intervals [51]; however, the noise associated with gained SE is possibly generated by the fluctuation of high-order biomolecular assembly, and thus, the source of heterogeneity in gene expression might be different. While SE can form phase-separated transcription hubs containing multiple enhancers and/or promoters that might enable the higher diffusion rate of active enhancers and consequently, higher possibility of genomic DNA interactions [52], the current analytical method to analyze SE ignores individual enhancer properties [53] and nonlinear interactions. Therefore, a dissection of each enhancer contact may be required to understand how these differences arise in the context of SE.
To answer this question, we tried to utilize a CRISPR-Cas9 system with double-flanking gRNA and enhanced specificity S. pyogenes Cas9 [54] to induce long deletion in regions with the highest co-accessibility score change in the CD83 promoter ( Fig 5F). Unfortunately, we could not confirm the deletion of the regions within CD83 SE, which may be attributable to the possibly incomplete genome database of the chicken reference genome used for the gRNA design (GRCg6a). We induced a long deletion in the enhancer region outside of the CD83 SE. However, the obtained cell clone demonstrated low survival and impaired proliferation rates. Therefore, we could not confirm the changes in gene expression dynamics upon deletion of this enhancer region.
In contrast, SE has been reported to induce both switch-like gene expression [8] and enhanced heterogeneity [23]. A mathematical model has been used to explain that the threshold response is generated by the number of multivalent interactions mediated by transcriptional assembly within SE [24]. Therefore, genes showing threshold response might also show higher gene expression heterogeneity. However, our analysis revealed that this is not the case. Although we observed that there is a higher proportion of genes with Fano factor above the median in genes with higher Hill coefficient (1 < N < 5), there is no clear relationship between heterogeneity and threshold gene expression (S9B Fig). These results indicate that although enhanced gene expression and transcriptional heterogeneity are controlled by genomic contacts between chromatin accessible regions, switch-like expression of genes is controlled by NF-κB binding to SE (Fig 4), and those contributions are different for an individual gene. With this regard, we believe that the velocity and duration of nuclear translocation of NF-κB, along with interaction with other transcriptional factors, might strongly affect those transcription responses.
Lastly, we noticed that sequence-based identification of SEs may not truly reflect the biologically equal functions of SEs for each target gene. For example, NFKBIA and CD83 mRNAs exhibited varying sensitivities to JQ1 (Fig 5H, right). Since the NFKBIA SE had constitutively open chromatin prior to anti-IgM stimulation (Fig 3C), the signal-dependent effect on SE regulation of NFKBIA appeared to be small. These two genes demonstrated large fold-changes in gene expression upon SE regulation. However, the molecular level of SE regulation in each gene varied, which might have caused differences in transcriptional heterogeneity. Therefore, our results suggest that the functional evaluation of SE should not solely rely on sequence analysis but also on additional experimental methods.
In conclusion, these results provide an alternative insight into the mechanism of SE-mediated gene regulation upon NF-κB activation. Especially in the generation of biologically relevant gene expression, heterogeneity in phenotypically identical cell populations by SE is comprehended through gained cis-regulatory interactions. We speculate that this might be the one of the mechanisms for non-genetic heterogeneity in response to environmental factors. Thus, we expect similar study to be applied in further studies on other signaling pathways regulated by SEs.

mKate2-BRD4S transposon plasmid construction
We engineered the PB-TA-ERP2-mKate2-BRD4S construct from two addgene clones (65378 [55] and 80477 [56]) and pmKate2-H2B (Evrogen, Russia). Overlap extension PCR was performed to amplify the mKate2 insert from the mKate2-H2B plasmid while adding the attB1 adapter and linker sequence. Another round of overlap extension PCR was performed to amplify the BRD4S insert from GFP-BRD4 while adding the attB2 adapter and linker sequence. A final round of fusion PCR was performed to fuse the fragments containing mKate2 and BRD4 to create an insert containing mKate2-linker-BRD4S. The primers used for the overlap extension PCR are listed in S3 Table. A BP reaction using BP Clonase II (Invitrogen, USA) was performed to clone the insert into pDONR221 (Invitrogen), creating the entry clone pENTR221-mKate2-BRD4S. The entry clone was amplified using NEB-Stable (New England Biolabs, USA). An LR reaction using LR Clonase II (Invitrogen) was performed with the destination vector PB-TA-ERP2 [56] to create the final expression vector.

Cell transfection
Cell transfection was performed using the Neon Transfection System (Invitrogen). A total of 2 μg of plasmid was mixed with 1 × 10 6 cells in 10 μL buffer R. Electroporation was performed at 1400 V for 30 ms with 1 pulse. For piggyBac co-transfection, a mass ratio of 1:3 (PB:pBase) of the plasmid was used.

Quantification of RelA foci in single cells
Prior to stimulation and subsequent live imaging, DT40 cells were concentrated 5× into a serum-free medium and incubated in a 39˚C incubator with 5% CO 2 . Cells were further transferred into an eight-welled coverglass chamber (Iwaki, Japan). An inverted microscope IX81 (Olympus, Japan) equipped with a CSU-X1 confocal scanner unit (Yokogawa, Japan) and oilimmersion objective (100×, NA 1.45) was used to obtain images. MetaMorph software (Molecular Devices, USA) was used to obtain 13 z-slices of stack images at 1 μm increments in the zdirection. The image resolution used was 512 × 512 pixels (1 pixel = 0.16 μm). The observation chamber was maintained at 39˚C during observation. FIJI ImageJ 1.52i (https://imagej.net/ Fiji/Downloads) was used to count the foci using a custom macro, where the diameter for foci detection was set to 0.96 μm.

Quantification using smRNA-FISH
Fluorescence-conjugated chicken GAPDH, CD83, and NFKBIA probes were generated for smRNA-FISH using the Stellaris Probe Designer (https://www.biosearchtech.com/support/ tools/design-software/stellaris-probe-designer) (Biosearch Technologies, UK) according to the protocols of Biosearch Technologies. The sequences of the probes are presented in S4 Table. A total of 1 × 10 7 DT40 cells in 600 μL RPMI without supplements were stimulated with anti-IgM or PBS (control) for 30 min before washing with PBS and resuspending in 1 mL suspension buffer. Procedures following cell fixation were performed as described in the manuals of Biosearch Technologies for suspension cells. The fixed cells were mounted using Vectashield (Vector Laboratories, USA), sandwiched in cover glasses (Matsunami Glass, Japan), and sealed with clear nail polish prior to imaging.
A DeltaVision Elite-Olympus IX71 (Olympus) fluorescence microscope equipped with a Coolsnap HQ2 camera (Photometrics, USA) and an oil-immersion objective (60×, NA 1.42) was used for image acquisition. SoftWoRx software (Applied Precision, USA) was used for image acquisition and deconvolution. FISH-quant v3 was used to quantify the RNA-FISH images [57]. The default parameters were used in the quantification.

Quantitative RT-PCR (qRT-PCR) analysis
Total RNA was collected from the DT40 cells using a NucleoSpin RNA kit (Macherey-Nagel GmbH & Co., Germany) and subjected to complementary DNA synthesis and quantitative PCR using a ReverTra Ace qPCR RT Kit (Toyobo Life Science, Japan) and KOD SYBR qPCR kit (Toyobo Life Science) according to the manufacturer's protocol. PCR cycling conditions were as follows: 40 cycles of 10 s at 98˚C, 10 s at 60˚C, and 30 s at 68˚C. The primers used for qRT-PCR are listed in S5 Table. Expression values (n = 3) were normalized to those of GAPDH. For the inhibitor treatment, IKK-16 and JQ1 were added 60 min before anti-IgM stimulation to investigate the dependences of CD83 and NFKBIA on SEs.

Single-cell RNA-sequencing analysis
The DT40 cells were stimulated with anti-IgM (0, 0.01, 0.1, 1, and 10 μg/mL) for 1 h and sorted using an SH800 cell sorter (Sony, Japan) with a 130 μm sorting chip to select single cells. RamDA-seq with oligo-dT primer single-cell RNA-seq method was used for cDNA preparation [31]. The samples were sequenced using HiSeq2500 (Illumina, USA).
Seurat v3.2.1 was used for clustering and differential gene expression analysis of the scRNA-seq data [32]. Prior to clustering, a quality check of the data was performed to remove cell outliers (total count � 1.5 million, detected genes � 8500, and mitochondrial gene count ratio < 0.04), resulting in 89, 92, 87, 92, and 93 cells used in the analysis of the 0, 0.01, 0.1, 1, and 10 μg/mL anti-IgM concentrations, respectively. Data normalization was performed using Baynorm [60] and the log scaling method in Seurat. Data were regressed based on the cell cycle scoring with the "CellCycleScoring()" function of Seurat and mitochondrial gene count ratio. The top 2000 variable features were used for dimensionality reduction and clustering with a resolution of 0.2. DEGs were extracted using the "FindAllMarkers()" function of Seurat for genes with adjusted P-value < 0.05, this lenient categorization was done to retain genes that exhibit heterogeneity. Pseudo-time analysis was performed by creating a principal curve using the Princurve v2.1.4 R package on the dimensionality-reduced projection with the Lowess smoother.

Single-cell ATAC-sequencing analysis
All protocols for generating scATAC-seq data on the 10x Chromium platform (10x Genomics, USA), including sample preparation, library preparation, and instrument and sequencing settings, are described below available here: https://support.10xgenomics.com/single-cell-atac. Prior to nuclei extraction, the DT40 cells were stimulated with 10 μg/mL of anti-IgM or PBS for 60 min. Cellranger-atac count v1.1.0 (https://support.10xgenomics.com/single-cell-atac/ software/pipelines/latest/algorithms/overview) was used with default options for performing quality checks and mapping the scATAC-seq data to the genome. The sequence files were downsampled to 250 million reads before running the Cellranger-atac count pipeline. The reference genome used was the ENSEMBL genome with the annotation file GRCg6a.96 (ftp://ftp. ensembl.org/pub/release-96/gtf/gallus_gallus/).

Identification of SEs and TEs
Peak calling and enhancer identification from ATAC-seq data were performed using Homer v4.10.4 (http://homer.ucsd.edu/homer/) using the bam files generated from the Cell Ranger pipeline. Tag directories were created for the bam file from each condition using the "make-TagDirectory" program with the "-sspe -single -tbp 1" option. Peak calling was performed using the "findPeaks" program with the "-style super -typical -minDist 5000 -L 0 -fdr 0.0001" option. This procedure stitches peaks within 5 kb and ranks regions by their total normalized number reads and classifies TE and SE by a slope threshold of 1. Peak annotation was subsequently performed using the "annotatePeaks.pl" program with the GRCg6a.96 annotation file. The consequent peak files were merged between each stimulation condition for the SE and TE peaks separately using the "mergeBed" program of bedtools. Peak annotation was performed for the second time for the merged peaks to create the final SE and TE peaks. ATAC foldchange was then calculated for both conditions for the merged peaks separately for SE and TE. Genes associated with both SE and TE were assigned only to the SE.

ATAC-seq data visualization
Mouse primary splenic B-cell ATAC-seq data was obtained from the accession number DRA009931 of DNA Data Bank of Japan (DDBJ) [8]. Mapping of data to produce BigWig files was performed using nf-core/rnaseq pipeline [61] v3.3 (default options) and the default GRCm38 genome. Visualization of both primary and DT40 data was performed using Gviz v3.1.3 [62].

Gene ontology analysis
Gene ontology analysis was performed using the function "enrichGO" of clusterProfiler v3.14.3 for gained SE and RNA upregulated genes [63]. The ENSEMBL gene id was converted to the mouse homolog gene id prior to enrichment analysis using the query "mmusculus_homolog_ensembl_gene".

Co-accessibility analysis
Peak calling was performed for co-accessibility analysis in Homer v4.10.4 using the "find-Peaks" program with the "-style super -typical -minDist 0 -L 0 -fdr 0.0001" option to identify peaks constituting SE and TE. The "FeatureMatrix()" function of Signac v1.0.0 was used to assign fragments from the "fragments.tsv" file previously filtered for cell barcoding to the bed file containing peaks. Cicero v1.3.4.10 [25] was used to calculate the co-accessibility scores [25] between the ATAC peaks using the reference genome GRCg6a.96. The "max_sample_windows" argument of the "distance_parameters()" function was set to 1000, and the "max_elements" argument of the "generate_cicero_models()" function was set to 500. The other options were set to default. Co-accessibility calculations were performed separately for both the stimulated and unstimulated cells. The final co-accessibility scores were determined as the differences between the coaccessibility scores of the stimulated cells and unstimulated cells, where the initial score was �0.1 in the stimulated cells. Gviz v1.30.3 was used to visualize the co-accessibility between genomic regions.

Calculation of Hill coefficient
The focus dose-response curve was fitted to the Hill function with a constant additive term accounting for basal activity. The parameters were optimized using all methods of the "optimx ()" function in the "optimx" R package, and the method with the best fit was selected [65]. The equation used was as follows: where N, Km, k, and a represent the Hill coefficient (cooperativity), the binding affinity of NF-κB to the enhancer region, the rate constant for foci formation, and the constant additive term accounting for basal activity, respectively. For fitting of mean gene expression fold change to Hill coefficient in the Fig 4C and 4D, we performed optimization to the Hill function with removed basal activity term for genes with a fold change above 0.2 between 0 and 10 μg/mL anti-IgM and also removed genes with Hill coefficient above 9 and below 0.3 as these genes are either overfitted or do not fit at all. We then categorized genes with having high (5 < N), medium (1 < N < 5), and low (N < 1) Hill coefficients.

Clustering of Fano factor change
Hierarchical clustering of Fano factor changes across anti-IgM doses was performed using the "hclust()" function of R with the method option "ward.D2" and the tree was cut at 4 branches.

Correlation coefficients
Correlation coefficients were calculated using Spearman's rank correlation test to minimize the effects of outliers.

Statistics and reproducibility
The qPCR data are presented as the mean ± standard deviation (SD), wherein N = number of biological replicates. Data were evaluated using the unpaired Student's t-test. The means were considered statistically significant at p < 0.05. Box plots represent the median (centerline), interquartile range (IQR; box limits), and 1.5× IQR for the whiskers.  21. P-value was calculated using one-way ANOVA with undersampling (n = 21), n.s.: not significant. (B) Venn diagram and bar plot of the number of genes with assigned Hill coefficient (N) below 1, above 5 and in between and subsequent Fano factor (10 μg/ml) below (low) and above the median (high). (C) Correlation plot of Fano factor and CV ratio between 10 ug/ml and 0 ug/ml anti-IgM stimulation. Correlation coefficients were calculated using Spearman's rank correlation test.