Histone remodeling reflects conserved mechanisms of bovine and human preimplantation development

Abstract

How histone modifications regulate changes in gene expression during preimplantation development in any species remains poorly understood. Using CUT&Tag to overcome limiting amounts of biological material, we profiled two activating (H3K4me3 and H3K27ac) and two repressive (H3K9me3 and H3K27me3) marks in bovine oocytes, 2-, 4-, and 8-cell embryos, morula, blastocysts, inner cell mass, and trophectoderm. In oocytes, broad bivalent domains mark developmental genes, and prior to embryonic genome activation (EGA), H3K9me3 and H3K27me3 co-occupy gene bodies, suggesting a global mechanism for transcription repression. During EGA, chromatin accessibility is established before canonical H3K4me3 and H3K27ac signatures. Embryonic transcription is required for this remodeling, indicating that maternally provided products alone are insufficient for reprogramming. Last, H3K27me3 plays a major role in restriction of cellular potency, as blastocyst lineages are defined by differential polycomb repression and transcription factor activity. Notably, inferred regulators of EGA and blastocyst formation strongly resemble those described in humans, as opposed to mice. These similarities suggest that cattle are a better model than rodents to investigate the molecular basis of human preimplantation development.

Keywords: cattle; epigenetics; genome activation; preimplantation embryo.

Figure 1. Profiles of histone modifications in bovine oocytes and preimplantation embryos

1. Principal components analysis (PCA) of CUT&Tag libraries.

2. Genome coverage of peaks identified in both biological replicates. Peaks classified as genic (overlapping 2 Kb promoter or gene bodies) or intergenic. For H3K4me3 and H3K27ac, genome coverage calculated for broad peaks until the 8‐cell stage, then for narrow peaks.

3. Normalized signal (counts per million; CPM) of one biological replicate per developmental stage and histone mark. Shaded regions correspond to noncanonical broad distributions (green, blue) and intragenic enrichment (red, orange). Viewing range from 0 to 1.5 CPM. Values exceeding maximum range indicated by pink bars.

4. Average normalized signal (CPM) for each histone mark at transcription start sites (TSS) and end sites (TES). Transcripts scaled to 2 Kb, 3 Kb upstream and downstream regions shown.

Data information: GV oocytes (GV), 2‐cell (2C), 4‐cell (4C), 8‐cell embryos (8C), α‐amanitin‐treated 8‐cell embryos (8A), morula (M), blastocyst (BL), trophectoderm (TE), and inner cell mass (ICM).

Figure 2. Characterization of partially methylated domains (PMDs) in GV oocytes

1. Clustering of PMDs based on normalized signal of H3K4me3, H3K27ac, H3K27me3, H3K9me3, and ATAC‐seq in GV oocytes. Signal normalized by CPM in 100 bp windows, with PMDs scaled to 3 Kb, and showing regions 1 Kb up‐ and downstream.

2. Representative gene track image of a bivalent PMD (cluster_1), overlapping three homeobox genes: CDX2, GSX1, and PDX1. One biological replicate shown per stage and mark. Viewing range from 0 to 1.5 CPM. Values exceeding maximum range indicated by pink bars.

3. Average epigenetic signal (CPM ± standard error) for PMDs belonging to each cluster across development.

Figure EV1. Identification of noncanonical broad H3K4me3 and H3K27ac relative to partially methylated domains (PMDs)

1. Comparison of peak sets from biological replicates identified using either “broad” or “narrow” peak calling parameters. Jaccard statistic measures base‐pair overlap between two peak sets, ranging from 0 (no overlap) to 1 (complete overlap).

2. Gene track image centered on the CDX2 locus, with representative PMDs predicted from GV oocyte CpG methylation data (DNAme). Highlighted regions correspond to PMDs and overlap with broad H3K4me3 and H3K27ac domains. Viewing range from 0 to 1.5 CPM. Values exceeding maximum range indicated by pink bars.

3. Average normalized signal (CPM) at PMDs for each histone modification and ATAC‐seq. PMDs scaled to 3 ± 1 Kb up‐ and downstream.

