Epigenomic Landscape of Human Fetal Brain, Heart, and Liver*

The epigenetic regulation of spatiotemporal gene expression is crucial for human development. Here, we present whole-genome chromatin immunoprecipitation followed by high throughput DNA sequencing (ChIP-seq) analyses of a wide variety of histone markers in the brain, heart, and liver of early human embryos shortly after their formation. We identified 40,181 active enhancers, with a large portion showing tissue-specific and developmental stage-specific patterns, pointing to their roles in controlling the ordered spatiotemporal expression of the developmental genes in early human embryos. Moreover, using sequential ChIP-seq, we showed that all three organs have hundreds to thousands of bivalent domains that are marked by both H3K4me3 and H3K27me3, probably to keep the progenitor cells in these organs ready for immediate differentiation into diverse cell types during subsequent developmental processes. Our work illustrates the potentially critical roles of tissue-specific and developmental stage-specific epigenomes in regulating the spatiotemporal expression of developmental genes during early human embryonic development.

The epigenetic regulation of spatiotemporal gene expression is crucial for human development. Here, we present wholegenome chromatin immunoprecipitation followed by high throughput DNA sequencing (ChIP-seq) analyses of a wide variety of histone markers in the brain, heart, and liver of early human embryos shortly after their formation. We identified 40,181 active enhancers, with a large portion showing tissuespecific and developmental stage-specific patterns, pointing to their roles in controlling the ordered spatiotemporal expression of the developmental genes in early human embryos. Moreover, using sequential ChIP-seq, we showed that all three organs have hundreds to thousands of bivalent domains that are marked by both H3K4me3 and H3K27me3, probably to keep the progenitor cells in these organs ready for immediate differentiation into diverse cell types during subsequent developmental processes. Our work illustrates the potentially critical roles of tissue-specific and developmental stage-specific epigenomes in regulating the spatiotemporal expression of developmental genes during early human embryonic development.
The accurate spatial and temporal expression of genes is crucial for mammalian embryonic development, and epigenetic regulation is important for controlling the precise tissue-or developmental stage-specific expression patterns of those genes (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16). Enhancers are one of the most important types of regulatory elements for controlling spatiotemporal gene expression patterns (11,(17)(18)(19)(20)(21)(22)(23). It is well known that distal regulatory elements, such as active enhancers, are normally marked by the combination of H3K27ac and H3K4me1, whereas poised enhancers are marked by only H3K4me1 (24). H3K4me3 normally marks the promoter regions of active/permissive genes, whereas H3K27me3 marks the promoter and gene body regions of repressed genes (25)(26)(27)(28)(29). To better understand the regulatory mechanisms during early human embryonic development, it is necessary to analyze the global profiles of these crucial epigenetic markers in the genomes of different human embryonic tissues. However, because of the scarcity of the materials that are available for analysis, a comprehensive chromatin immunoprecipitation (ChIP)-seq analysis of early human embryos has not yet been achieved. We attempted to isolate tissues from three representative organs from the three germ layers, namely the brain, heart, and liver, from human embryos at 12 weeks of gestation, and we systematically analyzed the epigenetic markers H3K27ac, H3K4me1, H3K4me3, and H3K27me3 using a ChIP-seq approach. Moreover, we also analyzed the bivalent H3K4me3/H3K27me3 domains in these three organs via sequential ChIP-seq. As a control, we also downloaded the ChIP-seq data for the corresponding markers in human embryonic stem cells (hESCs) 7 and those for the adult human brain, heart, and liver from the ENCODE database (30,31). We detected highly reproducible patterns of dynamic epigenetic changes in both the enhancer regions and other functional genomic regions.

Experimental Procedures
Informed Consent and Ethics Approval of the Study-The approval of the Reproductive Study Ethics Committee of Peking University Third Hospital (research license 2012SZ013) was obtained in this study for collection and use of the embryos from the donors. All the aborted fetuses without personal iden-tifiers were obtained after the donors signed a written informed consent. The gestational age was determined by the date of the woman's last menstrual period together with the ultrasound scans.
Isolation of the Early Human Embryos-Anonymous embryos at 12 weeks of gestation were used in this study. The fetal tissues were further identified and isolated after a mechanical dissection, and all the samples were serially washed and then immediately frozen in liquid nitrogen. The biological replicates were prepared from the same tissue of an individual embryo before tissue lysis.
ChIP-seq Library Preparation-The isolated tissues were first washed several times with ice-cold phosphate-buffered saline (PBS) to remove any excessive blood and other possible contaminant cells. The tissues were further minced with a homogenizer in tissue lysis buffer (0.32 M sucrose, 5 mM CaCl 2 , 3 mM, 0.1 mM EDTA, 10 mM Tris-HCl (pH 8.0), 1 mM DTT, and 0.1% Triton X-100) on ice. The collected lysate was cross-linked with 1% formaldehyde for 10 min at room temperature and then quenched with glycine. ChIP was performed using an EZ-Magna ChIP HiSens chromatin immunoprecipitation kit (Millipore) following the manufacturer's instructions. Briefly, the precipitated nuclei were resuspended in SCW buffer (part of the chromatin immunoprecipitation kit) containing protease inhibitor mixture and were then sonicated using a Covaris S2, yielding an average fragment size of ϳ250 bp. The chromatin was then immunoprecipitated using protein A Dynabeads (Invitrogen) with 5 g of ChIP-grade antibody overnight at 4°C. The beads were further collected and washed several times with high and low salt buffers, and the DNA fragments were then de-cross-linked with proteinase K at 65°C for 4 h and purified via ethanol precipitation. We used a pre-immune IgG as the negative control for a more accurate peak-calling model.
For the H3K4me3/H3K27me3 sequential ChIP, the sonicated chromatin was first incubated with 20 g of anti-H3K27me3 antibody at 4°C overnight, and the beads were then collected, washed, and eluted with 10 mM DTT. The eluate was diluted 20-fold with SCW buffer and was then used for the second round of chromatin immunoprecipitation using 2.4 g of anti-H3K4me3 antibody.
All of the ChIP assays were performed using the following antibodies: H3K4me3 (Abcam, ab8580), H3K4me1 (Abcam, ab8895), H3K27me3 (CST, catalog no. 9733), and H3K27ac (active motif, catalog no. 39133). Other groups previously validated all of the antibody lots for their specificity and efficacy.
The ChIP-seq libraries of the ChIPed DNA and the matched IgG control DNA were prepared using the NEBNext DNA library prep master mix kit, indexed for multiplexed runs, and sequenced using the Illumina HiSeq 2000/2500 platform at 2ϫ 100 bp.
RNA Extraction and RNA-seq Library Preparation-The total RNA was extracted from the frozen tissues of the early gestation embryos using an RNeasy mini kit (Qiagen) following the manufacturer's standard protocol. Only RNA with a relatively high quality (RNA Integrity Number (RIN) Ͼ7.0 on an Agilent 2100 Bioanalyzer platform) was used for the RNA-seq libraries. The RNA-seq libraries were constructed using an Illu-mina TruSeq RNA prep kit and sequenced using the Illumina HiSeq 2000/2500 platform with 100-bp pair-end mode.
Published Dataset Accession-In this study, we took advantage of the available ChIP-seq and RNA-seq data. To avoid any possible confounding biases because of the different pipelines used across labs, we downloaded the raw fastq sequences (if available) and incorporated the data into our analysis pipelines for comparison. All of the accession numbers for the published dataset that we used in this study are listed in supplemental  Table S1.
Sequencing Data Quality Control-All of the de-multiplexed sequencing reads were first cleaned to remove any artificial sequences, such as sequencing adapters introduced during the experimental process, and reads with more than 10% low quality bases were also discarded. For the ChIP-seq samples, we calculated the fraction of reads in peaks (FRiP scores), and only the samples with FRiP scores no less than 10% in this study (regardless of the sequential ChIP-seq samples) were kept for further analysis (see supplemental Table S1).
RNA-seq Data Processing-For the RNA-seq data, the cleaned reads were aligned to the hg19 GTF reference genome using TopHat (version 2.0.9) with the default parameters. The gene expression levels (FPKM) of each sample were calculated using the cuffinks, cuffquant, and cuffnorm suite (32,33) with the transcript annotations downloaded from the UCSC Genome browser (34,35). Only the RefSeq genes, with FPKM being more than 1 were used in the subsequent analysis.
ChIP-seq Data Processing-For the ChIP-seq data, the trimmed reads were mapped to the human genome (hg19 assembly) using BWA aligner (version 0.7.5a-r405) with the options "-i 15 -q 10 -t 4" (36). Only the unique mapped reads were retained for further analyses.
MACS (version 1.4) was used for the peak calling with the options "-f BED -g hs -n chip_output -B -S." The input mapped reads files for the individual replicates were merged as recommended by MACS (version 1.4.2) (37). Specifically, to obtain a precise result for the enhancer analysis using H3K27ac and H3K4me1, the two replicates were used for the peak calling, and the enhancer was defined as the union of the two peaks for each replication, which required coverage in both replicates and Ͼ50% coverage for the shorter peak.
For the broad domain analysis of H3K4me3, we used MACS (version 2.1.0) to call the broad peaks in the replicates using the options "-fix-bimodal -broad -f BED -g hs -n chip_output -B" (37,38). We filtered out the broad peaks that did not overlap with the regions that were 5 kb up/downstream of the transcription start sites (TSS) of the RefSeq genes. The broad H3K4me3 peaks that were longer than 95% of all of the called H3K4me3 peaks were interpreted as broad H3K4me3 domains. The broad domains were only retained for the subsequent analysis if they were both in the top 5% across the replicates or if one was in the top 5% in one replicate and more than one broad peak was called in the other replicate. The broad H3K4me3 domainassociated genes are defined as the TSS closest to the center of the corresponding broad domain.
For the bivalent H3K4me3/H3K27me3 domain analysis, H3K27me3, H3K4me3 and sequential H3K4me3/H3K27me3 markers were used to identify the actual status of the histone marker-enriched regions. Only the peaks that emerged within the technical replicates of the sequential ChIP of the fetal tissues were taken into consideration, and only the genes within 2 kb up/downstream of the TSS regions that were marked with H3K27me3, H3K4me3, and sequential H3K4me3/H3K27me3 peaks were identified as bivalent H3K4me3/H3K27me3 genes.
Annotation of the Genomic Regions-The exon, intron, and intergenic and intragenic regions were downloaded from the UCSC Genome Browser, and the promoters were defined as the regions Ϯ2 kb from the TSS. The promoters were further classified into the following three classes: high density CpG promoter, intermediate density CpG promoter, and low density CpG promoter, according to a previous publication (39).
Identification of the Enhancers-To identify the enhancer catalog, we took advantage of the H3K27ac and H3K4me1 ChIP-seq peaks, which are the hallmarks of putative enhancers. The H3K27ac-and H3K4me1-enriched regions are considered as distal regulatory element candidates if the peaks are located 2 kb away from the annotated TSS. The distal peaks were merged together if their peak summits were within 1 kb, and the peaks were filtered out based on their read density for the two replicates. To generate the final enhancer catalog, a scanning 2-kb window-based approach with a 100-bp step width was applied as described previously (20). The signal intensity of each window was calculated as the value reads per million bases (RPM). A representative 2-kb window with the strongest signal (highest RPM value) was selected as a regular enhancer, and the consecutive windows were merged together.
Clustering and Correlation Analysis-A non-overlapping 1-kb window was applied to perform the clustering and correlation analysis across the ChIP-seq datasets. Briefly, the genome was first binned into consecutive 1-kb windows, and the signal intensity of each window was calculated based on the normalized raw read counts (RPM). The Pearson correlation coefficient (r) was calculated between the pairwise data using the R package. Also, we used the recently published irreproducible discovery rate method to evaluate reproducibility across replicates. Briefly, we took the advantage of the reproducibility information provided by two "pseudoreplicates" generated by pooling and then randomly partitioning the reads from two replicates, and this information was further used to evaluate the consistency of two independent replicates as reported (see supplemental Table S1) (40).
Identification of the Active and Poised Enhancers and the Enhancer Transition Dynamics-The enhancer activity was defined based on the ratio of the H3K4me1 and H3K27ac enrichment levels, according to a previous publication (20). The distal regions enriched with H3K27ac signals were denoted as potential active enhancers, and these regions always have strong H3K4me1 signals as expected (data not shown). Conversely, the regions with H3K4me1 but without an H3K27ac signal were defined as potential poised enhancers. Strictly, we only use "K27acϩK4me1-positive" or "K4me1-only" DNA segments, instead of active or poised enhancers in most cases hereafter. For the enhancer transitions, we only focused on the K27acϩK4me1-positive DNA segments; the gain/loss of this type of segment across developmental stages is classified as dynamic enhancers. To define the stage-or tissue-specific enhancers, we classified the potential enhancer into active (K27ac-and K4me1-positive) and non-active (K27ac-negative) state. Only the enhancers in one stage are active, while in other developmental stages are non-active and further assigned as stage-specific or tissue-specific enhancers are also defined as these are active in one tissue and non-active in other tissues.
Identification of the Potential Enhancer-targeted Genes within the Topologically Associated Domains-We assumed that active enhancers or super-enhancers are more likely to regulate target genes (not the nearest genes in some cases) within a 100-kb up/downstream region. Enhancers tend to have higher correlations with their target genes than with non-targets; therefore, we systematically computed the correlation matrix between the enhancer signal intensities and expression levels of all of the possible genes within a 100-kb up/downstream region across different stages, and the genes with the highest correlations were selected as potential targeted genes of the given enhancers.
Gene Ontology Analysis of the Genes of Interest-The GO enrichment analysis was performed using DAVID (david.abcc. ncifcrf.gov) (41), and the top 10 non-redundant GO terms were represented with significant p values in a hypergeometric test.
Functional Enrichment and Transcription Factor-binding Motif Analysis-For the given stage-or organ-specific enhancers (or super-enhancers), the functional annotation enrichment analysis was performed using GREAT (42), which is based on the neighboring potential of targeted genes. The GREAT analysis was performed on line using the entire genome as the background and all of the default settings. The MGI phenotypic ontology was extracted and filtered to retain the terms that were significantly enriched. The transcription factor-binding motifs were called using HOMER with the following command: findMotifsGenome.pl input.bed hg19 output_dir -size 2000 -len 8 -S 100 (43).
Identification and Characterization of Super-enhancers-Some typical enhancers tend to cluster together to form superenhancers. The super-enhancers were identified using the ROSE algorithm, as described elsewhere (44,45), and were based on the regular constituent enhancers that have been previously identified. The average ChIP-seq raw read counts of the super-enhancer or regular enhancer regions were integrated as the total signal intensities, and the average signal densities of the super-enhancers or regular enhancers were the raw read counts normalized to the length of the corresponding enhancer region.
Evolutionary Comparison between Human and Mouse-The similar H3K27ac ChIP-seq dataset, including forebrain, heart, and liver, at three different mouse prenatally developmental stages (E11.5, E14.5, and E17.5) was downloaded (GSE52386) and re-analyzed (5). For direct comparison, we converted the mouse mm9 genome coordinates, together with their annotations to the human hg19 genome coordinates using lift-over tool (46), and only the mouse-to-human orthologous regions were used for further analysis. The principal component analysis and unsupervised hierarchical clustering were performed using R packages.  FEBRUARY 26, 2016 • VOLUME 291 • NUMBER 9

Analyses of Active Enhancers in Early Human Embryos-
First, we analyzed H3K27ac and H3K4me1, which are markers for putative active and active/poised enhancers, respectively, in the brain, heart, and liver of early human embryos at 12 weeks of gestation, which is morphologically equivalent to E12.5 to E13.5 in mouse embryos. For each tissue, we performed two technical replicates, and the results were reasonably reproducible, with more than 15,000 peaks in average passing the irreproducible discovery rate threshold (false discovery rate ϭ 1%) and more than 30% reads located in peak regions. Also, we binned the genome into consecutive 1-kb tiles, and find the RPM values were well correlated between replicates, with the average Pearson correlation at 85% (see "Experimental Procedures" and supplemental Table S1). After excluding H3K27ac from the TSS regions, we identified a total of 40,181 K27acϩK4me1-positive regions (potential active enhancers) in these three fetal-stage organs. In comparison, we identified 18,933 and 50,221 K27acϩK4me1-positive regions in hESC and human adult organs, from previously published raw data, respectively. For these K27acϩK4me1-positive regions in the fetal tissues, the majority showed tissue-specific patterns, with 8,557 being restricted only in fetal brain, 7,519 in fetal heart, and 2,757 in fetal liver (see "Experimental Procedures" and Fig.  1A). Moreover, a large number of the enhancers also showed developmental stage-specific patterns, with 6,848, 1,779, and 1,629 being stage-specifically present in the 12-week-old fetal brain, heart, and liver, respectively, but not in the corresponding adult organs (data not shown). Examples of representative enhancer regions are shown in Fig. 1B. To test whether these tissue-or stage-specific enhancers are functionally important, we explored the MGI phenotype database to examine the phenotypes of knock-out mice of their neighboring genes using the default parameters of the on-line GREAT tool (5,42). We found that these fetal stage-specific enhancers showed strong enrichment for the tissue-specific phenotypes of the corresponding organs as expected (data not shown). These findings suggested that the genes associated with these fetal stage-specific enhancers are functionally crucial for the development and physiological function of these tissues.
However, as we know, an enhancer can regulate a gene that is hundreds of kilobases away and sometimes not its nearest gene. To determine the most possible candidate gene regulated by an enhancer, we examined the correlation between the signal intensities of H3K27ac and H3K4me1 and the expression levels of genes within a 100-kb up/downstream region of the  The color key that changes from black to yellow indicates low to high enrichment scores, respectively.
K27acϩK4me1-positive regions in different organs and at different developmental stages. Organ-and stage-specific K27acϩK4me1-positive regions (candidate enhancers) were recognized (criterion see under "Experimental Procedures"), and among them, we identified 1,969 fetal brain-specific, 711 fetal heart-specific, and 811 fetal liver-specific regions ( Fig. 2; supplemental Table S2). The fetal brain-specific candidate enhancers are associated with brain development genes such as PAX6, TUBB2B, NEUROG1, NEUROG2, NEUROD2, and FEZF2. Similarly, the fetal heart-specific candidate enhancers were related to heart development genes such as TNNT2, MEF2D, HAND1, TEAD4, MYH7, and MYOG, and the fetal liver-specific ones are associated with liver-and hematopoiesisassociated functions, such as BLVRA, ALAD, FECH, HMBS, and UROS. Next, we performed a GO analysis and noted that the genes that are possibly associated with these fetal tissuespecific enhancers are clearly involved in the control of the biological processes required for the development and function of the corresponding organ (Fig. 2B). Moreover, by analyzing the sequence motifs of the known transcription factors that were enriched in these fetal-specific enhancers, we found that these enhancers showed strong differential enrichment for transcription factor-binding motifs across tissues and developmental stages (data not shown). These findings indicate that large numbers of the genes regulated by these tissue-and fetal stage-specific enhancers are functionally crucial for the development and physiological function of these fetal tissues.

Analyses of Super-enhancers in Early Human
Embryos-Recently, the super-enhancer has been identified to play a crucial role in cell fate determination during both embryonic development and tumor formation (44,45,47). Therefore, we investigated this type of enhancer region in the early human embryo. We identified 130, 51, and 100 potential super-enhancers in the fetal brain, heart, and liver, respectively (Figs. 3 and 4 and supplemental Table S3). The average sizes of these super-enhancers ranged from 20 to 50 kb, with the H3K27ac signal intensity being 5-8-fold higher than that of the average length of the regular enhancer regions (Fig. 3). The majority of the super-enhancers (82, 44, and 69% for the fetal brain, heart, and liver, respectively) showed tissue-specific and developmental stage-specific patterns, suggesting that super-enhancers are tissue-specific master regulators in early human embryos in      FEBRUARY 26, 2016 • VOLUME 291 • NUMBER 9 JOURNAL OF BIOLOGICAL CHEMISTRY 4391 vivo. For example, a fetal brain-specific 25-kb-long super-enhancer was identified across the SOX1 gene locus, which is a master gene for neurogenesis. A fetal heart-specific super-enhancer was also identified across the HCN4 gene locus, which encodes a member of the hyperpolarization-activated cyclic nucleotide-gated potassium channel family and is necessary for the cardiac pacemaking process. Similarly, a fetal liver-specific super-enhancer was identified across the ALAS2 gene locus, which is essential for definitive erythroblast development (Fig.  4). These findings indicate that these fetal stage super-enhancers potentially regulate a large set of master genes during early embryonic development in vivo.

Epigenomic Landscape of Human Organ Formation
Analyses of Broad H3K4me3 Domains in Early Human Embryos-The broad H3K4me3 domain has recently been proposed in hESCs (38,48). We asked whether this type of regulatory region also exist in the in vivo developing human organs. Using a similar filtering criterion, we identified 204, 645, and 276 broad H3K4me3 domains in the fetal brain, heart, and liver, respectively (supplemental Table S4). Our results showed that the H3K4me3 domains were largely nonoverlapping with the super-enhancers, which is consistent with the finding in hESCs and indicates that these two sets of regions may probably have distinct functions (Fig. 5, A and B). Furthermore, we found that some of the genes marked with the broad H3K4me3 domain were important developmental genes and displayed some disease phenotypes in the MGI database (Fig. 5). For example, the genes associated with the broad H3K4me3 domain in the fetal liver are enriched with phenotypes such as abnormal erythropoiesis and decreased hemoglobin content. The H3K4me3 broad domain-associated genes in fetal liver and brain are involved in some basic cellular metabolic processes, such as cell-cell signal, chromosome organization, and transcription. Taken together, these data suggested that although the genes associated with the broad H3K4me3 domains are different from those associated with super-enhancers, they are likely to have distinct functional properties in different tissues, which reveal the complexity of gene regulation in human early embryos.
Evolutionary Comparison between Human and Mouse-We systematically compared the H3K27ac ChIP-seq signal inten-   sity across the one-to-one orthologous regions in human and mouse. The principal component analysis showed that the enhancer H3K27ac signal was able to distinguish between both the organs and the species (Fig. 6A). It appears that the principal component 1 harbors the most differences between two species (species differences), whereas principal component 2 majorly  FIGURE 6. Enhancer comparisons between human and mouse. A, principal component analysis of the H3K27ac signal intensities between human and mouse in three different organs. B, Pearson correlation clustering analysis of these three organs. The color key from green to red indicates the correlation from low to high. Notably, the mouse genome coordinates (mm9) are converted to the human assembly (hg19) using a lift-over tool, and only the one-to-one orthologous regions were taken into considerations.
separates different organs away (tissue differences), which is consistent with a previous report (49). In the clustering analysis, the human and mouse tissues were clustered separately, and for species, the heart and liver are more closer than the brain (Fig. 6B). Together, the results indicate both the evolutionary conservations of the distal regulatory elements between these two species and the unique features of the human-specific distal regulation.
Analyses of H3K4me3/H3K27me3 Bivalent Domains in Early Human Embryos-The H3K4me3 and H3K27me3 bivalent domain has been well determined in embryonic stem cells (25, 50 -53) However, it is not known whether this type of region also exist in the in vivo developing early human organs. We performed sequential ChIP analyses of H3K27me3 and H3K4me3, and the results showed that 1,461, 1,976, and 780 genes were associated with bivalent H3K4me3/H3K27me3 modifications in their promoter regions in the human fetal brain, heart, and liver samples, respectively (Figs. 7 and 8 and supplemental Table S5). This demonstrated that the bivalent domains are widely present in early human embryos. The gene ontology analysis revealed that the fetal brain-specific bivalent domain-associated genes were strongly enriched for the neural developmental genes such as neuron differentiation (p value ϭ 3.9 ϫ 10 Ϫ29 , hypergeometric test) and pattern specification (p value ϭ 4.4 ϫ 10 Ϫ28 , hypergeometric test). Several representative genes (FGF5, SALL4, and GABRA4) were shown in Fig. 7A, The fetal heart-and liver-specific bivalent domains also occupy developmentally crucial genes, such as MIXL1, HEXB, CACNB4, FOXO1, KITLG, and SLC12A2 (Fig. 7, B and C). The genes associated with the fetal liver bivalent domains also showed enrichment in the enzyme-linked receptor protein signaling pathway, including FGFR3, ERBB3, and BCAR1, which are likely to play roles in response to the extracellular signal to promote cell proliferation or differentiation (data not shown). Collectively, these findings indicate that bivalent H3K4me3/H3K27me3 domains are prevalent in all three fetal organs shortly after organ formation and likely play roles in keeping the progenitor cells in these organs ready for immediate further differentiation into the diverse cell types necessary for development.

Discussion
Epigenetic regulation is important for controlling the precise tissue-specific and developmental stage-specific gene expression patterns in early human embryos. We comprehensively surveyed the epigenomes in 12-week-old human embryos by performing ChIP-seq analyses of H3K27ac, H3K4me1, H3K4me3, and H3K27me3. Moreover, because bivalent H3K4me3/H3K27me3 domains are considered to be important features of stem or progenitor cells, we also analyzed these domains using sequential ChIP-seq. These allow us to obtain a comprehensive landscape of the epigenomes of early human embryos shortly after organ formation. First, we identified more than 40,000 potential enhancers in these fetal tissues, and over 10,000 of them showed fetal-specific patterns, which greatly facilitates our categorization of human embryonic enhancers. Second, we identified a large number of super-enhancers in these fetal tissues, and 91% of them were specifically restricted to this fetal stage, which extends this type of distal regulatory element to early human embryos in vivo (Figs. 3 and  4). Third, we identified a total of 3,433 genes with bivalent H3K4me3/H3K27me3 modifications in these organs (Figs. 7 and 8). We have also made a direct comparison between human and mouse. It appears that the species difference is more prom-  . Representative examples of bivalent H3K4me3/H3K27me3 domains in the human 12-week-old fetal tissues. The ChIP-seq tracks showing the bivalent H3K4me3/H3K27me3 domains in the promoter regions of the representative genes that were analyzed using sequential ChIP-seq (red tracks). The enrichment of the H3K4me3 (green tracks) and H3K27me3 (blue tracks) peaks using regular ChIP-seq in the same loci across two replicates are also shown as controls. Overall, three fetal brain-specific (A), three fetal heart-specific (B), and three fetal liver-specific (C) bivalent H3K4me3/H3K27me3 domains are shown. FEBRUARY 26, 2016 • VOLUME 291 • NUMBER 9

Epigenomic Landscape of Human Organ Formation
inent than tissue differences. However, because the mouse data were downloaded from a previous publication (5), even though we have re-processed the data, this may also reflect a batch effect but not a biological effect. In summary, our work collectively offers some critical clues for dissecting the potential roles of tissue-specific and developmental stage-specific epigenomic features that regulate the spatiotemporal gene expression that occurs during organ formation in human embryos.  /development  abnormal embryonic tissue morphology  abnormal prenatal body size  abnormal neural tube closure  abnormal brain development  complete perinatal lethality  decreased embryo size  open neural tube  abnormal embryo size  abnormal artery morphology  abnormal craniofacial development  abnormal tympanic ring morphology  abnormal middle ear morphology  abnormal developmental patterning  abnormal brain ventricle/choroid plexus morphology  abnormal eye development  abnormal middle ear ossicle morphology  cyanosis  abnormal thymus morphology  abnormal nervous system development  abnormal telencephalon development  abnormal neuron differentiation abnormal olfactory lobe morphology abnormal olfactory bulb morphology abnormal telencephalon morphology abnormal limbic system morphology abnormal neuronal precursor proliferation abnormal hippocampus morphology abnormal temporal lobe morphology abnormal pallium development abnormal cerebrum morphology abnormal neuron number abnormal Ammon gyrus morphology premature neuronal precursor differentiation abnormal brain ventricle/choroid plexus morphology abnormal embryonic/fetal subventricular zone morphology abnormal dentate gyrus morphology skeletal muscle atrophy increased platelet ATP level thin myocardium compact layer muscle fiber atrophy decreased skeletal muscle mass shortened PR interval abnormal circulating creatine level abnormal fetal size aortic dissection decreased erythrocyte cell number abnormal hemoglobin abnormal hematocrit abnormal mean corpuscular volume anemia decreased hematocrit abnormal cell proliferation abnormal megakaryocyte progenitor cell morphology abnormal megakaryocyte morphology enlarged spleen abnormal hemoglobin content abnormal cell cycle increased erythroid progenitor cell number increased mitotic index decreased cell proliferation hemolytic anemia abnormal myelination tremors abnormal glial cell morphology abnormal neurite morphology abnormal myelin sheath morphology ventricle myocardium hypoplasia myocardium hypoplasia trabecula carnea hypoplasia abnormal oligodendrocyte morphology abnormal ventricle myocardium morphology myocardial trabeculae hypoplasia decreased systemic arterial blood pressure abnormal CNS glial cell morphology internal hemorrhage abnormal axon morphology motor neuron degeneration hepatic steatosis abnormal cardiovascular system morphology abnormal cardiovascular system physiology muscle phenotype abnormal heart ventricle morphology abnormal muscle morphology abnormal heart morphology abnormal cardiovascular development abnormal muscle fiber morphology abnormal heart layer morphology abnormal blood circulation abnormal blood vessel morphology abnormal myocardium layer morphology abnormal heart size embryonic lethality during organogenesis abnormal myocardial fiber morphology abnormal cardiac muscle tissue morphology abnormal vascular development increased heart ventricle size complete embryonic lethality during organogenesis cardiac hypertrophy liver/biliary system phenotype abnormal hepatobiliary system morphology abnormal liver morphology abnormal erythrocyte morphology abnormal erythropoiesis increased circulating triglyceride level increased triglyceride level abnormal circulating lipid level abnormal triglyceride level abnormal liver size abnormal lipid homeostasis abnormal hepatobiliary system physiology abnormal liver physiology abnormal lipid level abnormal myeloid leukocyte morphology abnormal cholesterol level decreased hemoglobin content abnormal phagocyte morphology