H2A.Z.1 Monoubiquitylation Antagonizes BRD2 to Maintain Poised Chromatin in ESCs

Histone variant H2A.Z occupies the promoters of active and poised, bivalent genes in embryonic stem cells (ESCs) to regulate developmental programs, yet how it contributes to these contrasting states is poorly understood. Here, we investigate the function of H2A.Z.1 monoubiquitylation (H2A.Z.1ub) by mutation of the PRC1 target residues (H2A.Z.1 K3R3 ). We show that H2A.Z.1 K3R3 is properly incorporated at target promoters in murine ESCs (mESCs), but loss of monoubiquitylation leads to de-repression of bivalent genes, loss of Polycomb binding, and faulty lineage commitment. Using quantitative proteomics, we ﬁnd that tandem bromodomain proteins, including the BET family member BRD2, are enriched in H2A.Z.1 chromatin. We further show that BRD2 is gained at de-repressed promoters in H2A.Z.1 K3R3 mESCs, whereas BRD2 inhibition restores gene silencing at these sites. Together, our study reveals an antagonistic relation-ship


INTRODUCTION
Pluripotent cells must translate signaling cues into rapid transcriptional responses during development to specify the many cell types in the adult. The histone H2A variant H2A.Z has essential but poorly understood roles in early metazoan development (Faast et al., 2001;van Daal and Elgin, 1992) and murine embryonic stem cell (mESC) differentiation (Creyghton et al., 2008;Hu et al., 2013). In mESCs, H2A.Z is enriched at H3K4me3-marked promoters of active genes and at silent but poised bivalent promoters marked by both H3K4me3 and H3K27me3 (Creyghton et al., 2008;Hu et al., 2013;Ku et al., 2012;Subramanian et al., 2013). Bivalent genes, which encode the majority of developmental regulators in mESCs, are silent yet maintain the capacity to be activated during differentiation (Surface et al., 2010;Voigt et al., 2013). These data suggest that H2A.Z acts as a mo-lecular rheostat for transcriptional output (Subramanian et al., 2015); however, we currently lack a mechanistic understanding of how H2A.Z regulates the balance between gene activation and repression.
H2A.Z facilitates interactions with specific chromatin-associated proteins to affect gene expression in a context-dependent manner (Draker et al., 2012;Fujimoto et al., 2012;Li et al., 2012). Moreover, histone PTMs, including monoubiquitylation, often serve as binding platforms for downstream effectors (Braun and Madhani, 2012). Here, we dissect the function of H2A.Z.1 monoubiquitylation (H2A.Z.1ub) by generating site-specific mutations in the C-terminal PRC1 target lysine residues (K120R, K121R, and K125R; denoted H2A.Z.1 K3R3 ). Using a transgenic mESC system, we show that H2A.Z.1ub is necessary for both proper induction of developmental programs and the maintenance of bivalent chromatin, including PRC2 recruitment. Furthermore, using a quantitative proteomics approach, we find that BRD2, a member of the BET family of transcriptional (legend continued on next page) activators, is enriched in H2A.Z.1 chromatin. BRD2 enrichment is significantly increased at bivalent genes in H2A.Z K3R3 mESCs coincident with gene activation, and BRD2 inhibition restores silencing of bivalent genes. Collectively, our study suggests that H2A.Z.1ub antagonizes BRD2 to maintain the balance between gene activation and repression in response to developmental cues.

