Re-evaluating the role of nucleosomal bivalency in early development

Nucleosomes, composed of DNA and histone proteins, represent the fundamental repeating unit of the eukaryotic genome1; posttranslational modifications of these histone proteins influence the activity of the associated genomic regions to regulate cell identity2–4. Traditionally, trimethylation of histone 3-lysine 4 (H3K4me3) is associated with transcriptional initiation5–10, whereas trimethylation of H3K27 (H3K27me3) is considered transcriptionally repressive11–15. The apparent juxtaposition of these opposing marks, termed “bivalent domains”16–18, was proposed to specifically demarcate of small set transcriptionally-poised lineage-commitment genes that resolve to one constituent modification through differentiation, thereby determining transcriptional status19–22. Since then, many thousands of studies have canonized the bivalency model as a chromatin hallmark of development in many cell types. However, these conclusions are largely based on chromatin immunoprecipitations (ChIP) with significant methodological problems hampering their interpretation. Absent direct quantitative measurements, it has been difficult to evaluate the strength of the bivalency model. Here, we present reICeChIP, a calibrated sequential ChIP method to quantitatively measure H3K4me3/H3K27me3 bivalency genome-wide, addressing the limitations of prior measurements. With reICeChIP, we profile bivalency through the differentiation paradigm that first established this model16,18: from naïve mouse embryonic stem cells (mESCs) into neuronal progenitor cells (NPCs). Our results cast doubt on every aspect of the bivalency model; in this context, we find that bivalency is widespread, does not resolve with differentiation, and is neither sensitive nor specific for identifying poised developmental genes or gene expression status more broadly. Our findings caution against interpreting bivalent domains as specific markers of developmentally poised genes.

In its original conception, the bivalency model posits that the combination of H3K4me3 24 and H3K27me3 represents a specific regulatory marker of developmentally staged genes. 25 Specifically, lineage commitment genes are thought to be held in a low-expression, 26 transcriptionally "poised" state by promoter nucleosomes bearing both H3K4me3 and 27 H3K27me3 16,18,21,22 . Upon differentiation, the bivalent domain "resolves" into a monovalent 28 state, and the associated gene is either transcriptionally activated or terminally repressed if 29 H3K27me3 or H3K4me3 is lost, respectively 16,18,21,22 . The elegance of this instructive model 30 inspired a host of follow-on studies that have suggested that bivalency is important in 31 differentiation [23][24][25][26][27][28][29][30] , embryogenesis 17,31-34 , genome architecture 22,35-38 , and oncogenesis [39][40][41][42][43] . 32 In the absence of unambiguous biochemical or functional validation 17,44,45 , these studies 33 have largely relied upon ChIP, with the vast majority of studies defining loci with independent 34 ChIP enrichment for H3K4me3 and H3K27me3 as bivalent domains. However, this analysis 35 cannot distinguish whether the two modifications coexist or represent two distinctly marked 36 subpopulations of alleles or cells. Further, because different ChIPs are normalized separately, 37 they exist on separate relative scales and cannot be quantitatively compared without internal 38 calibration [46][47][48] . As such, it is impossible to quantify the extent of bivalency at a given locus or to 39 measure its changes through differentiation. 40 To address the first problem, several studies have used sequential ChIP 16,29,49-51 , 41 measuring coexistence by using the eluent of an IP against H3K4me3 as the substrate for an IP 42 against H3K27me3 (or vice versa). However, these experiments used antibodies of unknown 43 specificity 16,29,49-51 , were uncalibrated, and were often undersampled 50,52 , precluding 44 quantification of the extent of modification. Moreover, many used relatively large chromatin 45 fragments in their pulldowns, making it difficult to determine whether modifications coexisted Metaprofiles of H3K4me3, H3K27me3, and bivalency at promoters identified as bivalent in 292 naïve mESCs (25% HMD threshold), tracked from naïve mESCs to primed mESCs to NPCs. 293 Heatmaps for bivalency are presented in Extended Data Fig. 3f. (g-h) Alluvial plots of 294 dominance and bivalency of genes from (g) naïve mESCs to NPCs or (h) primed mESCs to 295 NPCs. Bivalency [>25% HMD] can be subcategorized into dominance classes by independent 296 ICeChIP for the constituent marks, with H3K27me3 in excess (H3K27me3/H3K4me3 > e 1 ), 297 H3K4me3 dominant (H3K27me3/H3K4me3 < e -1 ), or intermediate ratios (no dominance). (i) 298 Bivalency metaprofiles for gene subsets indicated in panel (h) from -3kb to +3kb relative to the 299 TSS. **** p < 2.2 x 10 -16 . 300   Correspondence and requests for materials should be addressed to A.J.R. 1206 Reprints and permissions information is available at http://www.nature.com/reprints.  As each nucleosome has two H3 protomers, there are several different configurations of 437 bivalency that a bivalent nucleosome can theoretically adopt, each with a different avidity for 438