4. Percent of PMDs overlapped (by at least 1 bp) by broad peaks of a given histone modification.

5. Percent of broad peaks for a given histone modification overlapped (by at least 1 bp) by PMDs.

Figure 3. H3K9me3 and H3K27me3 mark gene bodies in pre‐EGA embryos

1. Average normalized signal (CPM) of histone marks and chromatin accessibility for each transcript cluster, defined by k‐means clustering (k = 3) based on normalized H3K27me3 and H3K9me3 signal (CPM) across development. Loci scaled to 2 Kb, regions 3 Kb upstream and downstream shown.

2. An example of “cluster_1” loci: the HOX cluster.

3. An example of “cluster_2” loci: GRK2 (a kinase). One biological replicate shown per stage and mark. Viewing range from 0 to 1.5 CPM. Values exceeding maximum range indicated by pink bars.

4. Gene body CpG methylation (%) of “cluster_1” and “cluster_2” genes throughout development.

5. For each cluster, comparison of change in gene body CpG methylation (GV oocytes versus 2‐ to 4‐cell embryos) to change in intragenic H3K27me3 signal (GV oocytes versus 2‐cell embryos). Trends summarized in red text.

Figure 4. Chromatin accessibility precedes establishment of canonical H3K4me3 and H3K27ac in morula

1. Average ATAC‐seq, H3K27ac, and H3K4me3 signal (CPM) at H3K27ac peaks which appeared during the 8‐cell to morula transition. Peaks classified as genic if they overlapped gene bodies or promoters (2 Kb upstream of TSS).

2. Average signal at TSS.

3. Normalized ATAC‐seq, H3K27ac, and H3K4me3 signal for one biological replicate per developmental stage at the LAPTM4A, NANOG, and DPPA3 loci. Shaded regions correspond to putative enhancers (E; blue) and promoters (TSS; orange). Viewing range from 0 to 1.5 CPM. Values exceeding maximum range indicated by pink bars.

4. Enrichment of pioneer transcription factor motifs in H3K27ac and ATAC‐seq peaks at each stage. From the GV to 8‐cell stage, broad H3K27ac peaks (*) were used for motif enrichment analysis, then narrow peaks were used.

Figure 5. Effect of transcription inhibition on the epigenetic profile of 8‐cell embryos

1. Normalized signal (CPM) in 4‐cell (4C), 8‐cell (8C), and 8‐cell embryos cultured in the presence of α‐amanitin (8A) at regions that gained, lost, or retained peaks during the 4‐ to 8‐cell transition. Peaks scaled to 500 bp (±500 bp).

2. Genomic coverage of regions that gained, lost, or retained peaks. Shading reflects the impact of transcription inhibition on peak status. Red indicates regions with a change in peak status in 8A embryos compared to 8C.

3. Average signal (CPM) of active epigenetic marks at all TSS in 2‐cell (2C), 4C, 8C, and 8A embryos.

4. Normalized H3K27ac signal at “cluster 3” TSS (±2 Kb), identified based on H3K27ac signal at TSS in 2C, 4C, 8C, and 8A embryos (Fig EV2), and a gene track view of representative locus, CDX2.

5. Normalized H3K27ac signal (CPM) at “cluster 1” TSS (±2 Kb) and a gene track view of a representative locus, OTX2. One biological replicate shown per stage and mark. Viewing range of gene tracks from 0 to 1.5 CPM. Values exceeding maximum range indicated by pink bars.

6. Enrichment of select transcription factor motifs in intergenic H3K27ac peaks that were gained or lost during the 4‐ to 8‐cell transition, and which were either sensitive or robust to inhibition of embryonic transcription.

7. Normalized expression of DUX4 (ENSBTAG00000049205) across development. Boxplots indicate the median and interquartile range (IQR), and whiskers span 1.5 times the IQR. Data points indicate biological replicates (n = 3 for GV, 4C, 8C, 8‐16C, 16C, BL; n = 4 for MII and 8A).

Figure EV2. Epigenetic and transcriptomic perturbations in transcription‐inhibited 8‐cell embryos