RESULTS
H2A.Z.1ub Is Dispensable for mESC Self-Renewal H2A.Z monoubiquitylation (H2A.Zub) is enriched in H3K27me3 containing nucleosomes in mESCs (Ku et al., 2012) and on the inactive X chromosome in human 293T cells (Sarcinella et al., 2007); however, the function of this H2A.Z modification is not known. Using tandem mass spectrometry (MS/MS), we confirmed that RING1B, a catalytic component of PRC1, monoubiquitylates three C-terminal lysines of H2A.Z (K120, K121, and K125) in mESCs (Figures S1A-S1H). We next investigated whether these lysine residues are subject to other modifications. By analyzing raw mass spectrometry data from mESCs (Ku et al., 2012), we confidently assigned signals to the unmodified and ubiquitylated forms of H2A.Z 118-127 . In contrast, we were unable to reliably detect acetylation or methylation of these sites (Figures S1F and S1G; Table S1). We also find that H2A.Zub levels decrease upon mESC differentiation (Figure S1I), coincident with decreased co-localization of H2A.Z and Polycomb complexes at target genes (Creyghton et al., 2008;Ku et al., 2012).
We focused further analysis on the H2A.Z.1 isoform because it is 20-fold more abundant than H2A.Z.2 in mESCs and essential for early development (Faast et al., 2001;Subramanian et al., 2013). To this end, mESC lines were generated that harbor either a doxycycline-inducible H2A.Z.1 transgene with a C-terminal YFP fusion (denoted H2A.Z.1 WT ) or an H2A.Z.1-YFP transgene with the three lysine residues mutated to arginine (K120R, K121R, and K125R; denoted H2A.Z.1 K3R3 ) (Figures 1A and S2A). These mutations resulted in near complete loss of H2A.Z monoubiquitylation in mESCs as determined by immunoblot of histone extracts ( Figures 1B and S2A). We next engineered our H2A.Z.1 WT and H2A.Z.1 K3R3 transgenic mESC lines to constitutively express a short hairpin targeting endogenous H2A.Z.1, resulting in $90% depletion of its mRNA and protein levels with minimal effect on transgene levels ( Figures S2A-S2C).
Similar results were observed with an independent hairpin (Figure S2B). Expression of the RNAi-resistant H2A.Z.1 WT -YFP transgene fully rescues H2A.Z.1 depletion (Subramanian et al., 2013), establishing a system for directly investigating the role of H2A.Z.1ub. Similar to induction of H2A.Z.1 WT , we found that expression of H2A.Z.1 K3R3 in H2A.Z.1-depleted mESCs did not appear to affect morphology, levels of pluripotency markers, cell-cycle dynamics, or response to DNA-damaging agents . These data suggest that H2A.Z.1ub is largely dispensable for mESC self-renewal and may play critical roles in early lineage commitment.
H2A.Z.1ub Is Enriched at Bivalent Promoters H2A.Z is enriched at discrete genomic sites in mESCs, including the majority of H3K4me3-marked transcriptional start sites (TSSs) as well as a small fraction of distal enhancers (Hu et al., 2013;Ku et al., 2012;Subramanian et al., 2013). Our prior work showed that H2A.Z.1-YFP incorporation is highly similar to endogenous H2A.Z as determined by chromatin immunoprecipitation sequencing (ChIP-seq) (Subramanian et al., 2013). Thus, we asked whether loss of H2A.Z.1 monoubiquitylation affected its chromatin incorporation. We found that H2A.Z.1 WT and H2A.Z.1 K3R3 are similarly enriched in chromatin by fractionation of mESCs lysates ( Figure S2H). We also show that H2A.Z.1 WT and H2A.Z.1 K3R3 nucleosomes contain overall similar global histone modification patterns ( Figure S2I). In addition, we found that H2A.Z.1 K3R3 displays similar dynamics to H2A.Z.1 WT in mESCs by fluorescence recovery after photobleaching (FRAP) (Figure S2J), unlike mutation of the divergent acidic patch, which leads to an increase in H2A.Z.1 dynamics (Subramanian et al., 2013). We further demonstrate that conditional ablation of RING1B did not alter H2A.Z incorporation at target genes ( Figure S2K).
We next performed ChIP-seq in H2A.Z.1 WT and H2A.Z.1 K3R3 mESCs using GFP antibodies. Both wild-type and mutant transgenes are similarly enriched at a large set of active (H3K4me3+) and bivalent (H3K4me3+, H3K27me3+) genes (Spearman correlation = 0.972; Figure 1C) as well as at a small number of H2A.Z-enriched enhancers ( Figure S2L). In total, 13,684 and 13,000 genes show H2A.Z.1 enrichment within 2 kb of a TSS in H2A.Z.1 WT and H2A.Z.1 K3R3 mESCs, respectively. As examples, H2A.Z.1 K3R3 occupies both active promoters including Thra as well as bivalent promoters including Nestin and genes (C) Heatmap of ChIP-seq in H2A.Z.1 WT and H2A.Z.1 K3R3 mESCs displays similar incorporation at active and bivalent target genes. H3K4me3 and H3K27me3 (Wamstad et al., 2012) are included for comparison. Heatmaps are centered on transcriptional start sites (TSSs) of all genes ranked by H2K27me3 and extend ±2 kb. Bivalent gene classification is from Subramanian et al. (2013). (D) Sequential ChIP using GFP followed by GFP or H2Aub (also recognizes H2A.Zub) antibodies. % input is calculated using 2 (Cp(Input)-Cp(ChIP)) *(%of total input). Error bars represent SD of triplicate reactions. (E) Distribution of the log2 fold change in expression in H2A.Z.1 K3R3 mESCs was plotted using cumulative density function of log 2 (H2A.Z.1 K3R3 /H2A.Z.1 WT ). Black trace represents all genes with RPKM of at least 1 in any sample and five unique reads in each sample. Active (yellow trace) and bivalent (red trace) gene classes that pass threshold were plotted independently. (F) Box plots represent the log2 fold change in expression of either H2A.Z.1 K3R3 or H2A.Z.1 KD mESCs relative to H2A.Z.1 WT mESCs at all (13,857 genes), active (11,010), bivalent (1,611), or K4me3-/K27me3-(61) genes. Classification of bivalent and active genes is from Subramanian et al. (2013). Center line represents median value; box limits indicate the 25 th and 75 th percentiles; whiskers extend 1.5 times the interquartile range from the 25 th and 75 th percentiles, and outliers are represented by dots. (G) Volcano plot of all genes that passed detection threshold were plotted by p value versus change in expression. (H) Top ten enriched categories using PANTHER Biological Process database are displayed for upregulated genes compared to all genes. See also Figures S1 and S2. associated with the HoxA cluster, in a pattern similar to H2A.Z.1 WT ( Figure S2M).
Our ChIP-seq data show that H2A.Z.1 K3R3 occupies promoters similar to endogenous H2A.Z.1 WT ; however, the localization pattern of H2A.Z.1ub is not known due to the lack of specific antibodies that distinguish H2Aub and H2A.Zub ( Figure S2I). To address this limitation, we performed sequential ChIP by first enriching for H2A.Z.1 nucleosomes in H2A.Z.1 WT and H2A.Z.1 K3R3 mESCs using a GFP antibody followed by re-ChIP with an antibody that recognizes H2A/Zub. H2A.Z.1 K3R3 mESCs that lack H2A.Z.1ub were used as a negative control. Sequential ChIP demonstrates that H2A.Z.1ub is largely enriched at bivalent promoters and low at active promoters ( Figure 1D). As a control, re-ChIP with GFP antibodies shows that both H2A.Z.1 WT and H2A.Z.1 K3R3 are incorporated similarly ( Figure 1D). Collectively, these data suggest loss of H2A.Z.1 monoubiquitylation does not affect global H2A.Z.1 incorporation or dynamics in mESCs.

