Dzip3 regulates developmental genes in mouse embryonic stem cells by reorganizing 3D chromatin conformation

In mouse embryonic stem (mES) cells, ubiquitylation of histone H2A lysine 119 represses a large number of developmental genes and maintains mES cell pluripotency. It has been suggested that a number of H2A ubiquitin ligases as well as deubiquitylases and related peptide fragments contribute to a delicate balance between self-renewal and multi-lineage differentiation in mES cells. Here, we tested whether known H2A ubiquitin ligases and deubiquitylases are involved in mES cell regulation and discovered that Dzip3, the E3 ligase of H2AK119, represses differentiation-inducible genes, as does Ring1B. The two sets of target genes partially overlapped but had different spectra. We found that Dzip3 represses gene expression by orchestrating changes in 3D organization, in addition to regulating ubiquitylation of H2A. Our results shed light on the epigenetic mechanism of transcriptional regulation, which depends on 3D chromatin reorganization to regulate mES cell differentiation.

which, like Ring1B, is an E3 ligase targeting H2AK119. However, even double knockout of Ring1A and Ring1B does not lead to the complete removal of ubH2A around the transcription start site (TSS) 6 . This result suggests the existence of additional site-specific factors that are involved in mediating ubH2A modifications and repressing a specific set of genes.
Dzip3 is known to be an E3 ligase targeting H2AK119 and was first identified as an RNA-binding RING-dependent ubiquitin ligase 7 . In C2C12 cells, Dzip3 shows nuclear localization and modulates specific histone modifications, rather than exerting global effects through interactions with Nco-R, HDAC1, and HDAC3 8 . It is of great interest whether Dzip3 contributes to gene regulation in ES cells.
For ES cells to differentiate, they need to switch on developmental genes by regulation of an epigenetic pathway involving deubiquitylases, which remove ubiquitin moieties from H2AK119. Until now, five deubiquitylases (USP3 9 , USP16 10 , USP21 11 , USP22 12 , and 2A-DUB [also known as MYSM1] 13 ) have been reported. Recently, USP16 and USP22 were shown to regulate the differentiation process in mouse ES (mES) cells 14,15 . However, whether other deubiquitylases also contribute to the gene regulation of ES cells has not been determined.
In this study, we tested whether H2A ubiquitin ligases and deubiquitylases are involved in the regulation of pluripotency and the differentiation process in mES cells and demonstrated that Dzip3 regulates developmental genes in mES cells by reorganizing 3D chromatin conformation.