Extended Data Figure Captions
ChIP pulldown with immobilized antibody. At one extreme, with the highest avidity, is the 439 symmetric cis-bivalency form, where both H3K4 and both H3K27 residues are trimethylated 440 (Extended Data Fig. 1e). This nucleosome has the most epitopes for antibody binding and will 441 thus have the highest avidity in pulldown reflected in apical pulldown efficiency (Extended Data 442 Fig. 1d). At the other extreme, with the lowest avidity, is the trans-bivalency form, where single 443 H3K4me3 and H3K27me3 marks decorate different histone tails (Extended Data Fig. 1e). This 444 has the fewest epitopes for antibody binding and will thus have no avidity in pulldown. 445 This poses a theoretical challenge in normalization and calibration of a ChIP study; 446 because we cannot separately measure trans-bivalency, symmetric cis-bivalency, nor any 447 intermediate states, it is impossible for us to definitively state whether a given locus with a given 448 HMD has relatively few nucleosomes that are symmetric cis-bivalent or whether it has relatively 449 many nucleosomes that are trans-bivalently modified. To accommodate for this limitation, we 450 include two different bivalent calibrants in our set of nucleosome standards: one that is 451 symmetric cis-bivalent and one that is trans-bivalent. The bivalency sequential ChIP can then be 452 normalized to either one of these standards, and because these two cases represent the limits of 453 pulldown avidity, normalization to these calibrants will define the theoretical "range" in which 454 true bivalency HMD (i.e. the proportion of nucleosomes with some bivalent configuration) exists 455 (Extended Data Fig. 1e). We note that, because the signal from calibration to these standards are 456 scalar multiples of each other, we cannot uniquely distinguish these two configurations in the genome. Absent any prior information about the dominant configuration of bivalency, the 458 proportion of bivalently modified nucleosomes at a given locus will exist in the range defined by 459 calibration to symmetric cis-or trans-bivalent standards (Extended Data Fig. 2a). 460 In practice, there are a few reasons why this is not a major concern. First, there is no mass 461 spectrometry evidence that H3K4me3 and H3K27me3 exist on the same histone tail, despite 462 specific enrichment for these marks and sensitive detection limits 49,81 , suggesting that 463 configurations other than trans-bivalency are at most, extremely minor in abundance. Second, 464 the scarcity of these cis-tail modifications is consistent with the bioschemical literature prior to 465 this work that suggests the biogenesis of these cis-tail modifications is enzymatically challenging 466 due to antagonistic allosteric effects (see Supplementary Note 4). Third, even if symmetric cis-467 bivalency does exist at some loci, for the purposes of tracking changes in bivalency across 468 differentiation, we can still observe an increase or decrease in bivalency by this calibration 469 method; we simply cannot precisely discern whether the effect is driven by nucleosomes 470 gaining/losing trans-bivalency, cis-bivalency, or some combination of the two. The overall 471 amount of bivalency would still increase or decrease in all those scenarios, and so long as our 472 choice of calibrant remains consistent, we can still measure that change regardless of the 473 calibrant that we use for our normalization. Therefore, though we have generated datasets using 474 both calibrants, we present analyses of our reICeChIP bivalency pulldowns calibrated to the 475 trans-bivalent standards. 476 477

