Reprogramming histone modification patterns to coordinate gene expression in early zebrafish embryos

Overviews of histone modification dynamics in early zebrafish embryos

We profiled genome-wide H3K4me3 and H3K27me3 occupancy in zebrafish gametes and early embryos (including 2-cell, 16-cell, 128-cell and 1 K-cell stages) with ChIP-seq. For each stage, data from two replicates (Additional file 1: Figure S1a) were combined with equivalent read pairs. Peaks were called and then annotated around the TSS to identify associated genes. Compared with published ChIP-chip data for sperm [10] and the 1 K-cell stage [15], our data uncovered many more genes and covered the majority of genes found by previous studies (Additional file 1: Figure S1b). We noticed that the peak widths (Additional file 1: Figure S2) and CpG dinucleotide density (Additional file 1: Figure S3) in promoter peak regions are higher than those in distal peak regions for both H3K4me3 and H3K27me3.

To provide an understanding of the global reprogramming progress of histone modifications in early zebrafish embryos, all genes marked by histone modifications were categorized into five groups (A-E) primarily according to their states in gametes (Fig. 1a). To our surprise, 23.9% of marked genes are in a co-occupied state in sperm (i.e., group A), and this proportion ranks highest among all of the investigated stages. Moreover, only 204 of those co-occupied genes, which are developmental regulators, e.g., shha (Fig. 1c), sox2 (Additional file 1: Figure S4), pou5f3 (Additional file 1: Figure S4) and hox clusters (Additional file 1: Figure S5), have a stable histone modification status throughout all of the investigated stages. Gene group C, which contains genes in an active state in sperm, represents nearly half (45.7%) of all marked genes. Another 22% of marked genes are embryo specific genes (group E). Repressed genes in sperm (group B) and oocyte-specific genes (group D) represent 5.3 and 2.8% of marked genes, respectively.

We also compared the histone modification dynamics in zebrafish with those of early mouse embryos (Additional file 1: Figure S6). The global patterns are approximately stabilized at the 2-cell stage, where the mouse zygotic genome initializes activation. Genes are divided into five groups depending on their states in gametes, and similar reprogramming characteristics are observed. In mouse, gene groups A to E represent 11.3, 15.9, 30.7, 30 and 12.1% of all marked genes, respectively. As in zebrafish, sperm is the main source of genes marked both by H3K4me3 and H3K27me3. For genes in group A, their co-occupancy is lost rapidly by conversion to an active state, except for 30 genes remain unchanged. The status of active genes (group C) is stabilized in the 2-cell stage. The largest difference between mouse and zebrafish is found in gene group D, which contains the oocyte-specific modified genes. Group D represents nearly as many genes as the active group (group C) in mouse, but only represents 2.8% in zebrafish. Furthermore, we examined the proportions of homologous genes in corresponding gene groups between zebrafish and mouse and found that the majority of genes in groups A and C are homologous (Fig. 1b).

Previous studies of DNA methylation have revealed larger methylation-level differences between gametes in mouse than in zebrafish. In our histone modification data, stage-wise dissimilarities (Additional files 2 and 3) also suggested that there are sperm-oocyte differences in histone modification status between mouse (0.66) and zebrafish (0.44). Taken together with the broad H3K4me3 domains reported in mouse oocytes but not in zebrafish oocytes, we speculate that the histone modification reprogramming processes in zebrafish and mouse share similarities but differ in detail.

Histone modification patterns are established prior to gene expression

Histone modifications could be indicators of gene expression activities. In our data, we found that gene expression and the gene histone modification intensity of H3K4me3 in promoter regions are globally positively correlated (Fig. 2a). Moreover, the correlation efficiency increased for H3K4me3 modifications as early embryo development proceeded (Additional file 1: Figure S7). However, for H3K27me3, there are weak but negative associations between histone modification intensity in promoter regions and gene expression (Additional file 1: Figure S7).