H2A.Z.1ub Regulates Developmental Gene Expression Programs
Given that the majority of enriched regions map to TSSs, we asked whether loss of H2A.Z.1 monoubiquitylation affects gene expression. We profiled the transcriptome of H2A.Z.1 K3R3 mESCs by RNA sequencing (RNA-seq) and found that bivalent genes were expressed at higher levels compared to H2A.Z.1 WT controls (p < 2.2 3 10 À16 , Kolmogorov-Smirnov test, two-sided) ( Figure 1E). In contrast, levels of active genes were largely unaffected in H2A.Z.1 K3R3 mESCs, similar to H2A.Z.1 depletion ( Figures 1E and 1F) (Hu et al., 2013;Subramanian et al., 2013). Because bivalent genes are lowly expressed (Mikkelsen et al., 2007), we applied stringent threshold criteria to reduce false positives in our analysis (expression of R1 reads per kilobase per million [RPKM] in at least one sample and at least five reads in every sample). A total of 9,667 active genes and 816 bivalent genes passed these criteria. Specifically, 370 and 107 of these genes are up-and downregulated, respectively, in H2A.Z.1 K3R3 mESCs using a cutoff of 1.5-fold change and p value % 0.05 ( Figure 1G). The set of 370 upregulated genes significantly overlaps the 816 bivalent genes (p < 2.65 3 10 À87 , hypergeometric test).
Using the gene annotation tool PANTHER (Huang et al., 2009a(Huang et al., , 2009b, we found that the upregulated genes function in cell communication, signaling, and development ( Figure 1H). In contrast, we did not find significant overlap with bivalent genes or enriched Gene Ontology (GO) terms among the downregulated genes. Expression changes were validated using a second independent H2A.Z.1 hairpin ( Figure S2N). Notably, although H2A.Z.1 depletion also leads to de-repression of bivalent genes (Hu et al., 2013;Subramanian et al., 2013), we observed overall higher expression of these genes in H2A.Z.1 K3R3 mESCs (Figure 1F; ANOVA, p < 0.0001), suggesting that H2A.Z incorporation is critical for gene activation in response to developmental signals.

H2A.Z.1ub Is Required for Proper mESC Differentiation
Because the precise regulation of bivalent genes is key for proper lineage commitment (Subramanian et al., 2015), we next examined the differentiation capacity of H2A.Z.1 K3R3 mESCs by allowing cells to aggregate into embryoid bodies (EBs), a process that leads to multi-lineage differentiation similar to the gastrulating embryo (ten Berge et al., 2008). Induction of the H2A.Z.1 WT transgene restores proper mESC differentiation as evidenced by appropriate expression of germ layer markers (Figures 2A-2C). In contrast, EBs generated from H2A.Z.1 K3R3 mESCs failed to undergo multi-lineage differentiation, as evidenced by H&E-stained sections displaying distinct differences in tissue representation relative to H2A.Z.1 WT (Figure 2A). In particular, H2A.Z.1 K3R3 EBs lack neuroepithelial structures and failed to activate the neural marker TUJ1 compared to H2A.Z.1 WT EBs ( Figure 2B). Additionally, we found that genes involved in neuroectoderm lineages (e.g., Sox1, Sox3, and Pax6) are not properly induced in H2A.Z.1 K3R3 EBs, whereas mesendodermal genes are highly expressed (Figures 2C and S3A). Hyperactivation of WNT signaling during EB formation promotes cardiac mesoderm and inhibits neuroectodermal differentiation (ten Berge et al., 2008). We found that WNT signaling genes including Wnt3A and Axin2 were upregulated in H2A.Z.1 K3R3 mESCs ( Figures S3B and S3C). H2A.Z.1 K3R3 EBs also display increased nuclear b-catenin levels, reduced binding of TCF3 at WNT target genes, and increased activity of a WNT signaling reporter compared to H2A.Z.1 WT controls ( Figures  S3D-S3F). Furthermore, treatment with the WNT antagonist KY02111 partially rescued their differentiation defects as shown by activation of the neuroectodermal genes Sox3 and Pax6 (Figure 2C). These results suggest that H2A.Z.1ub is necessary for coordinating specific transcriptional outputs in response to developmental cues.

H2A.Z.1ub Is Necessary for Maintenance of Bivalent Chromatin
The activity of PRC1 regulates bivalent gene expression, potentially through RNAPII pausing (Stock et al., 2007), and/ or recruitment of PRC2 Cooper et al., 2014;Kalb et al., 2014). However, we found that expression of highly paused genes in mESCs (Rahl et al., 2010) is not significantly affected in H2A.Z.1 K3R3 mESCs ( Figure 3A), consistent with the markedly lower levels of RNAPII at bivalent promoters compared to active genes ( Figure S4A) (Williams et al., 2015). Thus, we next asked whether H2A.Z.1ub, like H2Aub, contributes to maintenance of PRC2 at bivalent genes in mESCs. We performed ChIP-qPCR in H2A.Z.1 K3R3 mESCs and observed a significant decrease in the enrichment of the PRC2 component SUZ12 at bivalent promoters, as well as its corresponding mark H3K27me3 ( Figures 3B and 3C). In addition, PRC1 components including RING1B as well as KDM2B, a key factor responsible for recruitment of non-canonical PRC1 to CpG islands He et al., 2013;Wu et al., 2013) were reduced at bivalent promoters in H2A.Z.1 K3R3 mESCs ( Figures 3D and S4B). Conversely, loss of monoubiquitylation did not affect H3K27ac, H4ac, or H2A.Zac levels and led to a slightly increase in levels of H3K4me3, modifications associated with transcriptional activity ( Figures 3E and S4C-S4E). Together, these results suggest that H2A.Z.1ub plays a key role in maintaining bivalent chromatin and gene silencing in mESCs.