Supplementary Note 2 478
Throughout this study, we have defined gene promoters to be the region from 0 to +400bp 479 relative to the TSS, representing the +1 and +2 nucleosomes of each gene. These nucleosomes tend to be well-positioned 83 and, accordingly, are most likely to provide us with adequate read 481 depth to robustly quantify each histone modification. This definition is conservative; we find that 482 H3K4me3 and bivalent domains, which tend to be peak-like, have a median breadth of 550bp at 483 bivalent genes (Extended Data Fig. 2d). 484 The width of these domains raises an important point regarding the measurement of 485 histone modification density as a continuous variable. At a given nucleosome in a single allele of 486 a single cell, there are only three possible states for a histone modification: symmetric, 487 asymmetric or not present. However, nucleosome readers do not typically bind only a single 488 nucleosome at a single position; rather, the local density of the modification across multiple 489 nucleosomes is crucial in localizing these effectors through multivalent avidity-based 490 interactions 84-87 . Indeed, we find that the HMD across sequential nucleosomes relative to the 491 TSS is well autocorrelated (Extended Data Fig. 2e). This means that the interpretation of the 492 HMD across a multinucleosomal span becomes more nuanced; a given histone modification may 493 exist at one or more of those nucleosomes. Accordingly, despite the fact that a single nucleosome 494 is essentially ternary in whether it has a given histone modification or not (i.e. HMD of 0% or 495 100%), a region spanning multiple nucleosomes could have an intermediate HMD; it is this latter 496 quantity that is most relevant for the biological function imparted to the nearby genomic regions, 497 and this is the quantity we analyse through this work. The most important of these possibilities is low input depth. The ICeChIP datasets are 505 normalized to the input read depth at every genomic interval to accommodate for differences in 506 local nucleosome density when computing the HMD. However, this means that at regions that 507 are relatively nucleosome-depleted, there will be few reads in the input, meaning that the 508 denominator of the HMD computation is quite small (Methods). This increased Poisson noise in 509 these regions of low input can result in inflated apparent HMD beyond the physical limit of 510 100%. To accommodate for this, we can compute 95% confidence intervals for the HMD of each 511 modification at each genomic position, and these confidence intervals virtually always overlap 512 the physically possible range of HMD values (e.g., Fig. 1c). In naïve mESCs, only 0.5% of the 513 promoters have a bivalency HMD above 100%, and for the vast majority of these promoters 514 (86.1%), the 95% confidence interval error estimate ranges below 100%. The fact the apparent 515 bivalency HMD calibrated by trans-bivalent standards, is broadly constrained to less than 100% 516 further supports the idea that this choice of calibrant is appropriate and not inflationary 517 There are also several other possibilities that are more challenging to accommodate for. 519 First, some regions of the genome are known to be more artefact-prone for sequencing and 520 mapping 88 ; if the IP sample is enriched for these sequences relative to the input, then that could 521 be disproportionately represented in the IP and have an apparent HMD greater than 100%. 522 Second, the antibodies themselves could skew the apparent HMD. If the antibody is capturing 523 substantial off-target material, then that will result in systematic inflation of the IP, resulting in 524 an inflated HMD. Though ICeChIP barcoded nucleosome standards can help monitor off-target 525 pulldown of some nucleosome species, we can only measure the capture of the standards that we 526 actually have spiked into the experiment. If we do not have nucleosome standards available for a 527 potential off-target modification, then we cannot definitively state that the antibody is not 528 capturing that material. In this context, that is likely most important for H3K27me3 pulldowns; 529 though we cannot state this definitively due to the lack of H3K27me2 standards, it is plausible 530 that we are pulling down some amount of H3K27me2 with these IPs, resulting in slightly inflated 531 apparent H3K27me3 HMD. However, this may not be too problematic; H3K27me2 and 532 H3K27me3 are thought to be recognized by many of the same proteins and to have highly 533 similar functions 82 , so the conflation of the two -if present -likely does not pose a significant 534 problem in ascribing biologic function. 535 On a related note, at some loci, the bivalency HMD goes below 0%. In naïve mESCs, 536 8.8% of the promoters have a bivalency HMD below 0%, yet for the vast majority of these 537 promoters (90.8%), the 95% confidence interval error estimate ranges above zero. This is 538 because we employ in silico signal-correction for the bivalency dataset to remove signal that is 539 attributable to H3K9me3. In essence, we can measure the amount of H3K9me3 pulldown in our 540 bivalency ICeChIP dataset due to nucleosome standards employed, and we can separately 541 measure the H3K9me3 HMD by a highly specific IP. We can then a linear combination 542 correction matrix to remove the signal that is attributable to directly measured H3K9me3 at these 543 loci. This method can effectively reduce the impact of modest off-target binding H3K9me3, but 544 at some loci, will result in a subzero apparent HMD due to random sampling of read depth in the 545 two distinct pulldowns employed. 546 Finally, at some sets of gene promoters, the trans-bivalency HMD is shown to be greater 547 than the H3K4me3 or H3K27me3 HMD. This apparent discrepancy has a few possible reasons. 548 First, there is some nuance in the interpretation of HMD in the context of single-target ICeChIP 549 and reICeChIP. A nucleosome has two copies of each of its core histone proteins, including 550 histone H3. This means that there are two possible sites of modification on each nucleosome for 551 for each individual modification; if only one of those sites is modified, then that corresponds to 552 an HMD of 50% because only half the possible modification sites are actually modified. 553 However, this is different for the trans-bivalency HMD; by definition, only one trans-bivalency 554 modification pattern can exist on a given nucleosome at any given time. If two "trans-bivalent" 555 modification patterns existed on the same nucleosome simultaneously, then both H3K4 and both 556 H3K27 residues would be trimethylated -which is symmetric cis-bivalency. As such, if one 557 H3K4 and one H3K27 residue are trimethylated, then 100% of the possible trans-bivalency 558 configurations for the nucleosome of interest are satisfied, meaning that the trans-bivalency 559 HMD will be 100%. However, in this case, the H3K4me3 and H3K27me3 HMDs will only be 560 50% because only half the modifiable residues are actually modified. 561 The other caveat is that symmetrically modified nucleosomes will be pulled down more 562 efficiently than asymmetrically modified nucleosomes due to avidity effects, as can be seen in 563 the pulldown of symmetric vs. asymmetric H3K4me3 and cis-bivalency vs. trans-bivalency 564 (Extended Data Fig. 2), and observed previously 46 . This means that calibration to symmetric 565 nucleosome standards will have a larger denominator in computation of HMD and thereby yield 566 lower apparent HMDs; this can also contribute to the lower apparent HMD of H3K4me3 and 567 H3K27me3 relative to trans-bivalency. Accommodating for this phenomenon would require 568 detailed profiling of asymmetric H3K4me3 (which is currently difficult due to the low quality of 569 H3K4me0 antibodies), asymmetric H3K27me3 (which is not currently possible), and 570 distinguishing between trans-bivalency and cis-bivalency (which is also not currently possible).
However, as noted in Supplementary Note 1, so long as the method of calibration remains 572 consistent, increases in apparent HMD will still correspond to increases in the modification of 573 interest. Whether that increase in the target modification is due to asymmetric modification 574 becoming symmetric or due to new gain of the modification at a previously unmodified locus in 575 an instantaneous subpopulation remains unclear, but in both cases, modification density is still 576 being gained at that locus. As such, even with these caveats, we can still quantitatively compare 577 different datasets to each other as we use consistent calibration standards. 578 579