Prior to ZGA, gene expression levels are static, as they are mostly driven by maternal deposition, and expression levels increase upon the onset of ZGA [22]. To further characterize histone modification dynamics during ZGA, we examined histone modification patterns for genes that are expressed in the first wave in early ZGA of zebrafish embryos. Previously studies has uncovered a set of genes that expressed at the very beginning of ZGA [23]. Here we found that those early expressed genes largely consist of co-occupied and active gene groups (Fig. 2b) along with members from other groups. These early expression genes lack H3K27me3 modifications at the 1 K-cell stage (Additional file 1: Figure S1c). Histone modification patterns of key developmental regulators that are initialized at ZGA in zebrafish, including pou5f3, the soxB gene family and nanog, were examined, and we found that histone modification co-occupancy on pou5f3, sox19b and sox2 (Additional file 1: Figure S4) are stable across the investigated stages. However, this result did not hold true for nanog (Fig. 2e), which was reported to induce 70% of the first wave gene s [23], H3K4me3 enrichment has long been established early from gametes and embryos, but H3K27me3 modifications are removed after the 2-cell stage to match the pattern of sperm and then are restored after ZGA (48hpf stage). Another example of an early expressed gene with loss of repressive markers during ZGA is lft1 (Additional file 1: Figure S4). lft1 is an inhibitor of the nodal signaling pathway, which is responsible for mesoderm induction and L-R axis development [24, 25]. H3K4me3 enrichment is accompanied by loss of H3K27me3 enrichment around ZGA. However, for those genes induced by nanog, e.g., cldnb (Additional file 1: Figure S4), oep and nnr (Fig. 2d), H3K4me3 enrichment is initialized at the 1 K-cell stage, and a clear peak emerged at the dome stage. Overall, histone modification patterns of genes expressed in the first wave are converted to co-occupied or active states prior to gene expression (Fig. 2c). These findings are in line with the results from published data demonstrating that histone modification patterns are established prior to ZGA [15].

Distinct DNA methylation patterns in distal and promoter peaks

It is known that the DNA methylome of early embryos undergoes extensive reprogramming. To investigate the relationship between DNA methylation and histone modification reprogramming in the early embryo, we first examined the correlations between DNA methylation levels and histone enrichments around TSS and found that there are weak but negative correlations between them for both H3K4me3 and H3K27me3 (Additional file 1: Figure S8a). We then checked the histone modification status of the DNA methyltransferase gene family (Additional file 1: Figure S8b) and found that there are no repressive H3K27me3 peaks found in their promoter regions before ZGA, which coincides with their expressions [22]. Additionally, we noticed that there are dynamics in the proportion of all peaks overlapping with hypomethylated regions. The proportion first decreased in early embryo stages and then gradually increased in later stages (Fig. 3a). After examination of the DNA methylation levels in distal and promoter peak regions, we concluded that hypermethylated (average methylation level not less than 0.7) promoter peaks are responsible for the reduction in the proportion of hypomethylated peaks in early embryos, apart from significant differences in DNA methylation levels between distal and promoter peaks (Fig. 3b). Moreover, we also observed that unlike zebrafish embryos, there are minimal hypermethylated peaks in mouse, and DNA methylation levels decreased from gametes to early embryos for both distal and promoter peaks (Additional file 1: Figure S9). We wondered if the disappearance of these hypermethylated promoter peaks is related to DNA methylation reprogramming. We examined DNA methylation dynamics in hypermethylated promoter peaks found in the 16-cell and 128-cell stages. To our surprise, DNA methylation levels in those hypermethylated promoter peaks do not decrease in later stages (Fig. 3c). By further checking the histone modification dynamics for those genes with hypermethylated promoter peaks, we found that the histone modification patterns of associated genes are distinct to adjacent stages (Fig. 3d). Therefore, these hypermethylated promoter peaks mark stage-specific genes and the reprogramming process is independent of DNA methylation reprogramming.

Gene age preferences of histone modifications in early embryos

Previously, we observed that during DNA methylome reprogramming in zebrafish and mouse, evolutionarily young genes tend to undergo extensive reprogramming compared to older ones (unpublished data). We asked whether there are preferences in gene age in the context of histone modification states, as new genes evolved to render novel functions during species evolution. To address this question, a phylostratigraphic approach [31] was used. In this method, protein coding genes of zebrafish and mouse were assigned to a certain phylostratum which represented the gene age [32]. By comparing the proportion of genes assigned to a certain phylostratum in selected gene groups to the proportion in the whole genome, we could find over- and under-represented phylostratum in each group.

For the co-occupied gene group in zebrafish (Fig. 4c), the proportion of phylostrata 4 and 5 was significantly higher than that in the whole genome, while the proportions of the first two and last two phylostrata were significantly lower. In mouse, we find that the proportions of phylostrata 3–5 were higher and the proportions of phylostrata 1, 6, 12, 16, and 19 were lower than those in the whole genome (Fig. 4e). Genes originating from organisms that evolved during phylostratum 4 (opisthokonta to holozoa) and phylostratum 5 (holozoa to metazoa), which are multicellular species, are more likely to contain co-occupied modifications in both mouse and zebrafish.