Results
Dzip3 regulates developmental genes. ES cells have a distinct morphology, a small size with little cytoplasm, and tightly packed colonies with round or polygonal borders. Upon differentiation, the cells typically expand and flatten out, often losing their tightly packed appearance, which leads to expansion of the colony. To identify H2A ubiquitin ligase or deubiquitylase activities, which play important roles in mES cells, we examined changes in ES cell morphology after performing siRNA knockdown (KD) of eight proteins: Ring1B (also known as Rnf2), Rnf8, and Dzip3, which are ubiquitin ligases, and USP3, USP16, USP21, USP22, and Mysm1, which are deubiquitylases [16][17][18] .
First, we focused on morphological changes under pluripotent conditions (serum + LIF) 19 . KD of Ring1B and Dzip3 resulted in a decrease in tight packing of the cells and led to an expansion of the colony ( Fig. 1A and Supplementary Fig. 1b). KD of Dzip3 and Ring1B was confirmed by RT-qPCR analysis and western blotting ( Fig. 1B and Supplementary Fig. 2a), and the amount of H2A ubiquitylation (ubH2A) in the cell was determined by western blotting. The total amount of ubH2A in the cell decreased after Ring1B KD but not after Dzip3 KD (Fig. 1C), which suggests that Dzip3 functions by modulating a specific histone modification in the promoter region, rather than by global effects.
To determine the effect on gene transcription of Dzip3 KD in mES cells, we investigated the gene expression of pluripotency markers by RT-qPCR (Fig. 1D). However, transcriptional changes in these markers were not apparent. The transition from a pluripotent stem cell to a committed cell type is accompanied by stable silencing of pluripotency genes and activation of lineage-specific genes. We therefore next investigated lineage-specific gene expression after Dzip3 KD in mES cells. The transcriptional levels of several lineage-specific genes, such as Rhox6, Stra8, T (encoding the Brachyury protein), Acta1, and Eomes, were upregulated by Dzip3 KD. Significant numbers of genes were upregulated by simultaneous KD of Dzip3 and Ring1B (Fig. 1E). These results suggest that Dzip3, together with Ring1B, represses lineage-specific gene expression in mES cells.
Whole-transcriptome sequencing analysis shows overlap of Dzip3 and Ring1B target genes. To obtain insight into whole-transcriptome effects after Dzip3 KD in mES cells, we performed RNA-seq analysis. Lineage-affiliated gene expression is strongly repressed when mES cells are cultured in the presence of the two inhibitors PD0325901 and CHIR99021 (both components of 2i medium) plus LIF compared with culturing in the presence of fetal bovine serum (FBS) plus LIF. These results indicate that, in the presence of FBS plus LIF, the mES cell is in a metastable state rather than in a state exhibiting the inherent properties of pluripotent cells 20,21 . We confirmed that, after Dzip3 KD, the changes in colony morphology and expression levels of pluripotency markers were comparable between ES cells cultured with 2i medium plus LIF and serum plus LIF (data not shown). Therefore, we decided to perform whole-transcriptome sequencing analysis with cultured mES cells under 2i-medium-plus-LIF conditions.
Whole-transcriptome sequencing revealed that KD of either Dzip3 or Ring1B altered expression levels of a significant number of genes, but the genes affected by the former were somewhat different from those affected by the latter ( Fig. 2A and Supplementary Fig. 3a). However, there was a substantial set of genes whose expression was affected by both proteins. A Venn diagram (Fig. 2B) represents the overlap between the genes de-repressed by Ring1B KD and those by Dzip3 KD, showing that about one fifth of the genes whose expression was derepressed by Dzip3 KD were also derepressed by Ring1B KD (Fig. 2B). To gain insight into the molecular functions of genes that were repressed by both Dzip3 and Ring1B, we conducted gene ontology (GO) classification analyses, and the results suggest that the common targets of Dzip3 and Ring1B are mostly participants in developmental processes (Fig. 2C).
Some developmental genes that were upregulated by Dzip3 KD or Ring1B KD were selected based on RNA-seq results and validated by RT-qPCR (Fig. 2D). Neurod1 and Neurog1 are expressed in neuroectoderm, Rhox6 is involved in primordial germ cell differentiation, and Cdh2 is a neural cadherin. Neurod1 and Cdh2 were mainly upregulated by Dzip3 KD, while Neurog1 and Rhox6 were upregulated by both Dzip3 and Ring1B KD. N2B27 differentiation conditions caused a typical flattened morphology of ES cell colonies, but no noticeable changes were evident due to Dzip3 KD under differentiation-inducing conditions (data not shown). However, we found that, for some genes, Dzip3 KD not only derepressed their expression under pluripotent conditions but also elevated the magnitude of differentiation-associated induction of gene expression under differentiation-inducing conditions. These results suggest that Dzip3 represses developmental genes at the lineage-commitment stage in mES cells. Genome-wide ChIP analysis shows relationships between Dzip3-and Ring1B-binding loci.
To determine the genome-wide locations of Dzip3 and Ring1B binding sites and the positions of ubH2A modification in a genome-wide manner, we performed ChIP-seq with anti-Dzip3, anti-Ring1B, and anti-ubH2A antibodies under pluripotent conditions (2i + LIF). A Venn diagram was used to represent the 92 genes that overlapped between the set of genes bound by Dzip3 around their transcription start sites (TSS) and the set of genes upregulated by Dzip3 KD (Fig. 3A). To discover the spectrum of genes that are directly regulated by Dzip3 in mES cells, we analyzed these genes based on their gene ontology (GO) classification, which suggested that Dzip3 regulates developmental processes (Fig. 3B).
Next, we validated the ChIP-seq results by designing primers against the promoter regions of important developmental genes and performing ChIP-qPCR on these genes. We confirmed Dzip3 occupancy by comparing the Dzip3 signal around the promoter of Neurod1, Neurog1, Cdh2 and other genes in mES cells with control (NC) cells and cells subjected to Dzip3 KD (Fig. 3C) Unexpectedly, Ring1B was also enriched at Neurog1 loci. Moreover, Ring1B occupancy levels were decreased not only by Ring1B KD but also by Dzip3 KD. These results suggest that Dzip3 plays a role in localization of Ring1B and regulates gene expression together with Ring1B. Furthermore, ubH2A modification levels were not appreciably altered by Dzip3 KD around the promoter region of Neurog1, although an apparent decline was evident at the Cdh2 and Neurod1 loci. This result suggests that Dzip3 has the potential to repress gene expression not only by ubiquitylating H2A but also by other mechanisms.