H2A.Z.1ub Interacts with Specific Proteins in mESCs
Histone post-translational modifications often function to recruit downstream regulators (Braun and Madhani, 2012), and specifically, ubiquitin moieties often mediate downstream interactions with regulatory factors. Given that H2A.Z interacts with a distinct set of binding partners compared to H2A (Draker et al., 2012;Fujimoto et al., 2012), we reasoned that H2A.Z.1ub might mediate interactions with specific chromatin-associated proteins. To test this idea, we developed an approach that combines SILAC (stable isotope labeling of amino acids in cell culture) (Ong et al., 2002) with immunoprecipitation (IP) and mass spectrometry to quantitatively assess changes in chromatin associated protein interactions (SILAC-IP) ( Figure 4A). We first compared enriched proteins in mESC lines containing either an H2A.Z.1-YFP or H2A-YFP transgene (Subramanian et al., 2013) (Figures S5A  and S5B; Table S2). Proteins were labeled by culturing cells in isotope containing media prior to IP, resulting in greater than 93% labeling efficiency as measured by the incorporation of the heavy isotopes. SILAC-IP was performed on the soluble fraction of micrococcal nuclease-treated nuclei. Our approach allows for direct comparison between target proteins within one experiment and internally controls for nonspecific as well as YFP-associated interactions. To confidently identify H2A.Z.1specific interactions, we required representation of multiple peptides in replicate samples from paired isotope swaps (e.g., H2A.Z.1-YFP heavy and H2A-YFP light) ( Figure S5C). Differential enrichment was ranked by the fold change of ratio counts between paired samples and by a modified t-statistic for significance of reproducibility between replicates (Mertins et al., 2013;Sancak et al., 2013) (Figure 4B; Table S2). Validating our approach, we identified nearly all members of the H2A.Zspecific ATP-dependent deposition complex SRCAP among the H2A.Z.1-enriched proteins (Ruhl et al., 2006) (Figures 4B and S5D).
Given that many protein interactions are facilitated by specific domains, we asked whether particular motifs were enriched among the set of H2A.Z.1-interacting proteins (Franceschini et al., 2013). Our analysis identified the bromodomain as the most highly enriched motif (p value < 2e-3). Over 40 bromodomain-containing proteins exist in the mouse genome, whereas only 8 harbor a tandem bromodomain motif (Filippakopoulos et al., 2012). Five of these eight proteins are both expressed in mESCs and enriched in H2A.Z.1 chromatin, including two BET proteins (BRD2 and BRD3) and three proteins that harbor an additional WD40 motif (BRWD2, BRWD3, and PHIP) ( Figures  4B and 4C). Because BRD2 is the most highly expressed of the tandem bromodomain members in mESCs and it is necessary for murine development (Shang et al., 2009) (Figure S5E), we focused on validating its interaction with H2A.Z.1 using two independent approaches. First, we found that BRD2 is significantly   (Table S3). Second, BRD2 shows significantly greater association with H2A.Z.1 in chromatin compared to H2A (p value < 0.1, t test) by bioluminescence resonance energy transfer (BRET) (Deplus et al., 2013) (Figures S5F and S5G). In contrast, the BET protein BRD4 exhibited neither differential enrichment by SILAC-IP nor differential interaction by BRET ( Figure S5G)  We then asked if H2A.Z K3R3 alters the association with H2A.Zenriched proteins using our unbiased SILAC-IP approach. We found 212 differentially enriched proteins between H2A.Z.1 K3R3 and H2A.Z.1 WT SILAC-IPs ( Figure 4D; Table S4). Of these, 123 proteins including BRD2 are significantly enriched in H2A.Z.1 K3R3 chromatin compared to H2A.Z.1 WT SILAC-IPs. In contrast, other tandem bromodomain proteins as well as SRCAP components showed little difference in enrichment (Figures 4D-4F). BRD2 transcript and protein levels are similar between wildtype and mutant cell lines indicating that differential enrichment is not due to higher levels of BRD2 in H2A.Z.1 K3R3 mESCs (Figures S5H and S5I). Together, our quantitative proteomic analysis demonstrates that H2A.Z.1ub influences specific protein interactions including with BRD2.