For the active gene groups, genes originating in primitive unicellular organisms (phylostratum 1) and eukaryotic organisms (phylostratum 2) were significantly enriched. The proportion of genes that diverged from holozoa organism onward (phylostrata 5–14) are significantly lower in zebrafish (Fig. 4d). In mouse, which is similar to zebrafish, genes from first three phylostrata, which are the oldest genes, are more likely to contain active marks than the younger genes (phylostratum 5 onward) (Fig. 4f). In contrast to co-occupied and active gene groups, gene groups B, D and E contain more young genes. The gene age preferences for histone modification status were also confirmed in each stage (Additional file 8: Table S7and Additional file 9: Table S8). We further examined the gene age preferences for the active and co-occupied gene groups as well as all genes in each stage and found that opisthokonta, in which multicellular organisms originated, is the key boundary for gene age preference. Genes that diverged before opisthokonts are more likely to be favored by active modifications, while genes that arose between opisthokonta and metazoa are preferred by co-occupied histone modifications. In summary, gene age preference reflects the phylotypic phase [31] during the transition from unicellular to multicellular organisms.

Zebrafish stock

Wild-type zebrafish line TU, maintained in our own laboratory, was used in this study. Zebrafish were raised for three to twelve months under standard conditions (28 °C, 14 h light / 10 h darkness) before experimental materials were collected.

Zebrafish sperm collection

Adult male fish were first anesthetized in system water containing 0.016% Tricaine. Testes were dissected from adult male fish as described [37] and placed in 10 ml PBS buffer (filtered by 0.22um filter) for 10 min at room temperature. During that time, the testes were pipetted several times. Tissue fractions were removed with 75 nm filter units before they were transferred to 15 ml tubes. The tubes were centrifuged for 1 min at 300 g and then held still for 3 min. After that, the upper half of the liquid was transferred to a new 10 ml tube, which was filled with PBS to 10 ml. The last step was repeated again. The sperm purity was checked under a microscope, and additional steps were performed if the purity was not sufficient. Aliquot of approximately 1 million sperm per tube were made for one ChIP experiment. The tubes were centrifuged for 10 min at 10,000 g, and the supernatant was discarded. Cross-linking or snap freezing of the sperm sample in liquid nitrogen was performed before storage at − 70 °C.

Oocyte and early embryo collection

Approximately 5 min after spawning three washes with system water, unfertilized oocytes were collected by squeezing the abdomen of females. The squeezed females were transferred to separate tanks for recovery. Embryos were collected for 10–15 min after each round of spawning and washed with system water three times. Collected unfertilized oocytes and embryos were washed three times with system water to remove potential follicle cells or other contaminants and then were transferred to 90 mm glass petri dishes with approximately 10 ml system water before proceeding to next step. Embryos were grown at 28 °C and staged according to standard morphological criteria [38], Three minutes prior to the desired embryonic stage or immediately after the collection of oocytes, prewarmed (approximately 28 °C) pronase solution (2 mg/ml at final concentration) was added to oocytes and embryos to perform dechorionation. The dishes were placed on ice to prevent embryos and oocytes from developing further, the chorion was disrupted with a glass dropper, and then the embryos were washed three times with cold (4 °C) HBSS buffer to remove chorion fragments. One thousand embryos or oocytes (for one ChIP experiment) were transferred to 1.5 ml low-binding tubes with a glass dropper. One ml of cold (4 °C) deyolk buffer [18] were add to each tube after removing the HBSS buffer, and then the tube was rotated for 3 min at 4 °C. After centrifugation for 3 min at 1000 g, the supernatant was discarded. We proceeded directly to cross-linking steps or snap freezing the samples in liquid nitrogen before storage at − 70 °C. After collecting eggs and embryos, fish were raised in separated tanks until full recovery.

Cross linking

Tubes with pellets were filled with cold PBS buffer to 1 ml, 37% formaldehyde solution (sigma) was added to a final concentration of 1%, and then tubes were rotated for 10 min at room temperature. Cross linking was halted by adding fresh glycine solution to a final concentration of 125 mM, and tubes were rotated for another 5 min. The cross-linked cells were centrifuged for 10 min at 10,000 g and washed two times with cold PBS buffer. The supernatant was removed, and cell pellets were stored at − 70 °C or directly used in chromatin immunoprecipitation steps.

Chromatin immunoprecipitation and library preparation