Dzip3 regulates developmental genes by reorganizing 3D chromatin conformation.
To explore the molecular mechanisms by which Dzip3 represses gene expression, we focused on those of its target genes around which significant Dzip3 binding sites exist. A recent study revealed that about 25% of Ring1B-bound genes in ESCs possess prominent Ring1B binding sites (RBS) outside of their promoter regions (TSS ± 4 kb) and that Ring1B represses target genes by orchestrating 3D chromatin structural changes 22 . We hypothesized that Dzip3 also participates in regulating target genes by reorganizing 3D chromatin conformation. To assess the validity of this hypothesis, we performed a chromosome conformation capture (3C) assay. We first designed PCR primers for the 3C assay around the promoter region and the 3' end of the gene, as indicated by arrows ( Fig. 4A and Supplementary Fig. 4a). Undigested and unligated samples were prepared as negative controls. For Neurog1, Cdh2, Bspry, and Shank3, we detected signals implicating an interaction between the promoter region and the 3' end of the gene in the presence of the restriction enzyme Mse I and T4 DNA ligase, but not in their absence. Furthermore, the signal indicating chromatin conformation was decreased by Dzip3 KD (Fig. 4A and Supplementary  Fig. 5a). Therefore, we suggest that Dzip3 regulates gene expression by changing the local chromatin conformation.

Discussion
The results of this study indicate that Dzip3 represses differentiation-inducible genes in mES cells, possibly by regulating ubH2A and orchestrating changes in 3D chromatin structure. ubH2A has been reported to inhibit transcription initiation by interfering with H3K4 methylase, inhibiting the elongation of transcription by RNA pol II, or facilitating the recruitment of PRC2 and H3K27 trimethylation 8,11,23 . It was also previously reported that Dzip3 depletion didn't result in significant changes in the general level of ubH2A and that Dzip3 functions by modulating histone modifications at specific sites rather than globally 8 . In our study, KD of Dzip3 was not accompanied by a global decline in ubH2A modification, but instead, changes were restricted to gene promoters. However, our data demonstrated that binding of Dzip3 did not necessarily lead to a decline in ubH2A levels in some gene promoters, and therefore our results suggest that Dzip3 is also able to repress gene expression levels without ubiquitylating H2A. We assume that the ability of Dzip3 to reorganize 3D chromatin conformation, which was uncovered in this study, is closely linked to this H2A ubiquitylation-independent gene repression.
Recent studies on the distribution of PcG proteins in the mammalian genome revealed their preferential association with genes encoding developmental regulators, which exhibit dynamic changes in their spatiotemporal expression during development 24 . Given the results in this study, it is tempting to speculate that Dzip3 substantially contributes to such preferential Ring1B binding at the promoters of genes encoding developmental regulators.
PcG-mediated gene-silencing mechanisms remain elusive; however, it is obvious that there are several distinct mechanisms operating, which include modifications of histone tails (histone H3K27 trimethylation and H2AK119 mono-ubiquitination), condensation of chromosome segments, and mediation of the interactions between topologically distant regulatory elements. In this study, we propose that Dzip3 represses gene expression in conjunction with Ring1B by ubiquitylating ubH2A and inducing conformational changes in 3D chromatin structure (Fig. 4B).
Scientific RepoRts | 5:16567 | DOI: 10.1038/srep16567 RNA interference. mES cells (1.0 × 10 5 ) were seeded in 10-cm plates (previously coated with 0.5% gelatin) with mouse ESC medium. Forty-eight hours after seeding, the cells were transfected using Lipofectamine RNAiMAX reagent (Life Technologies) and Opti-MEM I Reduced Serum medium (Life Technologies), according to the manufacturer's protocol. The targeted siRNA or control siRNA (300 pmol, Life Technologies) were used for transfection. Cells were cultured with a change of medium every 24 hr and collected after 72 hr of siRNA transfection. The siRNA sequence information is shown in Table S1.  Supplementary Fig. 4a. (B) Models for transcriptional repression by Dzip3, which promotes interactions between the promoter and regions distal to its binding site. These changes in 3D chromatin structure repress transcription redundantly with Ring1B. Txn, transcription. Error bars, standard deviation.
Scientific RepoRts | 5:16567 | DOI: 10.1038/srep16567 RNA extraction and RT-qPCR. Total RNA was isolated using ISOGEN II reagent (Nippon Gene), according to the manufacturer's protocol. cDNA was created from 0.5 μ g total RNA using an oligo(dT) primer (Life Technologies), random hexamers (Takara), and M-MuLV reverse transcriptase (NEB). Real-time RT-PCR using reagents containing SYBR green was performed with an ABI PRISM 7900HT instrument (Applied Biosystems). Expression levels were compared with known standard samples and normalized to GAPDH. Primer sequences are shown in Table S2.