H2A.Z.1ub Inhibits BRD2 Recruitment to Bivalent Genes
We next performed ChIP-seq to analyze the global enrichment patterns of BRD2 in mESCs. BRD2 is enriched at 11,536 genomic regions in wild-type mESCs. Of these regions, 58.7% fall within 2 kb of a known TSS, leading us to focus on the role of BRD2 at promoters (Figures 5A and 5B). In total, BRD2 localizes to the promoters of 6,856 genes of which 84% also harbor H2A.Z.1 peaks. Moreover, BRD2 enrichment at promoters coincides with active chromatin marks including H2A.Zac, H3K27ac and H3K4me3 as well as RNAPII ( Figures 5C and S6A), consistent with the role of bromodomain-containing proteins in gene activation (Draker et al., 2012;LeRoy et al., 2008). In contrast, BRD2 is largely absent from bivalent genes, e.g., the HoxA locus ( Figures 5C and 5D). Indeed, genes co-enriched with both H2A.Z.1 and BRD2 are more highly expressed than H2A.Z-enriched genes that lack BRD2 (p value < 1e-142, t test) ( Figure 5E). These results suggest a role for BRD2 with H2A.Z.1 at active gene promoters.
Our proteomics data suggest that H2A.Z.1ub inhibits BRD2 recruitment. To test this idea, we performed ChIP-seq and observed a significant increase in BRD2 enrichment at bivalent promoters in H2A.Z.1 K3R3 compared to H2A.Z.1 WT mESCs (Figure 5F; median fold change 1.50, Wilcox test p < 2.4e-118, compared to all genes), whereas active genes showed only a slight gain ( Figure 5F; median fold change 1.09, Wilcox test p < 0.05, compared to all genes). Moreover, BRD2 levels were not significantly changed at either enhancers or super enhancers (two-sided t test, median fold change enhancers 0.97, super enhancers 0.95) ( Figure S6B). For example, BRD2 shows increased enrichment at the bivalent gene Mesp1 as well as across the HoxA gene cluster in H2A.Z.1 K3R3 mESCs ( Figure S6C), whereas active genes including Hira and Cul1 show no apparent difference in BRD2 levels ( Figure S6D). These results were independently validated by ChIP-qPCR ( Figure S6E). In total, we find that 1,128 genes display a 2-fold or greater change in BRD2 enrichment in H2A.Z.1 K3R3 mESCs across biological replicates (Table S5). GO analysis revealed that genes gaining BRD2 have roles in signal transduction and developmental processes ( Figure 5G) and showed a significant increase in expression compared to all genes (median change 1.3-fold, Wilcox test p < 9.2e-17). In contrast, genes that displayed >2-fold depletion of promoter BRD2 were largely unaffected in H2A.Z.1 K3R3 mESCs ( Figure 5H). BRD4, another BET family member, plays a key role in regulating gene expression and cell identity in ESCs and in many cancers (Asangani et al., 2014;Di Micco et al., 2014;Liu et al., 2014;Wu et al., 2015). In contrast to BRD2, BRD4 enrichment was more modestly affected at bivalent promoters in H2A.Z K3R3 mESCs (median fold change 1.29, Wilcox test p < 1.9e-17 compared to all genes) (Figures S6H and S6I). Moreover, BRD4 was not differentially enriched in our SILAC-IP or BRET analysis ( Figures 4B, 4D, 4F, and S5G; Tables S2 and S4). Notably, BRD2 and H2A.Z.1 both display a highly similar bimodal distribution at promoters by ChIP-seq, whereas BRD4 shows a distinct single peak ( Figure S6F). Furthermore, BRD4 is more highly enriched at enhancers, particularly at superenhancer clusters ( Figure S6G) (Whyte et al., 2013), whereas BRD2 and H2A.Z.1 are most highly enriched at promoters (Figures S6F and S6G). These data suggest H2A.Z.1 coordinates with BRD2 at promoters and that BET family members function independently to regulate gene expression programs in mESCs.

