Flipping between Polycomb repressed and active transcriptional states introduces noise in gene expression

Polycomb repressive complexes (PRCs) are important histone modifiers, which silence gene expression; yet, there exists a subset of PRC-bound genes actively transcribed by RNA polymerase II (RNAPII). It is likely that the role of Polycomb repressive complex is to dampen expression of these PRC-active genes. However, it is unclear how this flipping between chromatin states alters the kinetics of transcription. Here, we integrate histone modifications and RNAPII states derived from bulk ChIP-seq data with single-cell RNA-sequencing data. We find that Polycomb repressive complex-active genes have greater cell-to-cell variation in expression than active genes, and these results are validated by knockout experiments. We also show that PRC-active genes are clustered on chromosomes in both two and three dimensions, and interactions with active enhancers promote a stabilization of gene expression noise. These findings provide new insights into how chromatin regulation modulates stochastic gene expression and transcriptional bursting, with implications for regulation of pluripotency and development.


Response to Reviewer #3:
We thank Rev. #3 for his/her comments, which enabled us to clarify important points in our manuscript. The remarks of Rev. #3 are denoted in italics.
The revised manuscript does not address my core concerns, and contains a number of incorrect and unsupported statements. I repeat my major concerns below: 1. The response claim that this paper is the ' first transcriptome-wide study showing that differences in chromatin state can affect gene expression kinetics', and that none of the previous papers addressed the mechanistic nature of expression variation. These statements are incorrect, as has also been pointed out by the other reviewers. Both Kumar and Kolodziejczyk --in the same system --link expression heterogeneity to changes in chromatin state (in the latter case, by showing that differences exist across serum conditions, which have widely been linked to changes in H3K27me3 levels). The former manuscript ties changes in heterogeneity to miRNAproduction, and explicitly tests the relationship between H3K27me3 and expression heterogeneity.
The reviewer's statement about our earlier work (Kolodziejczyk et al. Cell Stem Cell 2015) is incorrect; this study has addressed the roles of pluripotency and cell cycle genes in mESCs across different conditions, however, histone modifications or H3K27me3 levels are not studied explicitly in any part of this study.
It is true that there may be implicit information in an indirect sense via cell culture conditions in these previous papers, but such global relationships do not provide a direct molecular insight as shown in our work here. To be specific, while Kumar et al. addressed the relationship between H3K27me3 levels and expression heterogeneity, this paper or previous papers do not focus on the mechanistic role of Polycomb in regulating transcriptional bursting. Furthermore, Kumar et al.
provided an analysis of serum/LIF-cultured cells, and we have shown that these cells are a mixture of multiple distinct biological states (Kolodziejczyk et al., Cell Stem Cell, 2015). This makes these samples unsuitable to study transcriptional bursting, as much of the variation in gene expression will be due to cell subpopulations, rather than bona fide transcriptional bursting that can be ascribed to chromatin state.
The LIF conditions under Oct4 selection that we focus on in this manuscript are in this sense ideal. With the present study, we do not wish to address the differences between naïve and less naïve states of differentiation, and instead our work here is focused on Polycomb repression and its impact on transcriptional cell-to-cell variation. Our paper is conceptually entirely distinct from the previous papers cited in this sense.
We had already described this in the previous revised manuscript. Now, in the revised version, we cite Kumar et al. in the discussion part of the manuscript as well, and have modified the text as follows: "Earlier work indicated that expression of Polycomb target genes negatively correlates with levels of H3K27me3, and suggested that dynamic fluctuations in chromatin state are associated with expression of certain Polycomb targets in pluripotent stem cells 38 ." Additionally, we have now modified the last paragraph of the introduction as follows: "Motivated by this hypothesis, we integrate states of histone and RNAPII modification from a published classification of ChIP-Seq data 5 with single-cell RNA-sequencing data generated for this analysis. The matched chromatin and scRNA-seq data sets allow us to decipher, on a genome-wide scale, how differences in chromatin state affect transcriptional kinetics." 2. I continue to find the kinetic analyses to be fundamentally flawed. The core of my argument is that the kinetic model used does not incorporate or account for technical noise in any way. The revision is clear that there are 20-40% technical noise, therefore the parameters inferred from this model are not meaningful.
The reviewer is right that our kinetic analysis does not account for technical noise, as we have already explained this in the previous rebuttal letter: "Since our data do not contain external spike-in molecules (the only way to incorporate technical noise in our kinetic model), we addressed this point by focusing on moderately or highly expressed genes with an expression cutoff of 10. The assumption is that technical noise for these genes is small enough to estimate kinetic parameters accurately." Additionally, in the previous rebuttal letter, we have also shown that the fraction of technical noise saturates to 20% at the mean normalised count of 100. Now, to make this point clearer, in the revised manuscript, we have also explained this in the Methods section of the manuscript as follows: "We should note that our kinetic analyses do not account for technical noise as our data do not contain external spike-in molecules (the only way to incorporate technical noise in our kinetic model).
Therefore, we addressed this point by focusing on moderately or highly expressed genes with an expression cutoff of 10. The assumption is that technical noise for these genes is small enough to estimate kinetic parameters accurately." While the response states that the authors 'do not change', by changing the cutoff from 10 reads to 100 reads, it is clear that there is a significant shift. Looking at the change in burst frequency in the original manuscript, there is an enormous shift with a highly signifcant p-value. Moving this cutoff (corresponding to a reduction in technical noise of 40-20%), the distributions appear nearly overlapping, and the p-value is now >0.01. To suggest that these results are unchanged is misleading and inaccurate -a key result ( Figure 1A) of the manuscript is heavily dependent on the selection of an inproper cutoff, and is flawed.
In the previous rebuttal letter, we have shown distributions of noise and burst frequency levels with a stricter cutoff of 100 mean normalized counts. However, we have only shown this considering a subset of genes (i.e. expression matched group where there are 262 genes per group). The reviewer has compared these distributions to Figure 2A and erroneously inferred that there is a significant shift.
In Figure 2A, we show the whole set of Active and PRCa genes (5,428 genes in total).
We apologize for not clearly highlighting this in the previous rebuttal letter. Now, in the revised version, we show distributions of noise and burst frequency levels across different expression cutoffs: from a cutoff of 10 mean normalized counts to 100. This is now directly comparable to Figure 2A.
This analysis implies that there is highly significant difference between Active and PRCa genes in terms of noise and kinetic parameters, and that our results are robust to changes in selection of  540, 479, 427, 374, 341, 304, 280 and 262 for expression cutoffs of 10, 20, 30, 40, 50, 60, 70, 80, 90 and 100, respectively.
To reiterate, the choice of 10 reads/cell is an inappropriate, cutoff, since when taking the total reads per cell into account, this represents ~1 molecule per cell. This can hardly be considered to be a moderately to highly expressed gene. Even multiplying this by 10 does not reach a threshold that can be considered to be immune from technical noise in scRNA-seq experiments.
Using our 2i mESC scRNA-seq data with ERCC spike-ins (Kolodziejczyk et al. Cell Stem Cell. 2015), we find that the cutoff of 100 normalized counts per cell corresponds to 16.57 transcripts per cell (Please see Rebuttal Figure 1 below). In our previous rebuttal, we have shown that the fraction of technical noise at the cutoff of 100 normalized counts saturates to 20% of technical noise on average. The 20% of technical noise can be considered as the minimum achievable technical noise level, and therefore the cutoff of 100 is a reasonable threshold minimizing the effect of technical noise. At this threshold, we have shown that our key result regarding Figure 2A  transcripts per cell according to our 2i mESC scRNA-seq data with ERCC spike-ins (Kolodziejczyk et al. Cell Stem Cell. 2015).
3. The response states that RNA FISH and single cell RNA-seq results consistently agree, and therefore there is no need for external validation. This is again an inaccurate statement. Numerous studies have described the extensive technical biases introduced by amplification of single cell RNA-seq data, which could explain the results the authors see here, but are not present in RNA FISH.
We disagree with this statement. In contrast, many studies have indicated that RNA FISH and We agree that there would be amplification biases in single cell RNA-seq data for lowly expressed genes, and we analyse sensitivity and specificity in depth in our power analysis of scRNA-seq protocols available online and accepted in principle in Nature Methods: http://biorxiv.org/content/early/2016/09/08/073692 In the work in the present manuscript, we are using a highly sensitive and specific SMART-seq protocol in Fluidigm C1 microfluidic platform, with a saturating sequencing depth of over a million reads per cell. This is one of the most sensitive and accurate scRNA-seq pipelines available, including microfluidic cell capture and full length transcript coverage. Moreover, we are comparing set of genes within this data set with comparable technical noise contribution, between PRC and non-PRC marked gene sets. In addition, our noise measure "DM" corrects for gene expression levels.