Statistical analysis. Statistical analysis was performed using formulae provided in Microsoft Excel.
The frequency of flattened colonies within the population of Dzip3 KD colonies was compared with the negative control (NC) colonies using a hypergeometric distribution analyzed using the HYPGEOMDIST function in Excel. Student's t-test was used to determine the significance level under the assumptions of two separate means with equal variance. The error bars represent the standard deviation (SD) of three independent experiments.

RNA-seq.
Total RNA was used for preparing an RNA-seq library. RNA quality was checked using the Agilent RNA 6000 Nano kit with an Agilent 2100 Bioanalyzer instrument (Agilent Technologies). RNA-seq libraries were prepared using the TruSeq Stranded mRNA LT Sample Prep kit (Illumina) and sequenced with the MiSeq system (Illumina). Samples were sequenced to a depth of approximately 3 million uniquely mapped reads per sample. Sequences were aligned to the mouse MM9 reference genome with the Illumina Analysis Pipeline, allowing one mismatch. Reads that could be uniquely mapped to a gene were used to calculate the expression level. Accordingly, the gene expression level was quantified by the number of uniquely mapped reads per kilobase of exon per million mapped reads (RPKM). To evaluate the correlation between the RNA-seq duplicate sample data sets, scatter plots were created using the Partek ® Genomics Suite. Sequence and gene ontology (GO) analysis were performed with the Partek ® Genomics Suite. Significantly enriched GO functional groups were defined as having an enrichment score equal to or greater than 3 (P value < 0.05), and each functional group was assigned with a GO enrichment score calculated using Fisher`s exact test. All RNA-seq data can be found online in the NCBI GEO SuperSeries GSE71884.

Heat map construction and hierarchical clustering. Heat map construction and hierarchical clus-
tering of the gene expression profiles were performed using the Partek ® Genomics Suite. Heat maps for each sample (negative control [NC], Dzip3 KD and Ring1B KD) were normalized by calculating the RPKM based on the sum of all reads found in the exon regions of that gene. One-way analysis of variance (ANOVA) was used to identify differentially expressed genes. By setting P < 0.1 and fold-change (FC) settings FC > 1.5 or FC > 2, we obtained lists of differentially expressed genes between NC and Dzip3 KD or NC and Ring1B KD cells. Hierarchical clustering of NC, Dzip3 KD, and Ring1B KD cells was performed using 894 genes that varied significantly among each groups with a statistical P value < 0.1. We subtracted the mean of the gene expression levels in the six paired samples, normalized each row in the data table, calculated the distance using Pearson correlation, and then used a "pairwise average-linkage" hierarchical clustering method for clustering. The scale was the standardized RPKM value.

Nuclear extraction and immunoprecipitation. For nuclear extraction, cells were washed with PBS
and incubated on ice for 30 min in buffer (10 mM Tris-HCl, pH 7.9, 10 mM KCl, 1.5 mM MgCl2, 1 mM DTT, 1 mM PMSF, 1 mM sodium metabisulfite). Next, cells were homogenized and centrifuged to pellet the nuclei, which were then suspended with buffer (25 mM HEPES-KOH, 0.02 mM EDTA, 10% glycerol, 0.01% NP40, 0.2 M KCl), sonicated (12% of maximum power, 10 sec, 3 cycles), and the supernatant collected as a nuclear extract. For immunoprecipitation, the supernatant was incubated with anti-mRing1B or anti-mDzip3 antibodies or purified rabbit IgG as negative control at 4 °C overnight. To precipitate the immune complexes, Protein A Sepharose 4 Fast Flow (GE Healthcare) was added, and the vessel rotated at 4 °C for 1 hr. The beads were washed, and the immune complexes were eluted in 50 mM glycine.
ChIP-seq. ChIP samples were prepared with a modification of our ChIP-qPCR method. Briefly, mES cells cultured in 2i + LIF were fixed with 1% formaldehyde and processed by fragmentation with micrococcal nuclease (Sigma-Aldrich). For fragmentation, fixed samples were washed with lysis buffer (50 mM HEPES-KOH, 140 mM NaCl, 1 mM EDTA, 10% glycerol, 0.5% NP-40, 0.25% Triton X-100), resuspended in micrococcal nuclease buffer (15 mM Tris-Hcl, 5 mM MgCl 2 , 1 mM CaCl 2 , 25 mM NaCl), and incubated with micrococcal nuclease at 37 °C for 20 min. For chromatin extraction, samples were sonicated (20% of maximum power, 12 sec.) in lysis buffer #2 (10 mM Tris-HCl, 100 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 0.1% Na-deoxycholate, 0.5% N-lauroylsarcosine) and rotated in 1% Triton X-100 to quench the sarkosyl. The supernatant was collected as chromatin. Input DNA fragmentation was confirmed by electrophoresis, and the amount of nucleotide was measured with a Nanodrop ND-1000 spectrophotometer (Thermo Scientific). The same amount of input DNA was mixed with antibody (anti-mRing1B, anti-mDzip3, or anti-ubH2A antibodies or purified rabbit IgG as negative control) and sepharose A beads to react and rotated at 4 °C overnight. The beads were washed with RIPA buffer (50 mM HEPES-KOH, 500 mM LiCl, 1 mM EDTA, 1.0% NP-40, 0.7% Na-deoxycholate) and incubated at 37 °C for 30 min in TE with added RNase (50 μ g/ml final concentration). For reverse cross-linking and elution, the beads were incubated at 65 °C for 10 hr in elution buffer (1% SDS, 50 mM Tris-HCl, 10 mM EDTA). Elution samples were incubated at 55 °C for 2 hr after adding proteinase K (200 μ g/ml final concentration). DNA from the eluates was purified by phenol-chloroform extraction and ethanol precipitation. ChIP-seq libraries were prepared from ChIP samples using the ChIP-Seq Sample Prep kit (Illumina). The resulting libraries were sequenced with the MiSeq sequencing system (Illumina). Samples were sequenced to a depth of approximately 30 million uniquely mapped reads per sample. Sequences were aligned to the mouse MM9 reference genome with the Illumina Analysis Pipeline, allowing one mismatch. To compensate for differences in sequencing depth and mapping efficiency, the data were normalized using MACS2 software, resulting in signal per million reads (SPMR). Peak calling was also performed using the Partek ® Genomics Suite ™ (PGS), with a statistical false discovery rate (FDR) value of 0.1. The genes having Dzip3 peaks within 100 kb upstream or 100 kb downstream of their TSS were determined using the PGS command "Find nearest genomic features". GO classifications were also assigned using PGS. To visualize the ChIP-seq results, we used the Integrative Genomics Browser (BioViz). All ChIP-seq data can be found online in the NCBI GEO SuperSeries GSE71884. Two-step cross-linking. For chromatin preparation, cells were washed on a plate three times with PBS and incubated with 2 mM disuccinimidyl glutarate (DSG, Thermo Scientific) in PBS for 30 min at room temperature. ES cells were washed in PBS three times, fixed for 10 min in 1% formaldehyde, neutralized with glycine (0.124 M), and stored in -80°C until use. All subsequent ChIP-qPCR steps were performed as described above.
3C-qPCR. 3C assays were performed as described previously 26 . Briefly, mES cells cultured in 2i + LIF were fixed with 1% formaldehyde and resuspended in lysis buffer. Nuclei were treated with Mse I restriction enzyme (NEB) and reacted with T4 DNA ligase (NEB). The genomic DNA was then reverse Scientific RepoRts | 5:16567 | DOI: 10.1038/srep16567 cross-linked and purified by phenol-chloroform extraction and ethanol precipitation. Each purified DNA concentration was quantified using a Nanodrop ND-1000 spectrophotometer to adjust the qPCR concentration to 50 ng/μ l. Quantitative PCR was performed using the KAPA SYBR FAST qPCR kit with the ABI PRISM 7900HT instrument (AQ method). To generate standard curves, duplicate samples with 30-, 150-, and 750-fold diluted DNA were used. Relative values for 3C were normalized using qPCR values of a GAPDH allele lacking restriction sites within its PCR product. Primer sequences are shown in Table S2.