H2A.Z.1ub Inhibits BRD2 Recruitment to Bivalent Promoters
To test whether H2A.Z.1ub prevents BRD2 recruitment at bivalent promoters, thereby inhibiting gene activation, we first treated mESCs with the pan-BET inhibitor JQ1 (Liu et al., 2014;Wu et al., 2015). After 24 hr of treatment with 100 nM JQ1, the increased expression of bivalent genes in H2A.Z.1 K3R3 mESCs was restored to basal levels ( Figure 6A). H2A.Z K3R3 mESCs exhibit de-repression of early mesoendoderm genes including Brachyury and Mesp1, whereas neuronal genes such as Pax6 and Sox1 are largely unaffected ( Figure 2C). Consistent with these observations, Pax6 and Sox1 were largely unaffected by addition of JQ1. We next investigated whether loss of gene silencing in H2A.Z.1 K3R3 mESCs can be rescued by direct BRD2 inhibition. Transfection with specific small interfering RNAs (siRNAs) resulted in significant BRD2 depletion after 48 hr ( Figure 6B). Upon BRD2 depletion, bivalent genes showed a marked decrease in expression in H2A.Z.1 K3R3 mESCs, whereas genes unaffected by JQ1 treatment remained unchanged ( Figure 6C). Similar results were observed with an independent BRD2 siRNA (data not shown). In contrast, BRD4 depletion led to a slight increase in bivalent gene expression in H2A.Z.1 WT mESCs and to a dramatic increase in expression in H2A.Z.1 K3R3 mESCs ( Figures 6D and 6E). These results are consistent with a direct role for BRD2 in regulating gene promoters, whereas BRD4 plays a pivotal role in regulating superenhancer-associated pluripotency genes such as Oct4 (Di Micco et al., 2014;Whyte et al., 2013).
H2A.Z.1ub is critical for gene silencing and maintenance of PRC2 at bivalent genes, leading us to ask whether BRD2 antagonizes PRC2 to reinforce gene activation. To test this idea, we analyzed SUZ12 occupancy at a subset of target genes upon JQ1 treatment or BRD2 depletion in H2A.Z.1 WT and H2A.Z.1 K3R3 mESCs. BET inhibition with 100 nM JQ1 for 24 hr led to an increase in SUZ12 enrichment at bivalent genes in H2A.Z.1 WT , with only a slight increase in H2A.Z.1 K3R3 mESCs, whereas active genes were largely unaffected ( Figure 6F). Similarly, direct BRD2 inhibition resulted in increased SUZ12 occupancy at bivalent genes in H2A.Z WT , but not in H2A.Z.1 K3R3 , mESCs, whereas SUZ12 levels at active genes were largely unaffected ( Figure 6G). Taken together, these results suggest H2A.Z.1ub is required for SUZ12 binding and functions to maintain the balance between gene silencing and activation of developmental programs in mESCs ( Figure 6H).