Supplementary Note 4 580
Intriguingly, the catalytic activity of the EZH2-PRC2 core complex on nucleosome substrates is 581 potentiated by pre-existing H3K27me3 89,90 , yet inhibited by H3K4me3, particularly when 582 symmetric 49,56,59 . Conversely, symmetric H3K27me3 has been reported to modestly inhibit 583 several of the human COMPASS-family complexes by qualitative assays, although only SET1 584 complexes were examined at the nucleosome level 60 . This presents a potential concern for our 585 data -if the enzyme complexes that install these marks are mutually antagonized by the 586 opposing mark, how might the widespread bivalency we observe arise? As the PRC2 effects are 587 well established with detailed quantitative enzymology 49,56,59 , which we recapitulate (data not 588 shown), we deployed more quantitative HMTase assays with a larger panel of relevant 589 nucleosomal substrates to evaluate the COMPASS/SET1B/MLL-family core complexes for 590 allosteric modulation by preexisting marks (Extended Data Fig. 6). The first way we evaluate different models for predicting DEGs is to compute the Bayes 602 Information Criterion (BIC). Though not definitive, this metric estimates whether addition of a 603 parameter to a model improves it more than would be expected from chance alone. When 604 comparing two models, the model with the lower BIC will tend to have more explanatory 605 parameters and/or fewer non-explanatory parameters than the model with the higher BIC. To this 606 end, if BIC increases when a parameter is added, then it can be interpreted that the parameter 607 being added contributes minimal additional explanatory power. Here, we find that adding 608 bivalency to a model increases the BIC, meaning that it is likely (though not definitively) not 609 contributing meaningfully more information in predicting DEG status in this differentiation 610 paradigm. 611 A more definitive way to evaluate model accuracy is to use hold-out cross-validation. In 612 this method, we split the set of all genes into two groups, one with 80% of the genes (the training 613 set) and one with 20% of the genes (the testing set). We then train our GLMs on the training set 614 and use the derived models to predict DEG status in the testing set. Hold-out cross-validation is a 615 highly effective way of testing whether a model is overfit or underfit upon addition or removal of a parameter. If model accuracy increases substantially, then that would suggest the parameter has 617 explanatory power over that provided by the other parameters. Conversely, if model accuracy 618 decreases substantially, then that suggests that the additional parameter causes overfitting. 619 Minimal changes in model accuracy suggest that the additional parameter contributes little to the 620 model over the existing parameters, positively or negatively. 621 There are two metrics we use to test the accuracy of the predictions in the testing set. octamer storage buffer (20 mM Tris-HCl pH 7.5, 2 M NaCl, 1 mM EDTA, 5 mM DTT, 55% 677 glycerol), and stored in -20°C. Concentration of octamer was measured spectroscopically using concentrator flow-through as a blank, ε280nm=44700 M -1 cm -1 , Moct≈108500 g 1 mol -1 . Octamers 679 were visualized using 18% separating (4% stacking) discontinuous Laemmli SDS-PAGE in 680 Mini-Protean gel running system (Bio-Rad) run for 70 minutes at 22mA, 200V max. 681 Asymmetrical H3K4me3, H3K27me3 octamers were done as above with the following 682 differences. Equimolar amounts of histone H2A, H2B, H3 and H4 were mixed in unfolding 683 buffer to the total of 1-2 mg, where 90% of histone H3 was trimethylated on Lys 4 or Lys 27 and 684 remaining 10% were unmethylated and had His6-tag at N-terminus with TEV cleavage site. 685 Octamers were reconstituted overnight by dialysis in 1000 volumes of phosphate refolding buffer 686 (50 mM NaxPO4 pH 7.5, 2M NaCl), at 4°C. Octamers were purified by S200 gel filtration 687 chromatography, and his-tagged octamers were isolated using cobalt-based immobilized metal 688 affinity chromatography Dynabeads magnetic particles. Octamers were incubated with magnetic 689 beads for 10 min at 4°C on rotator, followed by two 1 ml washes (50 mM NaxPO4 pH 7.5, 2 M 690 NaCl, 10 mM imidazole), and eluted with 50 ul of elution buffer (50 mM NaxPO4 pH 7.5, 2 M 691 NaCl, 250 mM imidazole, 1 mM EDTA, 1 mM DTT), the elution step was repeated 6 times, 692 fractions were characterized spectroscopically, pooled, diluted with 1 volume of octamer storage 693 buffer, and stored in -20°C. 694 Asymmetrical trans-bivalent H3K4me3-K27me3 octamers were prepared the same way as 695 symmetrical octamers with the following differences. Histones H2A, H2B, H4 and asymmetric 696 disulfide linked histones H3K4me3 -H3K27me3 were mixed 1.2 : 1.2 : 1 : 0.5 ratio. Remaining 697 steps were done as previously described but no reducing agents were used until octamer particles 698 were formed. 699 All other octamers were obtained from EpiCypher, Inc. 700 20°C. Nucleosome concentration was measured by densitometry of 2% agarose gels, 1x TBE (89 715 mM tris-base, 89 mM boric acid, 2 mM EDTA) run for 30 minutes in 5V/cm electrical field 716 gradient, followed by staining with 1x SYBR Gold (Invitrogen) for >30 minutes. Prior to 717 electrophoresis, nucleosomes were disassembled with 2 M NaCl, roughly 1 pmol of nucleosomes 718 were loaded per well and measured in triplicate against known quantity of free DNA of the same 719 size. For use as ICeChIP standards, the semi-synthetic nucleosomes were diluted to 1 nM 720 concentration using long-term storage buffer (10 mM Na•Cacodylate pH 7.5, 100 mM NaCl, 721 50% Glycerol, 1 mM EDTA, 1x RL Protease Inhibitor Cocktail, 100 µg/mL BSA(NEB)) and 722 stored at -20°C. 723