The bead-antibody complexes were prepared before cell lysis using 2 μg of H3K4me3 (Abcam ab8580) antibody or 2.5 μg of H3K27me3 (Millipore 07–449) antibody per 10 μl of protein A Dynabeads in low binding PCR tubes. Beads were washed two times with 100 μl ChIP buffer (10 mM Tris-HCl, pH 7.5; 1 mM EDTA; 0.5 mM EGTA; 0.1% Sodium Deoxycholate; 1% Triton X-100; 0.1% SDS; 140 mM NaCl; fresh proteinase inhibitor cocktail), then the total volume was brought up to 100 μl per assay. The bead-antibody complexes were incubated on a rotator at 4 °C for 2 h. The ChIP buffer was removed before adding the sonicated chromatin. Cell were lysed with 100 μl of lysis buffer (0.2% SDS; 10 mM Tris-HCl buffer, pH 8.0; 10 mM EDTA, pH 8.0; fresh proteinase inhibitor cocktail) on ice for 30 min. Lysed cells were sonicated in Bioruptor for 25 cycles (30s on and 30 s off per cycle) to achieve fragments that were 300–500 bp in length. The tubes were centrifuged for 10 min at 4 °C (> 10,000 rpm), and 100 μl of supernatant was transferred to the bead-antibody complex tubes. The tubes were rotated slowly at 4 °C overnight (8–10 h). The next day, the tubes were placed on a magnetic rack for one minute, and the supernatant was removed. The beads were washed two times with low salt buffer (10 mM Tris-HCl, pH 8.0; 1 mM EDTA; 0.5 mM EGTA; 1% Triton X-100; 0.1% SDS; 250 mM NaCl) and high salt buffer (10 mM Tris-HCl, pH 8.0; 1 mM EDTA; 0.5 mM EGTA; 1% Triton X-100; 0.1% SDS; 500 mM NaCl); one time with LiCl buffer (10 mM Tris-HCl, pH 8.0; 1 mM EDTA; 0.5 mM EGTA; 1% Triton X-100; 0.1% SDS; 300 mM LiCl). The supernatant was removed, and TE buffer was added (20 mM Tris-HCl, pH 8.0; 10 mM EDTA, pH 8.0) before the beads were transferred to another new tube. The supernatant was discarded, and 100 μl of elution buffer was added (20 mM Tris-HCl, pH 8.0; 5 mM EDTA, pH 8.0; 50 mM NaCl; 1% SDS; 2 mg/ml proteinase K). The tubes were placed on a thermal shaker operated at 68 °C and 1000 rpm for 6–8 h. DNA fragments were purified with AMPure XP beads (Beckman) directly before library preparation following Illumina’s instructions. For each sampled stage, two ChIP replicates were performed, and libraries were prepared individually for each replicate except for H3K27me3 ChIP with oocytes. For the H3K27me3 ChIP assay performed in oocytes, precipitated DNA fragments were pooled from two replicates to provide enough DNA for library preparation.

Date processing: read mapping, peak calling and annotations, visualization

ChIP-Seq data generated in this study or public ChIP-Seq data (dome and 48 hpf stages) used in this study were processed as following. Low quality and adapter-containing reads were trimmed using Trimmomatic [39]. High-quality reads then were fed into Bowtie2 [40] for mapping against the zebrafish genome (Zv9). Uniquely mapped read pairs were kept for downstream analysis. Duplicated reads were removed with Picard tools (http://broadinstitute.github.io/picard). The replicates of each stage were merged with equivalent read pairs and peaks were called using MACS2 [41] with default settings (except that ‘--broad’ was specified) and filtered with a qValue cutoff set to 0.05. For mouse, data for called peaks were directly downloaded from the NCBI GEO database. Peak annotations were done using the R package ChiPpeakAnno [42] against Ensembl gene annotations for zebrafish genome Zv9 or mouse genome mm9, assigning each peak to the nearest transcription start site. Those peaks that fell between 3 kb upstream and 3 kb downstream of the TSS were kept for further analysis. For data visualization, we adjust ChIP-seq data with internal reference to account for enrichment fold differences among stages. We choose peak that overlap with promoter of eef1a1b, which is a translation elongation factor that is reported stable in various developmental stages [43]. ChIP-seq track screenshots are plotted with Gviz [44] package in R environment.

GO, INTERPRO and gene Phylostratum enrichment analysis

Gene ontology and INTERPRO domain enrichment analyses were carried out using the DAVID web service [45, 46] with annotated Ensemble gene IDs. The enriched items were filtered with the FDR cutoff set to 0.01. Gene phylostratum information was acquired from previously published data [32], and enrichment analysis was performed with a hypergeometric test in R with the FDR cutoff set to 0.01 for determining significantly enriched phylostrata in each gene group.

Statistics test for significance

The Wilcox test was used in this study to determine the significant differences between groups unless otherwise specified in the manuscript. The significance level was uniformly set to 0.01.