1. A–D

   Normalized expression (CPM) of writers, readers, and erasers for (A) H3K4me3, (B) H3K27me3, (C) H3K9me3, and (D) H3K27ac in MII oocytes, 8‐ to‐16‐cell embryos cultured in the presence of α‐amanitin (8A), and control 8‐ to 16‐cell embryos (8C). Genes marked “*” were differentially expressed between 8A and 8C (adjusted P < 0.01, log2FC > 2).

2. E

   Genome coverage of 8A peaks that were not present in 4C or 8C embryos.

3. F

   K‐means clustering (k = 7) of TSS based on normalized H3K27ac signal (CPM) in 2C, 4C, 8C, and 8A embryos.

4. G

   Genome coverage of peaks gained or lost during the 4‐ to 8‐cell transition, subcategorized based on whether they were present in 8A embryos. Peaks classified as either genic (overlapping gene bodies or 2 Kb promoters) or intergenic.

Figure 6. Epigenetic shifts that underscore lineage segregation in the blastocyst

1. A

   Conserved mechanisms of TE and ICM specification and maintenance in mammals.

2. B

   Species‐specific differences in TE and ICM markers in morula and blastocysts. The timing of each stage is indicated by days post fertilization (dpf).

3. C

   Epigenetic profiles of key pluripotency genes and TE‐specific markers. Viewing range from 0 to 1.5 CPM. Values exceeding maximum range indicated by pink bars. All biological replicates shown.

4. D

   For each histone modification, proportion of peaks unique to the ICM, unique to the TE, or shared in common between ICM and TE that were already present in morula.

5. E

   Chromatin state predictions based on chromatin accessibility and histone modification data. Emission probabilities indicate the likelihood of a given mark occurring in a given state.

6. F–I

   Genome coverage of (F) polycomb repression, (G) strong active TSS and strong active elements, (H) weak active elements, and (I), active elements by developmental stage.

7. J

   Motif enrichment of selected regulators in strong active elements in M, ICM, and TE.

8. K

   Average normalized expression of selected regulators in 16‐cell embryos (16C) and blastocysts.

Figure EV3. Epigenetic signatures of blastocyst lineages

1. Epigenetic profiles of EOMES (mouse‐specific TE marker), PLAC8 (human‐specific TE‐marker), and the locus coding for pregnancy‐associated glycoproteins (PAG).

2. Normalized expression (CPM) of EOMES and PLAC8.

3. Genome coverage of peaks unique to ICM, TE, or M, or shared in common between ICM and TE. Peaks categorized based on presence or absence in morula.

Figure EV4. Dynamic activity of bovine‐specific long terminal repeats (LTRs)

1. Average enrichment of select chromatin states at LTRs across development, measured as log ratio of random to expected overlap of repeats with a given chromatin state (LR ± standard error). LTRs subcategorized based on repeat family (e.g., ERV1) and whether they were Bos taurus‐specific (BT‐specific).

2. Proportion of LTRs in each subcategory that were differentially expressed (DE) between the 4‐ and 16‐cell stages (adjusted P < 0.05, log2FC > 1).

3. Expression (z‐score) and chromatin state (LR) dynamics of highly abundant DE BT‐specific LTRs.

Figure EV5. Differential regulation of blastocyst lineages

1. Differentially regulated genes (adjusted P < 0.05, log2FC > 1) in ICM and TE based on H3K27me3 and H3K27me3 signal in 2 Kb promoters. Enriched functional terms and corresponding genes reported for genes that were differentially activated (DAG) in one cell type and differentially repressed (DRG) in the other.

2. Normalized expression of KLFs.

3. TFs with differential activity (P < 0.05) in ICM compared to M, based on footprinting analysis of ATAC‐seq data.

4. Average ATAC‐seq signal in ICM and M at JUNB footprints (JASPAR database motif MA0490.1).

5. Normalized expression of JUN and FOS factors.

6. Average ATAC‐seq signal in ICM and M at DUX4 footprints (JASPAR database motif MA0468.1).