ABSTRACT
The genome is transcriptionally inert at fertilization and must be activated through a remarkable developmental process called zygotic genome activation (ZGA). Epigenetic reprogramming contributes significantly to the dynamic gene expression during ZGA; however, the mechanism has yet to be resolved. Here, we find histone deacetylases 1 and 2 (HDAC1/2) can regulate ZGA through lysine deacetylase activity. Notably, in mouse embryos, overexpression of a HDAC1/2 dominant-negative mutant leads to developmental arrest at the two-cell stage. RNA-seq reveals that 64% of downregulated genes are ZGA genes and 49% of upregulated genes are developmental genes. Inhibition of the deacetylase activity of HDAC1/2 causes a failure of histone deacetylation at multiple sites, including H4K5, H4K16, H3K14, H3K18 and H3K27. ChIP-seq analysis exhibits an increase and decrease of H3K27ac enrichment at promoters of up- and downregulated genes, respectively. Moreover, HDAC1 mutants prohibit the removal of H3K4me3 by impeding expression of Kdm5 genes. Importantly, the developmental block can be greatly rescued by Kdm5b injection and by partially correcting the expression of the majority of dysregulated genes. Similar functional significance of HDAC1/2 is conserved in bovine embryos. Overall, we propose that HDAC1/2 are indispensable for ZGA by creating correct transcriptional repressive and active states in mouse and bovine embryos.
INTRODUCTION
The highly differentiated germ cells (sperm and egg) go through the fertilization process to form a totipotent zygote, which marks the beginning of a new life. However, the new life must complete a developmental process called maternal zygotic transition (MZT) or oocyte-to-embryo transition (OET) to acquire its developmental independence (Schultz et al., 2018; Vastenhouw et al., 2019). During this transition, the genome is transcriptional quiescent initially and the development control relies on maternal derived transcripts and proteins. These maternal factors gradually undergo clearance while the zygotic genome is activated (zygotic genome activation; ZGA). ZGA typically occurs in two consecutive waves, termed minor and major ZGA. The timing of ZGA varies among mammals. In mice, minor and major ZGA occurs in late one-cell and mid-to-late two-cell embryos, respectively. The gene expression pattern established during ZGA is required for setting up totipotency in early embryos and the ensuing development. Therefore, a central question in biology is the defining feature of chromatin states during ZGA and how it contributes to changes in gene expression.
A considerable epigenetic reprogramming is a remarkable feature of early stage embryogenesis (Eckersley-Maslin et al., 2018; Du et al., 2021; Xu et al., 2021). In particular, DNA methylations and various histone modifications are subject to dramatic remodeling, which is believed to ensure the correct gene expression program and ZGA. Indeed, aberrant DNA methylation and histone modifications in fertilized or cloned embryos have been closely linked to abnormal gene expression and developmental failure (Matoba et al., 2014; Gao et al., 2018). However, studies of molecular mechanisms of epigenetic reprogramming are traditionally limited by the scarce samples and a lack of tools for chromatin analysis in early embryos.
With the rapid advance of low-input chromatin profiling methods including micro-scale ChIP-seq and CUT&RUN (Brind'Amour et al., 2015; Dahl et al., 2016; Liu et al., 2016; Zhang et al., 2016; Skene et al., 2018; Kaya-Okur et al., 2020), how histone modifications patterns are passed, reprogrammed and established from oogenesis to early embryogenesis have been unveiled, and their correlation with gene expression profiles determined (Dahl et al., 2016; Liu et al., 2016; Zhang et al., 2016; Xu et al., 2021). However, the causal relationship between these histone modifications and gene expression, as well as the crosstalk between histone modifications have yet to be resolved.
Histone acetylation is widely considered a hallmark of active gene expression (Dahl et al., 2016). It has been shown that H3K9ac regulates the transition from transcription initiation to elongation (Gates et al., 2017). Importantly, insufficient H3K9ac acts as an epigenetic barrier in mouse somatic cell nuclear transfer (SCNT) embryos (Yang et al., 2021). During oocyte maturation, an aberrant increase of H4K16ac induces defective chromatin structure and impairs oocyte growth and development (Ma and Schultz, 2013; Zhang et al., 2014). As the most-characterized histone acetylation, H3K27ac is associated with increased release of RNA polymerase II into active transcription, and is typically accumulated at promoters and enhancers of active expressed genes (Creyghton et al., 2010). In zebrafish embryos, there is a global loss of H3K27ac marks at enhancers, while widespread de novo H3K27ac is acquired at promoters prior to ZGA (Zhang et al., 2018). Therefore, these results indicate the importance of histone acetylation reprogramming in early embryonic development.
Histone acetylation is directly controlled by histone acetyltransferases (HATs) and histone deacetylases (HDACs) (Sheikh and Akhtar, 2019). Mutations of these enzymes often lead to embryonic lethality in mice and are connected with human diseases (Sheikh and Akhtar, 2019). HDAC1 is the first HDAC to be characterized. To date, four classes of 18 HDACs have been identified in mammals. Among them, class I family members (HDAC1, HDAC2, HDAC3 and HDAC8) display the highest sequence homology to HDAC1 and are ubiquitously expressed, localized in the nucleus and exhibit high catalytic activity towards histone substrates. HDAC1 and HDAC2 (HDAC1/2) are the most highly expressed HDACs in mouse oocytes and embryos (Dahl et al., 2016). Genetic evidence reveals an overlapping and essential role for HDAC1/2 during oogenesis (Ma et al., 2012). Previous research demonstrates an important role for Hdac1 in preimplantation development in mice using a knockdown approach (Ma and Schultz, 2008). Moreover, our recent work established a redundant role of HDAC1/2 in mouse early embryogenesis (Zhao et al., 2020). However, the specific role of HDAC1/2 in ZGA remains unclear, given the technical difficulty in acutely removing maternal HDAC1/2.
In this study, we demonstrate that the lysine deacetylase activity of HDAC1/2 is required for ZGA in both mice and cattle by using a combination of a dominant-negative approach and a pharmacological approach. HDAC1/2 mutants lead to a large-scale dysregulation of gene expression during ZGA and a developmental arrest at the two-cell stage. Mechanistically, HDAC1 mutation results in a deviant genome-wide distribution of histone acetylation and prohibits the removal of the broad H3K4me3 domain by disrupting the expression of Kdm5 genes. Moreover, the developmental phenotype can be greatly rescued through injection of Kdm5b mRNA. We propose that HDAC1/2 are critical epigenetic modifiers that maintain the transcriptional states during ZGA through their lysine deacetylase activity.
RESULTS
HDAC1/2 are essential for ZGA in both mouse and bovine embryos
To study the functional role of HDACs in ZGA, we first validated the expression pattern of HDACs in mouse and bovine embryos. Transcriptomic analysis revealed that Hdac1/HDAC1 and Hdac2/HDAC2 were the most abundant among all Hdacs/HDACs in both mouse (Fig. S1A) and bovine early embryos (Fig. S1B). In particular, Hdac1 mRNA abundance increased gradually and visibly during ZGA in both species (Fig. S1A,B). In accordance with the mRNA change, the nuclear intensity of HDAC1 was increased during ZGA (Fig. S1C,D). These results suggest that HDAC1/2 might play a role during ZGA.
To delve into the function of HDAC1/2 in ZGA, we treated the zygotes with FK228, a HDAC1/2-specific inhibitor (Furumai et al., 2002). Embryo culture results showed FK228-treated embryos arrested at the two-cell stage in mice (Fig. S1E). Concomitantly, FK228 treatment resulted in a developmental arrest in cattle, and cell number analysis with DAPI staining clarified that these embryos arrested at the 8- or 16-cell stage (bovine ZGA; Fig. S1F,G).
To probe the specific effects of the lysine deacetylase activity of HDAC1, we generated a mutant in the lysine deacetylase domain. The lysine deacetylase domain, which is conserved in all Class I HDACs, consists of a stretch of more than 300 amino acids, and a mutation in the amino acid sequence destroys its enzymatic activity (Hassig et al., 1998). We therefore constructed a HDAC1 mutant in which the residue histidine 141 (H141) was replaced by alanine; the mutant can abolish deacetylase activity of HDAC1 by competitive inhibition without destroying the interaction with HDAC1-associated proteins (Hassig et al., 1998). The resultant mutant is referred to here as H1MU, while the wild type is referred to as H1WT. We microinjected early zygotes with mRNA for H1WT and H1MU. H2O was injected in control groups (Fig. 1A). Nearly all the embryos in H1MU group suffered developmental arrest at the two-cell stage (87%) while the control embryos exhibit robust development, with 80% developing to blastocyst stage (Fig. 1A-C).
HDAC1 and HDAC2 are two homologous histone deacetylases that exhibit functional redundancy in oogenesis and preimplantation development (Ma et al., 2012; Zhao et al., 2020). We therefore asked whether these two proteins play a redundant role during ZGA. The Hdac2 mutant was also generated by employing a mutation in the catalytic domain (Kurita et al., 2012). Likewise, the mutation also caused embryonic arrest at the two-cell stage (Fig. 1D). Overall, these data revealed that HDAC1/2 are crucial for ZGA in both mouse and bovine embryos.
HDAC1 mutation results in an aberrant gene expression pattern during ZGA
To gain insight into developmental failure of H1MU embryos at the molecular level, we measured the level of phospho-Ser2 (Ser2P), a hallmark of RNA polymerase II activity, to examine whether the inhibition of HDAC1 affect the genome-wide transcriptional program. However, there was no visible difference between H1MU and H1WT embryos (Fig. S2).
Furthermore, we collected early two-cell embryos from H1WT and H1MU groups, and late two-cell stage embryos from H2O, H1WT and H1MU groups, and performed RNA-seq (Fig. 2A). Two independent replicates of RNA-seq samples from each group displayed high correlation (Fig. S3A). In addition, Hdac1 was significantly increased in both H1WT and H1MU groups, further confirming a robust overexpression efficiency (Table S2).
H1MU has no significant impact on minor ZGA, as the analysis revealed only 77 and 60 upregulated and downregulated genes, respectively (Padj≤0.05; fold change ≥2 or ≤0.5), at early two-cell stage in H1MU versus H1WT groups (Fig. 2B). Moreover, at late two-cell stage, there are only 88 genes differentially expressed in wild-type groups compared with control groups, suggesting H1WT overexpression has a minor effect on the transcriptome during ZGA (Fig. S3B). In contrast, 6565 genes were differentially expressed in late two-cell embryos of H1MU compared with H1WT groups with 3898 and 2667 upregulated and downregulated genes, respectively (Fig. 2C and Table S2). These transcriptomic comparisons at the two different embryonic stages clearly indicated that HDAC1 specifically regulates gene transcription when major ZGA occurs.
Gene ontology (GO) enrichment analysis of upregulated genes revealed overrepresentation of genes involved in transcription and sequence-specific DNA binding (Table S3). Meanwhile, similar analysis of downregulated genes showed that a significant number of genes were involved in transcription and RNA splicing, ribosome biogenesis, and nucleotide metabolism, which are known markers of ZGA (Zeng and Schultz, 2005) (Table S4). We also noticed a drastic reduction in the expression of Smarca4 (also known as Brg1) (Bultman et al., 2006) and Nfya (Lu et al., 2016) (Table S2), two genes known for their putative function in regulating mouse ZGA.
To assess the molecular features of target genes of HDAC1, we then examined the expression pattern of the upregulated genes and downregulated genes during normal preimplantation development. We categorized the upregulated and downregulated genes into three main classes each using k-means clustering. Surprisingly, a considerable number of upregulated genes are both maternal (Cluster 1) and developmental genes (Cluster 3), which are supposed to be silent, whereas the majority of downregulated genes are ZGA genes (Cluster 5; Fig. 2D-E). These results suggest the enzymatic activity of HDAC1 is required not only for transcriptional silencing but for transcriptional activation during ZGA.
To address directly whether there is functional redundancy of HDAC1/2 in regulating gene transcription, we compared differentially expressed genes (DEGs) caused by H1MU and H2MU. Results indicated that the majority of upregulated and downregulated genes are common to both groups (Fig. 2F and Fig. S3C), suggesting a compensatory role for HDAC1/2 in transcriptional regulation. We also observed that FK228 treatment resulted in similar changes in the transcriptome (Fig. 2G, Fig. S3D,E), indicating its specificity in inhibiting the lysine deacetylase activity of HDAC1/2. In agreement with the result of H1MU, we observed that the majority of upregulated and downregulated genes in H2MU or FK228-treated groups are maternal/developmental genes and ZGA genes, respectively (Fig. S3F-H). Furthermore, clustering analysis showed that the transcriptome of embryos in the H1MU, H2MU and FK228 group is distinguishable from those of wild-type early two-cell embryos, which demonstrated that the observed transcriptional changes are not due to developmental delay (Fig. S3E). Overall, these data confirmed the importance of HDAC1/2 for ZGA.
The lysine deacetylase activity of HDAC1/2 is required for reprogramming of histone acetylation during ZGA
Given the deacetylase activity of HDAC1/2, we sought to determine whether HDAC1/2 mutation affects the reprogramming of histone acetylation. Immunofluorescence results showed a decrease of H3K27ac, H3K14ac, H3K18ac, H4K5ac and H4K16ac in H1WT embryos compared with those in H2O-injected controls (Fig. 3A and Fig. S4A-D), confirming the effective overexpression of HDAC1. Importantly, we observed a robust increase of H3K27/K14/K18ac and H4K5/K16ac in H1MU embryos (Fig. 3A-B and Fig. S4A-D). We validated that similar changes in H3K27ac were found in FK228-treated group and H2MU groups (Fig. 3C-D and Fig. S4E-F). Hence, the developmental failure observed upon inhibition of HDAC1/2 might be ascribed to the aberrant histone acetylation reprogramming that occurs during ZGA.
The defective ZGA accompanied by hyperacetylation described above prompted us to explore the correlation between histone acetylation and ZGA. For this purpose, we measured the reprogramming of H3K27ac, which is a well-established marker of active chromatin in somatic cells. In mice, relatively strong signals of H3K27ac were detected in both maternal and paternal pronuclei at zygote stage and its nuclear intensity drastically diminished from the early and late two-cell stages, concurrent with mouse major ZGA (Fig. 3E). In cattle, H3K27ac is still relatively bright at the four-cell stage, when minor ZGA occurs; however, it underwent a sharp decrease at the 8- and 16-cell stage, coinciding with bovine major ZGA (Fig. 3F). These data are in agreement with previous studies in mice (Santenard et al., 2010) and such coincidence with ZGA has also been reported in pigs (Zhou et al., 2014), suggesting that this global loss of histone acetylation is conserved in mammals.
HDAC1 mutation results in aberrant distribution of histone acetylation in early embryos
To explore the correlation between hyperacetylation and dysregulated gene expression upon the inhibition of HDAC1, we performed Ultra-Low-Input-NChIP-seq (ULI-NChIP-seq) of H3K27ac at late two-cell stage and obtained high-quality reproducible data (Fig. 4A,B, and Fig. S5A). The protocol used here is highly efficient as our data in control groups faithfully recapitulates the previously published H3K27ac ChIP-seq data (Dahl et al., 2016) with a high correlation (R=0.80; Fig. S5B,C). Genome distribution of H3K27ac changed significantly in H1MU embryos. In particular, H3K27ac enrichment at promoter and distal regions were decreased and increased, respectively (Fig. 4C). By scanning the genome with a 5 kb sliding window, we identified H3K27ac-lost and -gained regions in H1MU embryos. In consistent with our immunofluorescence results, H3K27ac-gained regions covered more genome than H3K27ac-lost regions (Fig. 4D).
To establish whether the epigenetic changes associated with the transcriptomic differences, we profiled H3K27ac signal in promoters of upregulated and downregulated genes. H3K27ac accumulated in both promoters and gene bodies of upregulated genes in the H1MU group (Fig. 4E and Fig. S5D). On the contrary, the promoters and gene bodies of downregulated genes displayed reduced H3K27ac signal (Fig. 4F and Fig. S5E). These results suggest that the disorder of gene expression pattern caused by H1MU is correlated with the aberrant histone deacetylation.
DUX is upregulated upon HDAC1 mutation
To identify sequence motifs for DNA-binding proteins that were potentially enriched in aberrant H3K27ac-deposited regions, we scanned the entire genome within a 1 kb window to calculate H3K27ac signals, and annotated motifs at the top 2% H3K27ac-increased and -decreased regions by HOMER. Notably, DUX was significantly enriched in the H3K27ac-increased regions (Fig. 5A). However, these motifs were not detected in H3K27ac-decreased regions. Altogether, these results suggest that H1MU influences expression of active and repressed genes through distinct regulatory mechanisms.
DUX (DUX4 in humans) is a putative transcriptional factor involved in regulation of mammalian ZGA (De Iaco et al., 2017, 2020; Hendrickson et al., 2017). Dux is supposed to be transiently expressed during early to mid two-cell stage and rapidly disappears at late two-cell stage (Fig. 5B) (Guo et al., 2019). However, although there was a distinct decrease of Dux in late two-cell embryos of H1MU groups, the amount was still greater than that in control embryos (Fig. 5B). Zscan4s and ERVLs are identified as target genes of DUX, and can be robustly expressed in mouse embryos overexpressing DUX (Guo et al., 2019). We found the expression of Zscan4 genes and ERVLs (MERVL-int and MT2_mm) was also increased in H1MU embryos (Fig. 5C-E). Therefore, it is highly possible that the insufficient repression of Dux is one leading factor causing the upregulation of repressed genes, including developmental and maternal genes.
It has been reported that DUX4 function as a DNA-binding protein that recruits EP300/CBP through its C-terminals to activate expression of nearby genes (Choi et al., 2016). Intriguingly, the average mRNA abundance of upregulated genes in late two-cell stage H1WT embryos is lower than those of downregulated genes (Fig. 5F), suggesting the timely loss of DUX in late two-cell embryos is crucial to maintain the low expression of those upregulated genes.
HDAC1 mutation inhibits the removal of the broad H3K4me3 domain during mouse ZGA
H3K4me3 is an established hallmark of permissive promoters that can tether chromatin remodelers and transcriptional regulators. Genome-wide distribution of H3K4me3 is unique in early embryos for its broad domains established during oogenesis and removed upon ZGA (Dahl et al., 2016; Liu et al., 2016; Zhang et al., 2016). Next, we wondered whether H3K4me3 reprogramming was affected in H1MU embryos. We observed a notable level of H3K4me3 at early two-cell stage in control groups and no difference was found in H1MU groups (Fig. S6A). Furthermore, H3K4me3 signal can be barely seen in control embryos, while the intensity was obviously increased in H1MU embryos at late two-cell stage, suggesting the removal of broad H3K4me3 domain is blocked (Fig. 6A,B). Similar to the effect of H1MU, we also found that the erasure of the broad H3K4me3 domain was suppressed when H2MU or FK228 treatment was employed (Fig. S6B-D).
We then analyzed the global dynamics of H3K4me3 in H1MU and control embryos by ChIP-seq. We validated the H3K4me3 ChIP first and found there was high correlation between our data in control groups and the previously published H3K4me3 ChIP-seq data (Zhang et al., 2016) (R=0.81; Fig. S7A,B). The two biological replicates for both groups exhibit high correlation, indicating a high reproducibility (Fig. S7C,D). We identified 34,590 H3K4me3-gained domains and 5890 H3K4me3-lost domains in H1MU embryos, which covered 4.5% and 0.7% of the genome, respectively. Moreover, H3K4me3-gain regions accounted for nearly 50% at distal regions, whereas H3K4me3-lost regions accounted for 25% at distal regions, further suggesting a block in removal of the broad H3K4me3 domain (Fig. 6C).
H3K4me3 reprograming is crucial to ensure the gene expression program during ZGA. We therefore sought to determine whether H3K4me3 involved in transcriptomic changes in H1MU embryos. Surprisingly, H3K4me3 was reduced at promoters of both up- and downregulated genes in the H1MU group (Fig. 6D-G). However, H3K4me3 signals were enhanced in the gene body regions of upregulated genes whereas they were reduced in gene body regions of downregulated genes (Fig. 6D-G), suggesting that the dysregulated gene expression could be ascribed to the sustained broad H3K4me3 domain. Altogether, the enzymatic activity of HDAC1 is crucial for the immediate removal of broad H3K4me3 domain upon ZGA.
Broad H3K4me3 domains are removed by KDM5 family members during mouse ZGA (Dahl et al., 2016; Liu et al., 2016; Zhang et al., 2016). Interestingly, GO analysis of downregulated genes revealed genes involved in histone H3-K4 demethylation (Table S4). Indeed, mRNA abundance of all three Kdm5 genes (Kdm5a, Kdm5b and Kdm5c) was dramatically decreased in H1MU embryos (Fig. 6H). Kdm5b is zygotic transcribed from early to late two-cell stage (Fig. S7E) whereas H3K27ac is accumulated at Kdm5b (Fig. 6I). In contrast, H1MU caused a reduction in H3K27ac enrichment at Kdm5b (Fig. 6J), suggesting H3K27ac regulates transcription of Kdm5b.
KDM5B is a key mediator of the HDAC1 MU phenotype
We next asked whether ectopic injection of Kdm5b mRNA could rescue the H3K4me3 pattern and the developmental outcome of H1MU embryos (Fig. 7A). Because Kdm5b normally begins its transcription at the two-cell stage, we thus performed microinjection of Kdm5b mRNA into two individual blastomeres of early two-cell stage embryos. Results showed that microinjection of Kdm5b resulted in a significant decrease in global H3K4me3 in H1MU embryos, indicating a recovery of H3K4me3 reprogramming (Fig. 7B). As a consequence, the majority of the co-injected embryos could overcome the two-cell block (Fig. 7C,D), with 46% developing to the four or eight-cell stage and 14% even developing to morula stage (Fig. S7F). RNA-seq analysis revealed that the expression of 53.99% of downregulated genes and 42.21% of upregulated genes in H1MU groups could be at least partially corrected by ectopic expression of Kdm5b (Fig. 7E), and the rescued genes were significantly enriched in major ZGA (Fig. 7F), suggesting that the accurate distribution of H3K4me3 is required for ZGA. Remarkably, expression of Dux dropped to normal levels when KDM5B was overexpressed, and expression of its target Zscan4 genes was also reduced greatly (Fig. 7G-H), suggesting that HDAC1 regulates DUX through KDM5B. Together, these findings implicate KDM5B as a key molecular mediator of HDAC1 function in controlling H3K4me3 reprogramming.
The function of HDAC1 on ZGA is conserved in mouse and bovine embryos
Recent studies in epigenome reprogramming have revealed remarkable differences among mammals (Hanna et al., 2018; Xia et al., 2019). To determine whether the effect of HDAC1 on transcriptional regulation during ZGA is conserved, we collected DMSO and FK228-treated eight- to 16-cell embryos and performed RNA-seq. We identified 1990 and 2794 upregulated and downregulated genes, respectively (Fig. 8A,B and Fig. S8A,B). Seventy-four percent (2061/2794) of the downregulated genes are supposed to be expressed at the eight- to 16-cell stage (Fig. 8B-C and Fig. S8B-C), indicating a ZGA defect. Strikingly, KDM5 genes were repressed in FK228-treated embryos (Fig. 8D), which was accompanied by an increase in H3K4me3 during ZGA (Fig. 8E,F). Furthermore, we microinjected H1MU mRNA into bovine zygotes (Fig. S8D) and found similar results to the FK228 experiment. Most of the H1MU embryos were arrested at the eight- to 16-cell stage (Fig. S8D,E), with an obvious increase in H3K27ac and H3K4me3 compared with control embryos (Fig. S8F,G). In summary, these results suggest that HDAC1/2 play a conserved role in the regulation of gene expression patterning during ZGA.
DISCUSSION
A fundamental question in developmental biology is what factors trigger ZGA in mammals. Accumulative evidence has revealed a growing list of proteins involved in the process (Bultman et al., 2006; Lu et al., 2016; De Iaco et al., 2017; Hendrickson et al., 2017). However, the molecular mechanisms underlying ZGA have yet to be determined (Mal et al., 2001; Humphrey et al., 2008). Here, we demonstrate that the lysine deacetylase activity of HDAC1/2 are crucial for establishing correct gene expression profile during ZGA. Mechanistically, HDAC1/2 regulate gene expression likely through controlling the reprogramming of histone acetylation and H3K4me3. Importantly, HDAC1 acts as an upstream factor of KDM5B to regulate expression of target genes. We propose that HDAC1/2 are not only involved in the development of the transcriptional repressive state for repressed genes but also in the development of the transcriptional active state for ZGA genes. To our knowledge, HDAC1/2 represent the first proteins that hold a dual role in transcriptional regulation during ZGA.
We demonstrated that the lysine deacetylase activity of HDAC1/2 is responsible for development progression throughout ZGA. In agreement with our data, previous reports show that HDAC1 is zygotically expressed from the two-cell stage and only HDAC1 is sensitive to α-amanitin among the HDACs expressed in two-cell embryos (Zeng and Schultz, 2005; Dahl et al., 2016). Upon fertilization, HDAC1 seems to be one of key molecular players directly involved in the development of the transcriptional repressive state for maternal and developmental genes. Indeed, inhibition of HDAC1 leads to upregulated expression of 3898 genes and the majority of them belong to maternal and developmental genes.
Interestingly, we further identified that Dux was aberrantly upregulated upon HDAC1 mutation. As an intron-less and multi-copy gene that encodes a double-homeobox transcriptional factor involved in major ZGA, Dux is transiently transcribed during minor ZGA and its timely removal is crucial for early embryogenesis (Guo et al., 2019), of which the biological significance remains unclear. It is possible that the DUX removal is crucial to inhibiting the transcription activity for those upregulated genes induced by HDAC1 mutation. Indeed, human DUX4 can recruit EP300 through its C-terminal domain and therefore is involved in histone acetylation (Choi et al., 2016; Bosnakovski et al., 2019). However, further investigation is needed to explain how HDAC1/2 regulate Dux expression.
Histone acetylation has long been linked with the transcriptional permissive state. We observed a sharp decrease of global H3K27ac in both mouse and cattle. This phenomenon also occurs at other histone acetylation sites during ZGA in mice, cattle and pigs (Santenard et al., 2010; Zhou et al., 2014), suggesting that histone deacetylation is a conserved event that establishes the chromatin state before ZGA. Hence, we speculate that histone deacetylation is one upstream inducer of the major ZGA as it occurs before the major ZGA and earlier than H3K4me3 reprogramming (early to late two-cell stage), which appears as a ZGA regulator (Dahl et al., 2016). Indeed, mutation of HDAC1/2 elicits a dramatic downregulation of ZGA genes. These genes account for 24% of all ZGA genes. Interestingly, H3K27ac was decreased at promoters of these genes after H1MU injection, in contrast to the deacetylase activity of HDAC1. Thus, HDAC1 may regulate transcriptional activity of these active genes in an indirect manner, which awaits additional investigation.
Crosstalk of different epigenetic modifications makes up a precise regulatory network for gene expression (Eckersley-Maslin et al., 2018; Xu and Xie, 2018). Kdm5 genes were downregulated in embryos injected with mutant Hdac1 mRNA and caused inadequate removal of broad H3K4me3 in gene body and intergenic regions. Timely removal of broad H3K4me3 is viewed to be associated with mouse ZGA, but Kdm5 gene-knockdown embryos can develop beyond ZGA, although a set of ZGA genes was slightly downregulated at the two-cell stage (Dahl et al., 2016; Liu et al., 2016), suggesting that abnormal H3K4me3 alone plays limited roles in ZGA. However, the inhibition of HDAC1 could cause insufficient removal of histone acetylation and H3K4me3 concurrently, and led to developmental arrest at the two-cell stage. Surprisingly, the injection of exogenous Kdm5b mRNA could partly rescue the development of HDAC1 mutant embryos, suggesting that the role of HDAC1/2 in ZGA is at least partly mediated through regulation of H3K4me3 reprogramming.
Previous studies have shown substrates of HDAC1/2 include not only histones but also non-histones (Luo et al., 2000; Nalawansha et al., 2017). Here, we found that the activity of HDAC1/2 is required for deacetylation at specific lysine sites of histones. In particular, we determined the changes in H3K27ac given its well-established association with active gene expression. However, we cannot rule out the possibility that the inhibition of HDAC1/2 could potentially affect other substrates, including non-histones, which warrants further investigations in the future.
Altogether, we demonstrate an indispensable role for the lysine deacetylase activity of HDAC1/2 for ZGA in both mice and cattle. HDAC1/2 are required for the proper reprogramming of histone acetylation and H3K4me3. The chromatin reprograming is crucial for the development of both a transcriptional repressive state for silent genes and a transcriptional active state for ZGA genes.
MATERIALS AND METHODS
Ethics statement
All procedures involving laboratory animals were carried out in accordance with the guidelines for the care and use of laboratory animals, and approved by Zhejiang University.
Mouse embryo collection and in vitro culture
All experimental mice were raised and housed in a temperature-controlled room (22-25°C) with a relative humidity of 60% in a 12 h light/dark cycle in the Laboratory Animal Center at Zhejiang University. All mice were provided with access to food and water ad libitum. Eight to 10-week-old female BDF1 (C57BL/6×DBA/2; Beijing Vital River Laboratory Animal Technology) mice were super-ovulated by an intraperitoneal injection of 7.5-10 IU PMSG (Sansheng Pharmaceutical, Ningbo, China). Two days later, 7.5-10 IU hCG (Sansheng Pharmaceutical) were administrated and mated with BDF1 male mice. 17-19 h after hCG injection, mice were sacrificed and zygotes were collected from the swollen upper part of the oviduct. Zygotes were incubated in the hyaluronidase (Sigma-Aldrich) solution to remove the cumulus cells and subsequently cultured in KSOM (Millipore) micro-drops under mineral oil at 37°C in 5% CO2.
In vitro production of bovine embryos
Procedures of bovine in vitro production, including in vitro maturation (IVM), in vitro fertilization (IVF) and in vitro culture (IVC), were carried out routinely in accordance with previously published studies (Wang et al., 2020; Shi et al., 2021). Briefly, cumulus-oocyte complexes (COCs) with more than three layers of cumulus cells were collected from ovaries derived from a local slaughterhouse. IVM was conducted with Medium-199 (M4530) supplemented with 10% FBS (Gibco), 1 IU/ml FSH (Sansheng Biological Technology), 0.1 IU/ml LH (Solarbio, Beijing, China), 1 mM Na Pyruvate (Thermo Fisher Scientific), 2.5 mM GlutaMAX (Thermo Fisher Scientific) and 10 μg/ml gentamicin. The IVM condition was 38.5°C under 5% CO2 in humidified air for 22-24 h. IVF was performed by incubating COCs (60-100 COCs per well in four-well plates) with spermatozoa (1∼5×106), which were purified from frozen-thawed semen by using a Percoll gradient in BO-IVF medium (IVF Bioscience). IVF conditions were 38.5°C under 5% CO2 for 9-12 h. Cumulus cells were then removed from putative zygotes by pipetting up and down in Medium-199 (M7528) supplemented with 2% FBS (Gibco-BRL). Embryos were cultured in BO-IVC medium (IVF Bioscience) at 38.5°C under 5% CO2 in humidified air until use.
FK228 treatment experiments
Romidepsin (FK228, Depsipeptide) dissolved in dimethyl sulfoxide (DMSO, Sigma-Aldrich) was added to KSOM (for mice embryo) or BO-IVC medium (for bovine embryo) at a final concentration of 50 nM. In the control group, an equivalent amount of DMSO (0.1%) was added.
In vitro transcription and microinjection
Wild-type or mutant cDNA for Hdac1, Hdac2 and Kdm5b were subcloned into T7-driven vectors. The Hdac1 mutant (H141A) and Hdac2 mutant (H141A) were constructed as described previously (Ito et al., 2002; Kurita et al., 2012). All sequences were validated by Sanger sequencing before use (Sangon, Shanghai, China). To prepare mRNAs for microinjection, the plasmids were linearized and in vitro transcribed, capped, and poly(A) tailed using T7 mMESSAGE mMACHINE Ultra Kit (Life Technologies), according to the manufacturer's manual. RNA was recovered and purified by phenol:chloroform extraction and the integrity validated by gel electrophoresis.
mRNAs were microinjected into the cytoplasm of zygotes 20-22 h post hCG injection using a Piezo-drill (Eppendorf) and Eppendorf transferman micromanipulators. Hdac1 wild-type/mutant mRNA (500 ng/µl) or Kdm5b mRNA (1 μg/μl) was loaded into a microinjection needle and a constant flow was adjusted in order to achieve successful microinjection. Around 10 pl of mRNA was microinjected into the cytoplasm of zygotes or two-cell blastomere.
Immunofluorescence
Embryos were collected, washed and fixed for 10 min at room temperature in 4% formaldehyde in PBS (mouse, 10 min; cattle, 30 min). Embryo permeabilization was performed by treating with 0.5% Triton X-100 in PBS for 30 min in mice and 1 h in cattle. Embryos were then blocked in 10% fetal bovine serum in PBS for 1 h, and incubated in primary antibody for at least 1 h or 4°C overnight. After three washes with 0.1% Triton X-100 in PBS, fixed embryos were incubated with secondary antibodies for at least 1 h. DNA was counterstained with DAPI to locate the nuclear region. Embryos were imaged with Zeiss LSM 880 confocal microscope. All the antibodies used in the present study are listed in Table S1.
RNA-seq library preparation and sequencing
Mouse early two-cell and late two-cell stage embryos (50 embryos/sample, n=2 biological replicates/group) were collected at 24 h and 36 h post-fertilization, respectively. Bovine eight- and 16-cell stage embryos (30 embryos/sample, n=2) were collected at 72 h post fertilization. Total RNA was extracted using Arcturus PicoPure RNA Isolation Kit (Life Technologies) according to the manufacturer's manual. Then mRNAs were separated using oligo(dT)25 beads. Sequencing libraries were constructed by using NEB Next Ultra RNA Library Prep Kit for Illumina (E7530) based on the instruction. Briefly, mRNAs were fragmented and reverse transcribed. The cDNA library underwent end repair, poly(A)-tailing, adaptor ligation and PCR amplification for 12-15 cycles in order to prepare the sequencing libraries. Paired-end 150 bp sequencing was performed on a NovaSeq (Illumina) platform by Novogene.
ULI-NChIP-seq library preparation and sequencing
Mouse late two-cell embryos (120 embryos/sample, n=2) were collected at 48 h after hCG injection. The zona pellucidae of the embryos were removed with Acid Tyrode's solution, then the embryos were washed three times in 0.5% bovine serum albumin (Millipore) in DPBS (Gibco) before flash-freezing in liquid nitrogen. ULI-NChIP was performed according to the published protocols with slight modifications (Brind'Amour et al., 2015). One microgram of H3K4me3 (Cell Signaling Technology, 9751) or H3K27ac (Active Motif, AM39133) antibody was used for each immunoprecipitation reaction. ULI-NChIP-seq libraries were generated using the NEB Ultra DNA Library Prep Kit (E7645) according to the manual. Paired-end 150 bp sequencing was performed on a NovaSeq (Illumina) platform by Novogene.
RNA-seq data alignment and quantification
The raw sequencing reads were trimmed with Trimmomatic (version 0.39) (Bolger et al., 2014) to remove low-quality reads and adaptor sequences. Clean reads were then mapped to mm10 (mouse) or ARS-UCD1.2 (bovine) with Hisat2 (version 2.1.0) (Kim et al., 2015). The raw read counts of genes were calculated with featureCounts (version 1.6.3) (Liao et al., 2014). For quantification of repetitive element, repeat annotations were downloaded from UCSC genome browser, and the raw counts of repetitive element were counted with BEDTools (version 2.30.1) (Quinlan, 2014). Gene expression values were normalized to FPKM with Cufflinks (version 2.2.1) (Trapnell et al., 2012) for heatmap, line plot and bar plot visualization. The raw counts of repetitive elements were normalized with DESeq2 (Love et al., 2014) and then used for box pot visualization.
Differential expression analysis and functional annotation
Reads counts from FeatureCounts (genes) or BEDTools (repetitive elements) were used for differential expression analysis with DESeq2. The differentially expressed genes or repetitive elements between groups were identified when fold change >2 or <0.5 and adjusted P-value <0.05. The differential expression of repetitive elements was further validated by TEtranscripts python package (version 2.2.1) (Jin et al., 2015). Functional annotation and enrichment analysis of differentially expressed genes was performed with the Database for Annotation, Visualization and Integrated Discovery (DAVID) (Huang et al., 2009a; Huang et al., 2009b).
Classification of gene sets and clustering analysis
RNA-seq data for mouse and bovine oocytes and embryos was obtained from Jingyi Wu et al. (2016) and Alexander Graf et al. (2014), respectively. The FPKM matrixes were generated with Cufflinks as described in the RNA-seq data alignment and quantification section. Then the non-zero FPKM of genes was used for k-means clustering, which was performed with R based on Euclidean distance. After k-means clustering, genes with similar expression pattern were classified as maternal, minor ZGA, major ZGA or developmental genes. Clustering of differentially expressed genes was also performed as described above. For hierarchical clustering of mouse late two-cell embryos in different groups, euclidean distance of samples was calculated based on FPKM matrixes, and hierarchical clustering was performed with hclust (dist, method=‘average’) in R (http://www.rproject.org).
ChIP-seq data processing
The raw sequencing reads were trimmed with Trimmomatic (version 0.39) (Bolger et al., 2014) to remove residual adapter sequences and low-quality reads. Then the clean reads were aligned to mm10 using Bowtie2 (version 2.3.5) (Langmead and Salzberg, 2012) with default parameters. Alignments with low mapping quality were discarded by SAMtools (version 1.7) (Li et al., 2009), and PCR duplicates were removed with Picard (version 2.23; https://broadinstitute.github.io/picard/). To visualize H3K4me3 and H3K27ac signal in IGV genome browser, the genome was binned into 50 bp windows and RPKM for each window was calculated using bamCoverage function from DeepTools (Ramirez et al., 2014).
Correlation between biological replicates
The genome was binned into 2000 bp windows using the ‘makewindows’ utility of BEDTools (Quinlan and Hall, 2010), and the reads counts of each window were calculated and normalized to RPKM. Then the RPKM value for each sample was used to calculate Pearson correlation coefficient and draw scatter plots with R (http://www.rproject.org).
ChIP-seq peak calling
The peak number is affected by sequencing depth, and equal numbers of ChIP and input reads result in the best performance of peak callers (Chen et al., 2012; Jung et al., 2014). Therefore, we merged alignments of biological replicates and randomly subsampled them to equivalent depth using SAMtools before peak calling. The peaks were called by MACS2 (version 2.2.7.1) (Zhang et al., 2008) using the following parameters: -B -p 1e-5 --nomodel --broad --extsize 73. The peaks with a fold enrichment of less than two were removed. The filtered peaks were then annotated using ChIPseeker (Yu et al., 2015).
Identification of H3K27ac and H3K4me3-gained or -lost regions
To compared H3K27ac or H3K4me3 signals between groups, the mouse genome was scanned using a sliding window of 5 kb and step size of 1 kb, and then RPKM for each window was calculated. Next, H3K27ac or H3K4me3 signals were compared parallelly between H2O and Hdac1 MU. The H3K27ac and H3K4me3-gained or -lost regions were identified with following threshold: log2 (fold change) >1.5 for gained regions or <−1.5 for lost regions, and sum of RPKM in H2O and Hdac1 MU >1. Then the selected regions were merged with BEDTools if they were overlapped, and the genome coverage of merged regions was calculated by genomecov mode of BEDTools.
Motif analysis
The mouse genome was first scanned using a sliding window of 1 kb and step size of 1 kb, and then RPKM of H3K27ac for each window was calculated. Next, the ratio of H3K27ac signal in Hdac1 MU to H2O was calculated, and windows with the top 2% of ratios were retained as the most increased H3K27ac regions. On the other hand, windows with the top 2% of ratios of H2O to Hdac1 MU were retained as the most decreased H3K27ac regions. The selected increased or decreased regions were merged with BEDTools if they were overlapped. The merged regions were used for motif identification by Homer (version 4.11; http://homer.ucsd.edu/homer/motif/) with the parameters: -size given -p 10 -len 8.
Statistical analyses
Three biological replicates were conducted unless stated. Two-tailed unpaired Student's t-test were performed to compare differences between groups. Immunofluorescent intensities were measured with ImageJ. Figures were generated using GraphPad Prism 8 (GraphPad Software). P<0.05 indicates statistical significance. Results are shown as mean±s.e.m.
Acknowledgements
We thank Dr Li Shen and Dr Xudong Fu for their helpful discussions, and we thank Dr Zhe Zhang from Yuchun Pan's lab for their help on bioinformatic analysis.
Footnotes
Author contributions
Conceptualization: Y.D., K.Z.; Methodology: Y.D.; Formal analysis: S.L., P.Z., Y.S.; Investigation: Y.D., S.L., P.Z., L.X., L.W., Y.S., L.L.; Resources: H.W.; Data curation: S.L.; Writing - original draft: Y.D.; Supervision: S.W., K.Z.; Project administration: K.Z.; Funding acquisition: K.Z.
Funding
This project was supported by the National Natural Science Foundation of China (31672416, 31872348 and 32072731 to K.Z.; 31941007 to L.L. and S.W.; 32072939 to H.W.), the Zhejiang Provincial Natural Science Foundation (LZ21C170001 to K.Z. and LY19C180002 to H.W.) and the China Postdoctoral Science Foundation (2020M671742 to L.L.).
Data availability
Previously published datasets used in the present work include H3K27ac ChIP-seq (GSE72784; Dahl et al., 2016), H3K4me3 ChIP-seq (GSE71434; Zhang et al., 2016), mouse RNA-seq (GSE66582; Wu et al., 2016) and bovine RNA-seq (GSE52415; Graf et al., 2014). Our RNA-seq and ULI-NChIP-seq data have been deposited in GEO under accession number GSE182555.
Peer review history
The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.200854.
References
Competing interests
The authors declare no competing or financial interests.