ICeChIP -input preparation 725
ICeChIP has been done as previously described 46,54 . Briefly, 10 7 -10 8 plate adherent cells were 726 released using Accutase (Millipore), quenched with complete medium and collected (500 rcf, facility. Data analysis was performed as previously described 46 . Briefly, reference genome was 845 modified to contain sequences of semi-synthetic nucleosome barcodes. Reads were mapped to 846 GRCm38/mm10+barcodes reference genome, using Bowtie2 102 , end-to-end, sensitive preset. 847 Subsequently, SAM files were filtered to reject unmapped and unpaired reads, as well as 848 fragments with length > 200bp and Phred quality score < 20. Paired reads were merged into 849 single interval of the fragment. BEDTools 103 were used to calculate bedgraphs of the genome 850 coverage. IP and input bedgraphs of genome coverage were subsequently merged into multiple 851 entry interval file of genome coverage. Number of DNA fragments coming from semi-synthetic 852 nucleosome standards were counted for each barcode and IP efficiency was calculated for each 853 histone mark. 854 where, n is the identity of the n-th barcode assigned to the specific histone mark. 856 where, IP and input refer to the depth of the genome coverage at any given position of IP and 858 input for specific histone mark. An IP efficiency can be interpreted as maximal yield of the IP for 859 a given histone mark, or 100% HMD. However, there is a number of factors that can lead to 860 apparent HMD values greater than 100% including: uncertainty due to random sampling of 861 sequenced fragments, in that case standard deviation is approximately equal to square root of 862 sequencing depth, or off-target capture by antibody, where the off-target histone mark or 863 combination of histone marks have greater IP efficiency than the antibody intended target mark. 864 For all analyses, the HMD averaged over the N+1 and N+2 nucleosomes (taken to be 0 to 865 +400bp into the gene body) was employed as representative of the promoter-this captures the 866 most substantial H3K4me3 and H3K27me3 enrichment. 867 Genomic browser views were made using IGV. Heatmaps and gene ontology analysis was made 868 using Homer Software 104 . Further analysis and sectioning of data was conducted in R using the R 869 code provided in Data Availability. Plots were made using ggplot2 and Microsoft Excel. For the 870 bivalency tracks, in silico antibody off-specificity correction was performed as previously 871 described for H3K9me3 46 . ICeChIP-qPCR was performed with previously described primers 46 . 872 873

Analysis of External Data 874
Bisulfite sequencing data was obtained from GEO series accession number GSE41923, dataset 875 accession numbers GSM1027571, GSM1027572, GSM1027573, and GSM1027574. 876 Methylation count files were obtained for each dataset and lifted to mm10. The average 877 methylation for each promoter was then calculated for the 0 to +400bp region relative to the TSS 878 of Refseq promoters using BEDTools. 879 Bulk RNA-seq data was obtained from GEO series accession numbers GSE108832 and 880 GSE65697, dataset accession numbers GSM2913929, GSM2913930, GSM2913931, 881 GSM1603282, GSM1603283, GSM1603284, GSM1603285, GSM1603286, and GSM1603287.  Fig. 1 | Workflow and evaluation of reICeChIP-seq. (a) Schematic of reICeChIP-seq. The recombinant α-H3K4me3 Fab 304M3-B achieves high affinity by "clasping" the histone tail between two Fab molecules 53 , a binding mode readily achieved by multiple copies of the Fab presented on a bead, but not by the Fab in solution. Thus, protease cleavage not only elutes nucleosomes from the beads but also likely from the Fab complex.  Dom.  Fig. 3b. (e) Distribution of bivalency HMDs at all Refseq promoters in three cell states, zoomed to below 125% HMD. Overall, 99.5% of naïve promoters, 87.3% of primed promoters, and 91.6% of NPC promoters have an HMD below 100%. Full plot in Extended Data Fig. 3a. (f) Metaprofiles of H3K4me3, H3K27me3, and bivalency at promoters identified as bivalent in naïve mESCs (25% HMD threshold), tracked from naïve mESCs to primed mESCs to NPCs. Heatmaps for bivalency are presented in Extended Data Fig. 3f. (g-h) Alluvial plots of dominance and bivalency of genes from (g) naïve mESCs to NPCs or (h) primed mESCs to NPCs. Bivalency [>25% HMD] can be subcategorized into dominance classes by independent ICeChIP for the constituent marks, with H3K27me3 in excess (H3K27me3/H3K4me3 > e 1 ), H3K4me3 dominant (H3K27me3/H3K4me3 < e -1 ), or intermediate ratios (no dominance). (i) Bivalency metaprofiles for gene subsets indicated in panel (h) from -3kb to +3kb relative to the TSS. **** p < 2.2 x 10 -16 .             Endpoints were established at 180 min by kinetic evaluation to be sensitive to difference in activity for this panel. Signal is corrected for background and no nucleosome substrate activity. Error bars represent standard deviation. Fig. 7