DISCUSSION
H2A.Z is incorporated at both active and bivalent promoters in mESCs and appears to regulate contrasting gene expression states. While H2A.Z acetylation is widely associated with gene (G) GO analysis of genes gaining BRD2 at their promoter as determined by DAVID (Huang et al., 2009a(Huang et al., , 2009b. (H) Box and whisker plot of the median fold change in expression of genes gaining BRD2 >2-fold at promoters (p value 6.3e-11, un-paired t test) and losing BRD2 (not significant). Box and whiskers defined as in (E). See also Figure S6.   activation from yeast to mammals (Bruce et al., 2005;Hu et al., 2013;Ku et al., 2012;Millar et al., 2006), how this variant contributes to gene silencing is poorly understood. Using a defined system to directly test the function of the three C-terminal PRC1 monoubiquitylated lysine residues, we show that H2A.Z.1ub is necessary for both the maintenance of bivalent chromatin in mESCs and for the appropriate induction of developmental programs.
Non-canonical PRC1 containing KDM2B catalyzes H2A ubiquitylation and appears necessary for PRC2 recruitment and transcriptional repression both in vitro and in vivo; however, a direct role for H2Aub in this process was not established in these studies Cooper et al., 2014;Kalb et al., 2014). We show that loss of H2A.Z.1ub leads to disruption of poised, bivalent genes and to faulty mESC differentiation, suggesting that H2A.Z.1ub contributes to transcriptional repression of lineage programs in mESCs. The observation that mutation of the catalytic activity of PRC1 components RING1A/B does not affect chromatin compaction at target regions despite the partial loss of gene repression and H3K27me3 enrichment Eskeland et al., 2010;Illingworth et al., 2015) suggests neither H2Aub nor H2A.Zub contributes directly to chromatin compaction. Notably, RING1B catalytic mutants progress further in development than RING1B knockout embryos, similar to the partial defects in lineage commitment we observe upon differentiation of H2A.Z K3R3 mESCs. Together, these data suggest that RING1B-mediated monoubiquitylation is also an important PRC1 activity. In contrast, recent work showed that the ubiquitylation activity of the RING1A/B homolog, SCE, is not required for gene repression in Drosophila (Pengelly et al., 2015). Thus, dissecting how histone variant incorporation and PRC1 contribute to gene repression and chromatin compaction will lead to a better understanding of how these pathways coordinate developmental programs and may explain the different observations in Drosophila and mammalian systems.
Loss of H3K27me3 alone is not sufficient to confer transcriptional activation (Marks et al., 2012), suggesting that activation is mediated through a distinct mechanism. Our data suggest that modulation of H2A.Z.1 by PTMs is critical for its interaction with downstream activators such as BRD2. In support of this model, studies show that BRD2 co-localizes with H2A.Z.2 in melanoma cell lines to activate cell-cycle genes (Vardabasso et al., 2015) and that H2A.Z is necessary to recruit BRD2 to activate target genes in response to hormone stimulation in LNCaP cells (Draker et al., 2012). Thus, it is possible that the H2A.Z is monoubiquitylated at the poised hormone-responsive genes prior to stimulation.
How BRD2 drives transcription remains unclear, as it lacks the p-TEFb interacting motif of BRD4 (Jang et al., 2005).
In vitro, BRD2 directly promotes transcription on a chromatinized template by acting as a chromatin remodeler (LeRoy et al., 2008). Additionally, BRD2 can bind both TBP and Mediator, both of which are normally absent from bivalent promoters (Peng et al., 2007;Pinz et al., 2015). Thus, future studies will be directed toward understanding how H2A.Z.1ub inhibits gene activation through antagonism of BRD2. More broadly, H2A.Z homologs in both yeast and C. elegans coordinate with BET proteins to regulate gene expression states (Shibata et al., 2014;Zhang et al., 2005), suggesting that coordination between H2A.Z and BRD2 is an evolutionarily conserved mechanism for controlling transcriptional output at inducible genes. Together, our work establishes a previously unknown role for H2A.Z.1ub in maintaining the transcriptional balance of developmental programs during early lineage commitment.

EXPERIMENTAL PROCEDURES
Growth of mESCs V6.5 (129SvJae and C57BL/6; male) mESCs were plated with irradiated murine embryonic fibroblasts (MEFs) and cultured using standard conditions on gelatinized tissue culture plates as described (Subramanian et al., 2013). For WNT signaling inhibition, mESCs were treated with DMSO or the WNT signaling inhibitor KY02111 (10 mM) for 48 hr and then aggregated to form EBs in the absence of leukemia inhibitory factor (LIF) and in the presence of the inhibitor. Cells were treated with 100 nM JQ1 or DMSO as a control. mESCs were transfected with specific siRNAs to Brd2 and Brd4 (Origene) using 10 ml DharmaFect 1 (GE Healthcare). Expression was assayed 48 hr after the start of siRNA transfection. Detailed protocols are included in Supplemental Experimental Procedures.
Embryoid Body Differentiation mESCs were aggregated to form EBs by plating them in Corning Ultra-Low Attachment Tissue Culture Plates (Corning) at a density of 100,000 cells/ml in mESC media lacking LIF. Media was changed every other day, and EBs were cultured for the indicated number of days.

