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

  1. Johannes N. Wibisana
  2. Takehiko Inaba
  3. Hisaaki Shinohara
  4. Noriko Yumoto
  5. Tetsutaro Hayashi
  6. Mana Umeda
  7. Masashi Ebisawa
  8. Itoshi Nikaido
  9. Yasushi Sako
  10. Mariko Okada

Reviewed by

Read the full article See related articles

Listed in

Log in to save this article

Abstract

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.

Article activity feed

  1. Version published to 10.1371/journal.pgen.1010235
  2. Review Commons

    Note: This rebuttal was posted by the corresponding author to Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Reply to the reviewers

    Reviewer 1

    Reviewer #1 (Evidence, reproducibility and clarity (Required)):

    *The manuscript by Wibisana et al. describes an impressive set of experiments that analyse the NFkB response at the single-cell level, using a variety of cutting-edge techniques (live cell imaging, single-cell RNA-seq, single-molecule RNA FISH, and single-cell ATAC-seq) in chicken DT40 B-cells.

    In the fist half of the paper, the authors perform a detailed characterization of the cell-to-cell variation arising from a homogeneous stimulation with various doses of anti-IgM. They observe that the NFKB TF RelA forms clear nuclear 'foci' upon stimulation in DT40 cells: this was anecdotally shown in a different cell-type by the same authors in ref 7, but (to my knowledge) has never been systematically studied. This allows them to quantitatively analyse the foci formed in response to stimulation, and they show that this is dose-dependent, heterogeneous and biomodal, and exhibits properties of cooperativity. In parallel, the authors analyse the resulting stimulus-driven changes in gene expression, first using single-cell RNA-seq, and then, elegantly, using RNA FISH, which allows them to directly compare the number of RelA foci to gene expression in individual cells. Like the RelA foci, they find that cell-to-cell gene expression is heterogeneous and bimodal (this has been described before). Interestingly, though, they are able to show that individual stimulus-responsive genes exhibit distinct patterns of cell-to-cell hetereogeneity: they can categorize 4 clusters of responding genes according to different patterns of cell-to-cell variation at distinct stimulus doses, and moreover they show that while the heterogeneity of NFKBIA arises due to bimodal expression levels, that of CD83 is simply due to broad variation between cells. Although focused on NFkB, there is a lot of information here with some important (and non-intuitive) implications that could apply to many other stimulus-driven or developmental responses that exhibit heterogeneous patterns of gene expression. A more in-depth analysis of the single-cell datasets would certainly be very worthwhile and fruitful.

    In the second half of the paper, the authors attempt to use their single-cell data, alongside ATAC-seq genomic analyses, to draw inferences about how or whether the model genes NFKBIA and CD83 are regulated by super-enhancers (SEs). Both of these genes are associated with SEs that gain accessibility upon stimulation (recapitulating the authors' findings in ref 8 in a different cell-type), and the CD83 promoter exhibits co-accessibility with two regions within an adjacent SE. The authors also show that both genes are sensitive to treatment with 1.6-HD, a compound that disrupts liquid-like condensates (a characteristic that has been reported for SEs), and CD83 is sensitive to an inhibitor of Brd4 (which has been associated to SE function). However, while these findings could be considered to be suggestive of regulation by SEs, they are clearly not definitive (nor do the authors claim so).

    Finally, the authors show (figure 4a-c) that while the level of stimulus-driven gene upregulation correlates with co-accessibility with both SEs and typical enhancers (TEs), the cell-to-cell heterogeneity of gene expression correlates only with co-accessibility with SEs. This would agree with a model in which SE-regulated gene regulation may generally impart heterogeneous or switch-like gene expression. *

    **Specific comments**

    *• The experiments are adequately presented, and the authors indicate that not only the sequencing data but also the analysis code is available. Nevertheless, the methods section is rather terse, and could benefit from more detail to understand the various analyses, particularly concerning the analyses of SEs in figures 3 and S7, where it is often difficult to understand how peaks or genes are categorized. *

    Response: We thank the Reviewer for pointing this out and we agree that the Methods section was not described in detail, particularly in how the SEs were analyzed and categorized. Therefore, we have added more details on how SEs were categorized in the Methods section as follows:

    “ 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 “makeTagDirectory” 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 fold-change was then calculated between both conditions for the merged peaks separately for SE and TE. Genes associated with both SE and TE were assigned only to the SE.”

    Similarly, we have added more details for other analyses in the Method section and the main sentences.

    • The imaging, scRNA-seq and RNA-FISH experiments are well-presented, although the supplementary figures 4 and 5 include key results that would merit inclusion within the main figures. *

    Response: We thank the Reviewer for this comment. We have included supplementary figures 4b and 5d in the main figures (new Fig. 2g) since both of these figures represent the raw data revealing the differences between smFISH counts and RNA-seq derived gene expression.

    • It is strking that although all the conclusions about SEs are drawn almost exclusively from analysis of ATAC-seq data, no raw ATAC-seq data is directly shown in any figure (even in the browser snapshots of figure 4d & e). It would be important to show the actual ATAC data from which the inferences of figures 3 and 4 are drawn, especially so that it is possible to visualize the implication of a particular 'ATAC fold-change' or of 'ATAC-gained enhancers'. * Response: We have added a browser snapshot of the ATAC-seq data, presenting the super-enhancer region assigned to both CD83 and NFKBIA (new Fig. 3c).

    Reviewer #1 (Significance (Required)):

    *• This manuscript can be considered as a follow-up of the authors' previous paper (Michida 2020, ref 8), here focusing on cell-to-cell heterogeneity rather than on the overall magnitude of the stimulus-induced response. Overall, the experiments are well-performed and bring new data to an interesting angle of gene regulation. However, the analyses presented do not seem to fully exploit the data, and the authors do not manage to present any strong conclusions, particularly relating to the possible involvement of super enhancers. *

    Response: To strengthen our conclusions about the possible involvement of super-enhancers in regulating heterogeneity, we performed additional analyses on the properties of the SE including the number of transcription factors, NF-κB and PU.1 binding motifs and the length of the enhancers, according to a previous report (Michida et al., 2020, Cell Rep). This was also conducted to confirm whether the ATAC-seq-based SE identification method presents results consistent with those provided by H3K27Ac-ChIP-based methods utilized in the previous study (Michida et al., 2020, Cell Rep). SEs revealed longer genomic length (new Supplementary Fig. 8a) and this length was positively correlated with the ATAC signal (new Supplementary Fig. 8b). Furthermore, gained and lost SE revealed a correlation with enhanced gene expression upregulation and downregulation, respectively, compared to TE (new Fig. 3g). We also demonstrated that SE-regulated genes have a higher Fano factor change, which is consistent with the state of an SE whether it is gained or lost (new Fig. 5a, 5b). For binding motif analysis, we observed a slightly higher PU.1 motif density at SEs (new Supplementary Fig. 11), corresponding to the results of the previous study (Michida et al., 2020, Cell Rep). Interestingly, only the density of NF-κB and not PU.1 was correlated with ATAC signal change in SE (new Fig. 4a), suggesting that those SEs were controlled by nuclear translocation of NF-κB.

    As a mechanism to produce gene expression heterogeneity in phenotypically identical cells, we observed that co-accessibility, which has been reported to be concordant with genomic contacts is correlated to Fano factor change, indicating that gene expression heterogeneity possibly stems from cis-regulatory interactions. NF-κB activation has been reported to increase the heterogeneity in some genes and is attributed to the accumulation of Ser5p RNAPII (Wong et al., 2018, Cell Rep). Additionally, Ser5p RNAPII has been reported to accumulate at enhancer regions (Koch et al., 2011, Nat Struct Mol Biol), and that the accumulation of RNAPII is suggested to assist in gene expression activation through enhancer-promoter contact (Thomas et al., 2021, Mol Cell). Our results support these conclusions since co-accessibility or putative cis-regulatory interactions correlate to Fano factor changes. SE can form phase-separated transcription hubs containing multiple enhancers and/or promoters, which may enable the higher diffusion rate of active enhancers; therefore, it may induce a higher possibility of genomic DNA interactions (Gu et al., 2018, Science). In contrast, the enrichment of TATA motif has also been proposed to generate transcriptional heterogeneity (Faure et al., 2017, Cell Syst). Therefore, we examined this possibility with our data. However, we observed a higher occurrence of TATA box in genes associated with lost SE (new Supplementary Fig. 18) which might have caused gene expression heterogeneity in unstimulated cells. This heterogeneity might be due to the differences in Pol II loading intervals (Tunnacliffe & Chubb, 2020, Trends Genet) however the noise associated with gained SE is possibly generated by the fluctuation of high-order biomolecular assembly. Therefore, we believe that the source of heterogeneity in these conditions were different.

    Additionally, we performed Hill function analysis to reveal the threshold behavior of gene expression in our analysis since previously gained SEs were associated with threshold gene expression (Michida et al., 2020, Cell Rep). In this study, we presented that threshold behavior in gained SE is related to motif density of NF-κB (Fig. 4d), however, threshold behavior does not seem to be related to heterogeneous gene expression.

    Following these results, we concluded that NF-κB activated SE has two closely related but distinct functions for gene control: (1) enhanced heterogeneity and fold-changes and (2) switch-like expression. These are controlled by different mechanisms stemming from chromatin status: (1) frequency of cis-regulatory genomic interactions possibly mediated by phase separation and (2) cooperative binding of NF-κB to DNA. These differences were well represented by expression profiles of CD83 (higher heterogeneity and weak bimodal expression) and *NFKBIA *(lower heterogeneity and strong bimodal expression).

    • For instance, the existence of multiple gene clusters that exhibit distinct patterns of heterogeneity implies that switch-like gene activation occurs on a per-gene basis, rather than corresponding to an all-or-nothing activation of individual cells. This would be an exciting finding, and the authors have the data to test this. Likewise, the division of heterogeneous gene expression into bimodal (like NFKBIA) or unimodal (like CD83) distributions could be a nice paradigm if systematically applied to the other 1335 differentially-expressed genes identified by the authors. * Response: We appreciate this comment. Following your comment, we analyzed the relationship between heterogeneity and bimodality (switch-like expression or high Hill coefficient) for the remaining genes. We observed that SE having a high Hill coefficient contained a higher number of NF-κB motif in SE (new Fig. 4), indicating that cooperative binding of NF-κB to DNA shaped non-linear gene expression profiles as we indicated in a previous paper (Michida et al., 2020, Cell Rep). Additionally, as described in the earlier section, we observed that heterogeneity arises from cis-regulatory genomic interaction. We compared these gene groups and observed that these properties were not completely shared (new Supplementary Fig. 15), indicating that bimodality and heterogeneity originated from different mechanisms. We assume that those differences are mediated through a combination of chromatin accessibility and the biophysical properties of NF-κB.

    *In contrast, although the authors try to use their data to investigate gene regulation by SEs, these inferences are all somewhat indirect, and the authors themselves do not manage to draw any definitive conclusions. * Response: We appreciate this comment. We performed the additional computational analysis and carefully interpreted the data. Additionally, we have now concluded that SEs have two major biological functions: (1) gene expression heterogeneity, which is mediated via cis-regulatory interactions (Fig. 5) and (2) bimodal gene expression, which is mediated by NF-κB binding (new Fig. 4). The latter finding has also been reported in a mouse primary B cell, albeit the mechanism causing heterogeneity was a novel conclusion of this study.

    I feel that the authors are under-selling their data here. As-is, the data represents more of a resource than a study with a clear message, but I believe that with more in-depth analysis the authors could make a much more significant advance, particularly concerning the cell-to-cell heterogeneity of gene expression. I would be very enthusiastic to review the same data again with a more detailed analysis, which I believe would enormously improve the manuscript. Response: We appreciate this comment. As described in this report and the revised manuscript, we performed a considerably detailed computational analysis and gained several novel insights to answer the question regarding the functional roles of SE. We are grateful to learn that gene expression patterns may be estimated from ATAC-seq profiles and that they may even be controlled. We hope that this Reviewer would observe the scientific value of our study and provide us with your valuable feedback on our revised manuscript.

    Reviewer #2

    Reviewer #2 (Evidence, reproducibility and clarity (Required)):

    *Imaging and single cell sequencing analyses of super-enhancer activation mediated by NF-κB in B cells" by Wibisana et al. examined the relationship between super-enhancers, NF-κB nuclear aggregation, and target gene regulation. The authors have generated a large amount of data from fluorescent microscopy, scRNA-seq, scATAC-seq, smRNA FISH. While this is an impressive dataset in terms of diverse technically advanced methods employed, it is not clear what to take as a main conceptual advance. What could be the functional implications of observed cell-cell variability in B cell transcriptional responses to environmental stimuli? In addition to this general point, the following are specific comments that could improve the manuscript. *

    • In Figure 2, smRNA FISH foci of CD83 and NFKBIA are quantified as # of spots per cell (Supplementary figure 5). But it is difficult to see in Figure 2 the colocalization of any mRNA spots with RelA foci. Ideally, it will be convincing to show by DNA FISH that these target loci are indeed located within NF-κB occupied super-enhancer puncta. Even with the current RNA FISH data, some colocalization analysis could have been performed. * Response: In Figure 2, we were unable to perform accurate colocalization analysis with the current smFISH data as the probes used by us map to exons. Moreover, we have also previously performed DNA-FISH; nevertheless, it was difficult to assess co-localization between the DNA and RelA proteins secondary to the degradation of RelA-GFP proteins. Therefore, we decided to perform intronic smRNA-FISH, which can be used to pinpoint the site of active transcription (Levesque and Raj, 2013, Nat Methods). The results, along with the quantification results, are presented in the new Fig. 2f.
    • Supplementary Figure 5a shows lower correlations of # GFP-RelA foci to CD83 transcripts in comparison to NFKBIA. Even though the foci and smRNA FISH spots are derived from high resolution imaging data, we should remember that any snapshot measurements have limited information content for gene regulatory relationships. Live cell studies (for example, from the groups of Suzanne Gaudet, Kathryn Miller-Jensen, and Myong-Hee Sung) have shown that time-integrated measures (e.g. maximum fold change and area under the curve of RelA signaling time course in single cells) are better correlates to transcriptional output of target genes (Lee REC et al 2014 Mol Cell; Wong VC et al 2019 Biophysical J; Sung MH et al. 2014 Science Signal; Martin EW et al. 2020 Science Signal). *

    Response: We thank the Reviewer for this valuable comment. One of the reasons for a lower correlation between GFP-RelA foci and CD83 transcripts compared to NFKBIA may be the difference in expression timing of CD83 and NFKBIA and the timing of nuclear localization of GFP-RelA. RelA localizes in the nucleus 10−30 mins after cell stimulation, and NFKBIA is an early responsive gene, however, CD83 is expressed later (new Supplementary Fig. 17). Therefore, this time difference possibly affects correlation accuracy. Although we agree that high-throughput time-course measurement of RelA-GFP combined with smFISH measurements, such as that reported in Wong VC et al., 2019, will be ideal, it is technically difficult since DT40 are suspension cells and the smFISH protocol requires multiple washing and centrifugation steps. Thus, with this experimental setup, we were unable to perform the time-course analysis.

    Nonetheless, we measured the time-course foci formation at the same single-cells (new Supplementary Fig. 1b) and observed that it effectively represents Figure 1a, which is a snapshot of the population dynamics of RelA foci across time. Additionally, the observed dynamics, which revealed a steep initial increase and slight decrease with time, effectively recapitulates the previous reports (Lee et al., 2014, Mol Cell; Wong et al., 2019, Biophys J).

    In our analysis, we performed imaging analysis to demonstrate that NF-κB foci formation is switch-like, and this formation might be involved in the formation of phase-separated condensates enhancing DNA to DNA contact. The number of foci may depend upon the intracellular concentration of NF-κB, and fold change in the RelA signal may be correlated with gene expression as previously reported (Lee et al., 2014, Mol Cell; Wong et al., 2019, Biophys J). However, there is another report presenting that promote/enhancer proximity is not related to gene expression (Alexander et al. 2019, eLife). Although we were unable to perform this analysis owing to the limitations stated above, we tried to find the relationship between RelA foci and gene expression by performing biochemical perturbations (Fig 1e-f, Fig 5h) and presented that these foci are related to gene expression.

    • The analyses have been performed using DT40 cells. In the Methods section, no description was provided about what type of B cells DT40 is, even though few outside of the field may not know that the cells were immortalized from chicken. This is an important consideration, because some nuclear bodies and genome organization features are different between host species and they also depend on whether the cells are primary or transformed. Because the authors do not discuss this point, it seems possible that the findings about NF-κB aggregates and super-enhancers may not necessarily hold true for primary B cells. *

    Response: We thank the Reviewer for pointing out these issues. We have added the following description on DT40 cells in the Methods section describing that DT40 cells are chicken bursal lymphoma cells.

    DT40 B lymphocytes have been widely used as a B cell model for studying B cell receptor signaling (Mori et al., 2002, J. Exp. Med.; Patterson et al., 2002, Cell; Saeki et al., 2003, EMBO J.) due to its high gene targeting efficiency. We also previously confirmed that anti-IgM stimulation induces the NF-κB signaling pathway in mouse primary splenic B cells and DT40 and that the signaling molecules and dynamics in these cells are well conserved (Shinohara et al., 2014, Science; Shinohara et al., 2016, Sci. Rep.; Inoue et al., 2016, NPJ Syst. Biol. Appl.). However, we understand the Reviewer’s concerns. Therefore, we have provided the track view of primary B cell ATAC-seq data to demonstrate that the chromatin accessibility changes upon anti-IgM stimulation in CD83 and NFKBIA were similarly observed in primary B cell data (new Supplementary Fig. 9b) and that the upregulation and association with SE of CD83 and NFKBIA were also observed in primary B cell (new Supplementary Fig. 9a).

    • Similarly, the GFP-RelA expressing DT40 cell generation should be described with more detail (beyond "provided by ..."). N-terminal or C-terminal fusion? Did the fusion construct contain an artificial promoter (e.g. CMV) or an upstream fragment of the genomic Rela locus (chicken or human)? Methods of transfection and cloning of stable lines? These choices affect the interpretation of the data, so they must be fully described and justified. *

    Response: We thank you for pointing this out. We have added the following details on the RelA-GFP construct in the Methods section:

    Mouse RelA-eGFP with eGFP on the C terminal was cloned into a pGAP vector containing Ecogpt resistance gene targeting endogenous GAPDH locus. This construct was further electroporated into wild-type cells and selected using Ecogpt to produce RelA-GFP-expressing DT40 cells.

    • DT40 cells were cultured in 39 degrees. Michael White and colleagues have shown that high temperatures can alter NF-kappaB dynamics and function (https://www.pnas.org/content/115/22/E5243). Did the authors try lower temperatures to ascertain that the NF-kB aggregates and other major findings are still observed in 37 degrees? *

    Response: We performed the experiments at 39 degrees to mimic the natural body temperature of chicken since DT40 cells were derived from chicken bursal lymphoma (Saribasak and Arikawa, 2006, Subcell Biochem). Previously, we cultured DT40 cells at 37 degrees and observed that the cell growth was inhibited, and thus, we believed that it was not ideal to perform experiments of DT40 cells at 37 degrees.

    Reviewer #2 (Significance (Required)):

    *It is not clear what to take as a main conceptual advance. *

    Response: Considering the original manuscript, we agree with the Reviewer on the lack of strong emphasis on the conclusions of our study. Therefore, in this revised manuscript, we have focused on the comprehensive mechanism of heterogeneity and switch-like activation in gene expression control. As we described in the comments to Reviewer #1, we performed an additional in-depth computational analysis on SE and TE. Consequently, we demonstrated that enhanced heterogeneity and expression fold-changes mediated by SE are defined by the number of cis-regulatory genomic interactions in open chromatin regions (Figure 5), however, switch-like expression (bimodal patterns) is determined by the number of NF-κB binding in SE (new Figure 4). The latter finding has also been reported in a mouse primary B cell in our previous study (Michida et al. 2020, Cell Rep.). However, the mechanism causing heterogeneity is a novel conclusion obtained in this study. We also concluded that these similar, albeit quantitatively and slightly different characteristics in gene control can be achieved through a combination of chromatin accessibility of host cells and biophysical properties of NF-κB molecule, which is involved in phase separation.

    What could be the functional implications of observed cell-cell variability in B cell transcriptional responses to environmental stimuli? *

    Response: We performed gene ontology analysis to reveal how the heterogeneously expressed genes (cluster 4) (Fig. 2d) presented enrichment for immune-related functions (Supplementary Fig. 5b). This result supports a previous study, which stated that variability in gene expression is related to function (Osorio et al., 2019, Cells).

    This discussion is incorporated in the manuscript as follows:

    “We observed that genes with an increased heterogeneity upon increasing stimulation dose are enriched with cell-type-specific immune regulatory genes (Supplementary Fig. 5b), supporting a previous report where heterogeneity in gene expression is tied to biological functions and may be used by cells as a bet-hedging or a response distribution mechanism (Osorio et al., 2019, Cells), where cells exhibit heterogeneity to enable 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.”

  3. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #2

    Evidence, reproducibility and clarity

    "Imaging and single cell sequencing analyses of super-enhancer activation mediated by NF-κB in B cells" by Wibisana et al. examined the relationship between super-enhancers, NF-κB nuclear aggregation, and target gene regulation. The authors have generated a large amount of data from fluorescent microscopy, scRNA-seq, scATAC-seq, smRNA FISH. While this is an impressive dataset in terms of diverse technically advanced methods employed, it is not clear what to take as a main conceptual advance. What could be the functional implications of observed cell-cell variability in B cell transcriptional responses to environmental stimuli? In addition to this general point, the following are specific comments that could improve the manuscript.

    1. In Figure 2, smRNA FISH foci of CD83 and NFKBIA are quantified as # of spots per cell (Supplementary figure 5). But it is difficult to see in Figure 2 the colocalization of any mRNA spots with RelA foci. Ideally, it will be convincing to show by DNA FISH that these target loci are indeed located within NF-κB occupied super-enhancer puncta. Even with the current RNA FISH data, some colocalization analysis could have been performed.
    2. Supplementary Figure 5a shows lower correlations of # GFP-RelA foci to CD83 transcripts in comparison to NFKBIA. Even though the foci and smRNA FISH spots are derived from high resolution imaging data, we should remember that any snapshot measurements have limited information content for gene regulatory relationships. Live cell studies (for example, from the groups of Suzanne Gaudet, Kathryn Miller-Jensen, and Myong-Hee Sung) have shown that time-integrated measures (e.g. maximum fold change and area under the curve of RelA signaling time course in single cells) are better correlates to transcriptional output of target genes (Lee REC et al 2014 Mol Cell; Wong VC et al 2019 Biophysical J; Sung MH et al. 2014 Science Signal; Martin EW et al. 2020 Science Signal).
    3. The analyses have been performed using DT40 cells. In the Methods section, no description was provided about what type of B cells DT40 is, even though few outside of the field may not know that the cells were immortalized from chicken. This is an important consideration, because some nuclear bodies and genome organization features are different between host species and they also depend on whether the cells are primary or transformed. Because the authors do not discuss this point, it seems possible that the findings about NF-κB aggregates and super-enhancers may not necessarily hold true for primary B cells.
    4. Similarly, the GFP-RelA expressing DT40 cell generation should be described with more detail (beyond "provided by ..."). N-terminal or C-terminal fusion? Did the fusion construct contain an artificial promoter (e.g. CMV) or an upstream fragment of the genomic Rela locus (chicken or human)? Methods of transfection and cloning of stable lines? These choices affect the interpretation of the data, so they must be fully described and justified.
    5. DT40 cells were cultured in 39 degrees. Michael White and colleagues have shown that high temperatures can alter NF-kappaB dynamics and function (https://www.pnas.org/content/115/22/E5243). Did the authors try lower temperatures to ascertain that the NF-kB aggregates and other major findings are still observed in 37 degrees?

    Significance

    It is not clear what to take as a main conceptual advance.

    What could be the functional implications of observed cell-cell variability in B cell transcriptional responses to environmental stimuli?

    Referee cross-commenting

    I concur with Reviewer #1's comments about systematic grouping of 1335 differentially expressed genes based on heterogeneity, and also about showing raw ATAC-seq data tracks and plots. We both commented that the study lacks a significant conclusion in its current form.

  4. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #1

    Evidence, reproducibility and clarity

    Summary of findings & key conclusions

    The manuscript by Wibisana et al. describes an impressive set of experiments that analyse the NFkB response at the single-cell level, using a variety of cutting-edge techniques (live cell imaging, single-cell RNA-seq, single-molecule RNA FISH, and single-cell ATAC-seq) in chicken DT40 B-cells.

    In the fist half of the paper, the authors perform a detailed characterization of the cell-to-cell variation arising from a homogeneous stimulation with various doses of anti-IgM. They observe that the NFKB TF RelA forms clear nuclear 'foci' upon stimulation in DT40 cells: this was anecdotally shown in a different cell-type by the same authors in ref 7, but (to my knowledge) has never been systematically studied. This allows them to quantitatively analyse the foci formed in response to stimulation, and they show that this is dose-dependent, heterogeneous and biomodal, and exhibits properties of cooperativity. In parallel, the authors analyse the resulting stimulus-driven changes in gene expression, first using single-cell RNA-seq, and then, elegantly, using RNA FISH, which allows them to directly compare the number of RelA foci to gene expression in individual cells. Like the RelA foci, they find that cell-to-cell gene expression is heterogeneous and bimodal (this has been described before). Interestingly, though, they are able to show that individual stimulus-responsive genes exhibit distinct patterns of cell-to-cell hetereogeneity: they can categorize 4 clusters of responding genes according to different patterns of cell-to-cell variation at distinct stimulus doses, and moreover they show that while the heterogeneity of NFKBIA arises due to bimodal expression levels, that of CD83 is simply due to broad variation between cells. Although focused on NFkB, there is a lot of information here with some important (and non-intuitive) implications that could apply to many other stimulus-driven or developmental responses that exhibit heterogeneous patterns of gene expression. A more in-depth analysis of the single-cell datasets would certainly be very worthwhile and fruitful.

    In the second half of the paper, the authors attempt to use their single-cell data, alongside ATAC-seq genomic analyses, to draw inferences about how or whether the model genes NFKBIA and CD83 are regulated by super-enhancers (SEs). Both of these genes are associated with SEs that gain accessibility upon stimulation (recapitulating the authors' findings in ref 8 in a different cell-type), and the CD83 promoter exhibits co-accessibility with two regions within an adjacent SE. The authors also show that both genes are sensitive to treatment with 1.6-HD, a compound that disrupts liquid-like condensates (a characteristic that has been reported for SEs), and CD83 is sensitive to an inhibitor of Brd4 (which has been associated to SE function). However, while these findings could be considered to be suggestive of regulation by SEs, they are clearly not definitive (nor do the authors claim so).

    Finally, the authors show (figure 4a-c) that while the level of stimulus-driven gene upregulation correlates with co-accessibility with both SEs and typical enhancers (TEs), the cell-to-cell heterogeneity of gene expression correlates only with co-accessibility with SEs. This would agree with a model in which SE-regulated gene regulation may generally impart heterogeneous or switch-like gene expression.

    Specific comments

    • The experiments are adequately presented, and the authors indicate that not only the sequencing data but also the analysis code is available. Nevertheless, the methods section is rather terse, and could benefit from more detail to understand the various analyses, particularly concerning the analyses of SEs in figures 3 and S7, where it is often difficult to understand how peaks or genes are categorized.

    • The imaging, scRNA-seq and RNA-FISH experiments are well-presented, although the supplementary figures 4 and 5 include key results that would merit inclusion within the main figures.

    • It is strking that although all the conclusions about SEs are drawn almost exclusively from analysis of ATAC-seq data, no raw ATAC-seq data is directly shown in any figure (even in the browser snapshots of figure 4d & e). It would be important to show the actual ATAC data from which the inferences of figures 3 and 4 are drawn, especially so that it is possible to visualize the implication of a particular 'ATAC fold-change' or of 'ATAC-gained enhancers'.

    Significance

    Significance

    This manuscript can be considered as a follow-up of the authors' previous paper (Michida 2020, ref 8), here focusing on cell-to-cell heterogeneity rather than on the overall magnitude of the stimulus-induced response. Overall, the experiments are well-performed and bring new data to an interesting angle of gene regulation. However, the analyses presented do not seem to fully exploit the data, and the authors do not manage to present any strong conclusions, particularly relating to the possible involvement of super enhancers.

    For instance, the existence of multiple gene clusters that exhibit distinct patterns of heterogeneity implies that switch-like gene activation occurs on a per-gene basis, rather than corresponding to an all-or-nothing activation of individual cells. This would be an exciting finding, and the authors have the data to test this. Likewise, the division of heterogeneous gene expression into bimodal (like NFKBIA) or unimodal (like CD83) distributions could be a nice paradigm if systematically applied to the other 1335 differentially-expressed genes identified by the authors.

    In contrast, although the authors try to use their data to investigate gene regulation by SEs, these inferences are all somewhat indirect, and the authors themselves do not manage to draw any definitive conclusions.

    I feel that the authors are under-selling their data here. As-is, the data represents more of a resource than a study with a clear message, but I believe that with more in-depth analysis the authors could make a much more significant advance, particularly concerning the cell-to-cell heterogeneity of gene expression. I would be very enthusiastic to review the same data again with a more detailed analysis, which I believe would enormously improve the manuscript.

    Reviewer field of expertise

    My expertise is in gene regulation and genomics. I am competent to review the implications of all parts of this paper, and all the technical aspects with the exception of the microscopy.

    Referee Cross-commenting

    I agree both with the specific points raised by reviewer #2, and also with the overall comment that - despite the large amount of data - the authors do not present any clear conceptual advance or tackle the functional implications of their results.

  5. Version published to 10.1101/2021.07.13.452147v2 on bioRxiv
  6. Version published to 10.1101/2021.07.13.452147v1 on bioRxiv