Histone Extracts and Immunoblot
Histone extracts were prepared as described previously (Subramanian et al., 2013). For immunoblot analysis, samples were resolved on SDS-PAGE gels, proteins were transferred to a PVDF membrane, blocked with 5% milk in PBST (0.1% Tween-20 in PBS [pH 7.4]), and blotted overnight with primary antibody in PBST. The presence of the antigen was detected by horseradish peroxidase (HRP)-conjugated secondary antibody. Antibodies listed in Supplemental Experimental Procedures.
(E) mRNA expression for H2A.Z.1 bivalent and active target genes after treatment with either a negative control siRNA or BRD4 siRNA, normalized as above. (F) SUZ12 ChIP-qPCR after treatment of H2A.Z.1 WT or H2A.Z.1 K3R3 mESCs with either 100 nM JQ1 or DMSO. Enrichment values are normalized to two gene desert regions. Error bars represent SE of a triplicate set of experiments. (G) SUZ12 ChIP-qPCR after treatment of H2A.Z.1 WT or H2A.Z.1 K3R3 mESCs with either BRD2 siRNA or negative control siRNA. Enrichment values were calculated as above.
(H) Model for BRD2 recruitment to bivalent genes. BRD2 is localized to active promoters that lack H2A.Z.1ub (indicated by red circles), whereas BRD2 is largely absent from bivalent genes enriched for monoubiquitylated H2A.Z.1. H2A.Z.1 K3R3 incorporation leads to BRD2 recruitment and gene activation, which further antagonizes PRC2 at promoters of developmental genes.
RNA Isolation, Real-Time qPCR, and RNA-Seq Total RNA was extracted using TRIzol (Invitrogen) and reverse transcribed using SuperScript III (Invitrogen) or M-MLV reverse transcriptase (Invitrogen) and random hexamers according to manufacturer protocols. qPCR reactions were performed with SYBR Green (KAPA Biosystems) and primers listed in Table S6. Relative mRNA levels were determined in triplicate for each transcript using the manufacturer's software (Advanced Relative Quantification with Roche Lightcycler 480 Software Version 1.5) using Tubb5 transcript levels for normalization. RNA-seq libraries were prepared as described elsewhere (Subramanian et al., 2013). Detailed protocols are described in the Supplemental Experimental Procedures.

Chromatin Immunoprecipitation
ChIP was performed with 10-25 million cells as described previously (Wamstad et al., 2012). Detailed protocols including site-specific qPCR and highthroughput sequencing as well as antibodies are described in the Supplemental Experimental Procedures.

SILAC-Based Immunoprecipitation and Mass Spectrometry
Cells were grown for 1 week in SILAC medium containing either R 0 K 0 or R 10 K 8 . Immunoprecipitation was performed using 10 9 cells for each condition. Native ChIP was performed on a modified version as described previously (Umlauf et al., 2004). Eluted proteins were separated by SDS gel. In-gel digestion and mass spectrometry were conducted essentially as described previously (Sancak et al., 2013). MS data were analyzed as in Cox and Mann (2008) and Cox et al. (2011). Detailed explanation can be found in Supplemental Experimental Procedures.

ACCESSION NUMBERS
The accession numbers for the sequencing data reported in this article are GEO: GSE53208 and GEO: GSE67944. The original mass spectra may be downloaded from MassIVE (http://massive.ucsd.edu) using the identifier MSV000079429. The data are accessible at ftp://MSV000079429:a@ massive.ucsd.edu.

SUPPLEMENTAL INFORMATION
Supplemental Information includes Supplemental Experimental Procedures, six figures, and six tables and can be found with this article online at http:// dx.doi.org/10.1016/j.celrep.2015.12.100.

ACKNOWLEDGMENTS
H2A.Z.1 WT -and H2A.Z.1 K3R3 -GFP vectors were a kind gift from Peter Cheung, Ontario Cancer Institute. Ring1b fl/À ; CreER T2 mESCs were a kind gift from Maarten van Lohuizen. GFP antibody was kindly provided by Iain Cheeseman. JQ1 was a kind gift from Jay Bradner. We thank Charles Lin and members of the MIT BioMicro Center especially Fugen Li, Vincent Butty, Jie Wu, and Stuart Levine for assistance with data analysis. We thank Danette Daniels and Jacqui Mendez from Promega for technical advice and reagents for the BRET assay. We thank Tanya Svinkina for assistance processing IP samples. Raw LC-MS data from Ku et al. (2012) were used with permission of Bradley Bernstein. We are also grateful to members of the L.A.B. lab for helpful and stimulating discussions. This work was supported in part by the NIH Pre-Doctoral Training Grant T32GM007287 (L.E.S. and P.A.F), the MIT School of Science Fellowship in Cancer Research (P.A.F.), core grant award P30-CA14051 from National Cancer Institute of the NIH, NHLBI Bench to Bassinet Program (U01HL098179), and the Massachusetts Life Sciences Center to L.A.B.