The accessible chromatin landscape during conversion of human embryonic stem cells to trophoblast by bone morphogenetic protein 4^†

Abstract

Human embryonic stem cells (hESCs) exposed to the growth factor bone morphogenetic protein 4 (BMP4) in the absence of FGF2 have been used as a model to study the development of placental development. However, little is known about the cis-regulatory mechanisms underlying this important process. In this study, we used the public available chromatin accessibility data of hESC H1 cells and BMP4-induced trophoblast (TB) cell lines to identify DNase I hypersensitive sites (DHSs) in the two cell lines, as well as the transcription factor (TF) binding sites within the DHSs. By comparing read profiles in H1 and TB, we identified 17 472 TB-specific DHSs. The TB-specific DHSs are enriched in terms of “blood vessel” and “trophectoderm,” consisting of TF motifs family: Leucine Zipper, Helix-Loop-Helix, GATA, and ETS. To validate differential expression of the TFs binding to these motifs, we analyzed public available RNA-seq and microarray data in the same context. Finally, by integrating the protein–protein interaction data, we constructed a TF network for placenta development and identified top 20 key TFs through centrality analysis in the network. Our results indicate BMP4-induced TB system provided an invaluable model for the study of TB development and highlighted novel candidate genes in placenta development in human.

Introduction

Cis-regulatory regions that regulate spatiotemporal gene expression are crucial to many processes, and often contribute to disease when disrupted [1,2]. The accessible chromatin is functionally related to transcriptional activity, since this remodeled state is necessary for the binding proteins [3]. DNase I hypersensitive sites (DHSs) are regions of chromatin that exposed the DNA and are sensitive to cleavage by the DNase I enzyme for its lost condensed structure. The combination of DNase I digestion and deep sequencing technology (DNase-seq) has been used to reveal chromatin accessibility [4] and underpinned the discovery of all classes of cis-regulatory elements in specific tissue or cell-type on a genome-wide scale [5,6].

The dysfunction of placenta development is highly correlated to the defects of pregnancy and fetal growth. Some of transcription factors (TFs) have been identified as key factors for the placenta development [7–9] in mouse. Yet, the majority of the cis-regulatory where the TFs bind and regulate the placenta development are not identified in human due to the lack of material in vivo [10,11]. Since Xu et al. demonstrated that bone morphogenetic protein 4 (BMP4) could induce human embryonic stem cells (hESCs) to differentiate efficiently to trophoblast (TB) lineage [12], multiple groups have used this system as a model for studying TB lineage specification in vitro. However, recent debate has challenged whether these BMP4-induction model created truly authentic TB [13]. Roberts's lab later demonstrated soundness of this model by analyzing gene expression profiling through RNA-seq technology [14,15]. Whereas, recently developed single-cell technology has revealed biological gene expression noise and the stochastic nature of molecular (RNAs, proteins, or metabolites) network regulations within a homogenous cell population [16–19].

Transcription factors have been identified as key players underling regulatory circuitry to control gene expression [20]. Chromatin profiling shows more precision on the definition of cell type than gene expression profiling [21]. In this work, by comparing the accessible chromatin landscape of the human H1 ESC to BMP4-induced TB which was derived from the ENCODE (encyclopedia of DNA elements) project [22], we identified 17 472 peaks with higher signals in TB than in H1. We then performed functional enrichment, motif scoring analysis, uncovered ontology, and TF motifs enriched within these peaks. We analyzed published transcriptome profiling data to identify the expression level of specific TF family members likely binding to these accessible chromatin domains. Finally, we incorporated protein–protein interaction (PPI) data to construct a TF network, which provides additional evidence for coregulation of genes with cis-regulatory regions, and highlights novel candidate genes that may be important for placenta development in human. A flow chart diagram illustrating the steps in the pipeline is shown in Figure 1. By comparing the accessible landscape of TB to that of H1, we were able to identify DHSs with higher signal at TB. We then performed a motif scoring analysis and discovered that TF motifs enriched within these DHSs. The subset of DHSs containing TF motifs is frequently associated with genes that are important to the invasion process in human.

Figure 1.

Experimental and analytical framework for illustrating the steps of this study. See text for a detailed description.

We analyzed transcriptome profiling to identify the specific TF family members likely binding to these DHSs. Finally, we incorporated PPI data to construct a gene network which highlights novel candidate genes that may be important to TB invasion in human.

Material and methods

Data resources

Processed peak files, RNA-seq FPKM data, and microarray data analyzed in this paper can be found in Table 1. The processed peak files of DNase-seq data were downloaded from website (http://www.genboree.org/epigenomeatlas).

Table 1.

The data sets from different laboratories analyzed in this paper. All of the data can be retrieved through the corresponding GEO number.

Cell or tissue	Species	Platform	GEO number	Reference
H1, BMP4-derived TB	Homo sapiens	Accessible chromatin landscape by DNase-seq	GSE16256	[37]
H1, BMP4-derived TB	Homo sapiens	Expression profiling by RNA-SEQ	GSE16256	[83]
H1, H9, BMP4-induced cells	Homo sapiens	Expression profiling by microarray	GSE10469	[91]
Preimplantation embryo	Homo sapiens	Expression profiling by microarray	GSE22032	[92]
Preimplantation embryo	Mus musculus	Expression profiling by microarray	GSE7309	[93]
TB stem cell (TSC)	Mus musculus	Expression profiling by microarray	GSE12985	[94]

Genomic distribution analysis

General analysis of DNase-seq data was mainly done using tools within the Cistrome pipeline [23]. Genomic location annotation and enrichment profiling were performed with cis-regulatory element annotation system (CEAS) [24] by calculating average read counts to generate average profile plots using the following parameters: Span 3000 bp, Profiling resolution 50 bp, Promoter/downstream lower interval 1000 bp, Promoter/downstream middle-interval 2000 bp, Promoter/downstream upper-interval 3000 bp, Bi-Promoter lower range 2500 bp, Bi-Promoter upper range 5000 bp, and Relative distance 3000 bp.

Identification of accessible chromatin regions

Handling, analysis, and visualization of tracks and scatter plots were done according to [25]. Trophoblast-specific DHSs were identified using peak-finding procedure, with “TB” as the positive sample and “H1” as the negative control using adaptive local thresholds.

Transcription factor binding site prediction and enrichment analysis

SitePro, an aggregation plot tool for signal profiling (version 1.0.0), was utilized to draw the average score profile around the center of H1- and TB-specific DHS. Motif analysis was performed by using the SeqPos Tool (version 1.0.0.) with a P-value cutoff of 0.001 and a scanning width of 600 bp.

By taking the peak locations as the input, SeqPos can find motifs that are enriched close to the peak centers. SeqPos can scan all of the motifs that we collected from JASPAR [26], TRANSFAC [27], Protein Binding Microarray (PBM) [28], and the human protein DNA interaction (hPDI) databases [29]. SeqPos can also find de novo motifs using the MDscan algorithm [30].

Gene expression analysis

Microarray data (Table 1) were processed following previous work [31].

Genomic region enrichment analysis (GREAT) ontology enrichment analysis

Genomic region enrichment analysis was performed using Genomic Regions Enrichment of Annotations Tool (GREAT; http://great.stanford.edu/) v2.0.2 [32], with default parameters, unless otherwise noted in the text. BED file which obtained from peak-finding procedure were take as input for GREAT analysis, which integrates 20 separate ontologies containing biological knowledge about gene functions, phenotype and disease associations, regulatory and metabolic pathways, gene expression data, presence of regulatory motifs to capture cofactor dependencies, and gene families.

Transcription factor network construction

H1- and TB-specific DHS-enriched TFs were used as input for string-db (http://string-db.org/) v9.1 [33]. We required medium–high confidence (STRING score of at least 0.5) for predicted interactions, which were visualized in the confidence view. Interaction enrichment was calculated using the string-db function.

Database for Annotation, Visualization, and Integrated Discovery

The function of TFs identified by SeqPos Tool was investigated by using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) v6.7 [34] based on the gene ontology (GO) annotations [35]. In addition, groups of genes associated with specific pathways, based on the Kyoto Encyclopedia of Genes and Genomes (KEGG), were analyzed together to assess pathway regulation during the reprogramming process.

Centrality analysis

To identify key transcriptional factors, we used centrality analysis. Herein, centrality is a concept which is used to identify the relative importance of a node in a given network. Generally speaking, a centrality measure is a function which assigns a numerical value C(v) to each node of a network, and there are many different definitions and formulas used to measure centrality. It is suggested that the combination of several measures is more effective. We thus engaged degree, eccentricity, closeness, centroid value, stress, shortest path betweenness, Collaborative Filtering (CF)-closeness and PageRank centrality measures in this study. All measures were performed using CentiBiN (http://centibin.ipk-gatersleben.de) [36]. Integrated analysis was ranked by the number of times they appeared in the top 20 key TFs according to each centrality measure.

Results

General features of the accessible chromatin landscape

Genome-wide DNase I sensitivity profiles of H1 human ES cells and BMP4-induced TB was derived from ENCODE project [37]. Density of DNase I cleavage sites for selected cell types provides a continuous quantitative measure of chromatin accessibility, shown for the whole genome region (Figure 2A; Figure S1A). The relative accessibility of DHS along the genome is highly consistent across cell types (Figure 2B; Figure S1B). Approximately 9.1% of DHSs localized at promoters (less than 1000 bp) and 14.6% localized within 3k bp of a promoter. The remaining 85.4% of DHSs localized more distally, and roughly evenly distributed between intronic and intergenic regions (Figure 2C). Near transcriptional start sites (TSSs), chromatin exhibits high accessibility in both H1 and TB (Figure S1D; Figure 2D).

Figure 2.

Characterization of DHS peaks in TB cell lines. (A) Distribution of TB DHS over chromosomes in TB cell lines. The blue bars represent the percentages of mappable regions in the chromosomes (genome background) and the red bars the percentages of the whole regions. These percentages are also marked right next to the bars. P-values for the significance of the relative enrichment of regions with respect to the genome background are shown in parentheses next to the percentages of the red bars. (B) Trophoblast DHS are distributed over the genome along with their scores or peak heights. The line graph on the top left corner illustrates the distribution of peak heights. The red bars plot TB DHS. The maximum peak height is 60 (the largest value on the X-axis of the line). The X-axis of the main plot represents the actual chromosome sizes. (C) The distributed of TB DHS are of over important genomic features. Trophoblast DHS tends to bind to promoters of genes to regulate the genes. (D) Average profiling of TB DHS within/near important genomic features. The panels on the first row display the average enrichment signals around TSS and TTS of genes, respectively. The middle panel (on the second row) represents the average signals on the meta-gene of 3 kb, which shows that TB DHS enriches in promoter and decreases towards the 3' end.

Trophoblast-specific DNase I hypersensitive sites are closed to trophoblast-related genes

By comparing the peaks between H1 and TB, we identified 17 472 TB-specific peaks (Figure 3A) and 64 402 H1-specific peaks (Figure S2A). Ratiometric heat map of the ratio (Figure 3B) and track of superimposed signal (intensity on Y-axis) (Figure 3C) between the “TB” and “H1” signals at the “peaks from TB using H1 as negative control” set of regions confirm the enrichment of DHSs in TB in these regions. Ratiometric heat map of the ratio (Figure S2B) and track of superimposed signal (intensity on Y-axis) (Figure S2C) between the “TB” and “H1” signals at the “peaks from H1 using TB as negative control” set of regions confirm the enrichment of DHSs in H1.

Figure 3.

Identification of TB-specific DHS. (A) Two-dimensional histogram of the “H1 reads” (on X-axis) and “TB reads” (on Y-axis), (B) the “H1 normalized” (on X-axis) and “TB normalized” (on Y-axis) parameters found in the TB-specific DHS. A bar illustrating the relationship between read counts and coloring can be seen on the right side of the plot. The values on the X-axis and Y-axis are presented on a log scale with a base of 10. The counts within each bin are log transformed with a base of 10 before controlling the color. (B) Ratio metric heat map of the ratio between the “TB” and “H1” signals at the “peaks from H1 using TB as negative control” set of regions. Regions were sorted according to decreasing pos reads. A bar showing the relationship between hue and ratio can be seen at the bottom of the plot. The density of the coloring depends on the highest intensity of the two samples at each position. (C) Track of the “H1” and “TB” superimposed signal (intensity on Y-axis) at the TB-specific DHS. Numbers correspond to chromosomal coordinates. Horizontal bar illustrates the size of the genomic region.

With the GREAT [32], the peaks were converted to mouse genomic coordinates. We found that the TB-specific DHSs were enriched in the placenta development-associated GO terms including morphological processes, mouse phenotypes, and Mouse Genome Informatics (MGI) expression data. In GO Biological Process ontology, the enriched terms are “regulation of cell aging” and “regulation of cellular senescence.” These terms correlated to placental health and pathology associated with cellular senescence [38]. For MGI, which provide integrated genetic, genomic, and biological data to facilitate the study of human health and disease [39], the enrichment is related to “blood vessel,” “amnion,” “trophectoderm” terms (Figure 4A). Table 2 lists genes related to the each term identified above. Our analysis truly contains errors and many terms do not have enough terms associated with. For example, genes such as ARID5B, CCND, DSP, FRS2, KAZN, etc. of the genes linked to TB are certainly not classical TB markers. This is because GO uses the true path rule. If a gene is annotated by a term, it is also implicitly annotated by all the parents of this term, up to the root [40]. Future work based on collaborations between the placenta biology community and members of the GO Consortium would lead to an increase in the number and specificity of TB- and placenta-related GO terms [41].

Figure 4.

Functional enrichment analysis from TB-specific DHS from GREAT. (A) Terms in GREAT for GO biological process and MGI mouse phenotype ontologies are shown. (B) Track of the “H1” and “TB” signal (intensity on Y-axis) at HAND1 and FOXO1. (C) Chromatin access at both megabase and fine scale.

Table 2.

Genes and related GO terms enriched in peaks from TB using H1 as negative control by GREAT.

Term	Related genes
Regulation of cellular senescence	CDK6, CDKN2A, HMGA1, HMGA2, ING2, NEK6, NUAK1, PNPT1, SIRT1, TWIST1, VASH1, ZKSCAN3
Trophectoderm	ARID5B, CCND3, DSP, FRS2, GATA2, HAND1, IGF2, KAZN, NEK6, SNAI1, TFAP2C
Blood vessel	FOXO1, GCOM1, JUN, KDR, MYZAP, SGK1, SOX17, TGM2

We observed the genes that are known to be important for placental development. For example, Heart and Neural Crest Derivatives Expressed 1 (HAND1) [42] (Figure 4B) and its enhancer exhibited higher accessibility in TB than in H1. The promoter of TF FOXO1, which has been reported essential for vascular growth [43] (Figure 4B), shows higher accessibility in TB than in H1. Figure 4C shows per-nucleotide DNase I cleavage density across a 50-kb region of chromosome 5 (positions 153860000–153905000) containing open reading frames as enhancers of HAND1. Our analysis reveals multiple putative enhancers in the regulatory domain of Hand1. Magnification shows positions of individual DNase I cleavage events, revealing DNase I footprints with known motifs for TB related to regulatory factors TFAP2C, GATA2, GATA3, and TEAD4 (Figure 4C).

We next employed microarray data published previously (Table 1) to analyze expression level of HAND1 in human ES to TB system and mouse trophoblast stem cell (TS) to trophoblast giant cells (TGC) system. Hand1 shows an unregulated tendency in both systems (Figure S3A; Figure S4). We also listed genes which may also be useful for scientific community. For example, placenta-specific 8 (PLAC8), which was recently identified and expressed in human oocytes and preimplantation embryos [44], shows an upregulation in both human ES to TB system and mouse TS to TGC system according to microarray data (Figure S3A and S4). We also identified the potential enhancer that promotes the expression level of PLAC8 in TB compared with H1 (Figure S3B).

Motif analysis reveals transcriptional code in trophoblast-specific DNase I hypersensitive sites

In order to determine whether a subset of TB-specific DHSs are coregulated by a common set of TFs, we searched SitePro for TF motifs that occur significantly more frequently in TB-specific DHS than in either of two background sets. We scanned a library of 1651 nonredundant motifs, requiring that predictions be conserved in mouse and human. We identified 104 TFs that are enriched at the center of TB-specific DHSs (Table S1). The TFs that are corresponded to the enriched motifs include Leucine Zipper, Homeodomain, Helix-Loop-Helix, GATA, DNA Polymerase-Beta, and Ets Domain family (Table S1). The most conserved binding sites are for those motifs recognized by JUND, JUNB, FOSL1, and JUN. Figure 5A shows the represented DNA motifs related to genes that are enriched at the center of TB-specific DHS motifs. We also noticed that the lost DHSs significantly enrich the motifs of key pluripotency factors (OCT4, SOX2, NANOG) during the transition from ESCs to TB [45] (Table S2). These results confirm the validation of our analysis. We also employed SeqPos to find de novo motifs using MDscan algorithm [46] and do screen based on curated cistrome motif database. This database includes human and mouse data and puts data from Transfac, JASPAR, UniPROBE (pbm), hPDI together, and it includes the motifs derived from ChIP-seq data. This analysis identified 41 TFs for H1-specific DHSs and 154 TFs for TB-specific DHSs (Table S4).

Figure 5.

Motif scoring and interactions in TF-enriched TB-specific DHS. (A) Logos of represented DNA motifs related to genes that are enriched to the center of TB-specific DHS motifs scan all of the motifs using SeqPos. (B) Gene expression level of TFs enriched for TB-specific DHS in H1 and four hESC-derived lineages [83], including trophoblast-like cells (TBL), mesendoderm (ME), neural progenitor cells (NPCs), and mesenchymal stem cells (MSCs). (C) The interaction network of TFs enriched for TB-specific DHS. The network consists of 641 interactions among 107 TFs.

KEGG pathway analysis of the 104 TFs indicates an enrichment of several TB differentiation related to pathway identified previously (Table 3), such as “Hippo signaling pathway” [47], “Osteoclast differentiation” [48], “Estrogen signaling pathway” [49], “MAPK signaling pathway” [50], “PPAR signaling pathway” [51], “TGF-beta signaling pathway” [52], “Adipocytokine signaling pathway” [53]. More detailed pathway information could be found in Table S5.

Table 3.

Signaling pathways enriched in 104 TFs by the KEGG pathway analysis (P < 0.05).

Term	P-value
Transcriptional misregulation in cancer	3.84E-16
HTLV-I infection	2.36E-07
Osteoclast differentiation	0.001159
Hepatitis B	0.001904
Hippo signaling pathway	0.002493
Estrogen signaling pathway	0.003051
TNF signaling pathway	0.003933
MAPK signaling pathway	0.004584
Inflammatory bowel disease	0.006322
Colorectal cancer	0.00662
Amphetamine addiction	0.007563
PPAR signaling pathway	0.010454
Pathways in cancer	0.010814
TGF-beta signaling pathway	0.014378
Thyroid cancer	0.015269
Cocaine addiction	0.036858

RNA-Seq data validate differential expression of predicted transcription factors

In order to determine TFs most likely binding to TB-specific DHSs, we focused on the expression levels for each member of these TF families. TFs highly expressed in TB compared with H1 have high potential to bind to TB-specific DHSs (Figure 5B). These expression-based predictions are promising, TFs BACH2, ESRRG [54], FOSB [55], GATA2 [56], GATA3 [56], MAF [57], PPARG [58], TFAP2A [59] have been implicated in placental development. We also analyzed the chromatin access landscape in selected ES and TB-specific TFs (Figure S5). ES-specific TFs (POU5F1, SOX2, NANOG) show open access genome in H1 compared with TB and TB-specific TFs (ETS1, PPARG), and H19 [60] shows open access genome in TB compared with H1 (Figure S5). Our analysis suggests that TF-related TB-specific DHSs are working together to regulate placenta development. Notably, BACH2 has been revealed as a “guardian” TF to regulate super enhancer-associated cytokines and cytokine receptors in CD4(+) T cells [61]. Future study could concentrate the function of this TF in human placenta development process.

Construction and centrality analysis for transcriptional regulatory network

We then used the TFs as input to the STRING database (Search Tool for the Retrieval of Interacting Genes/Proteins) [62], a database for building functional associations between proteins based on PPI information that can already be retrieved from a number of online resources. The functional network of TFs in H1 showed a central, highly interconnected area in which common pluripotency regulators such as POU5F1, NANOG, SOX2, and CCCTC binding factor (CTCF) were identified (Figure S6). The string-db network highlights a core set of interactions, which consists of multiple genes known to be important in placenta development such as GATA2, PPARG, TEAD4. Thus, the genes are highly enriched for interactions and are potentially functionally important (Figure 5C). After removing isolated nodes and separated links, we obtained a network with 641 interactions among 107 TFs (Figure S7). This protein interaction network provides a model for exploring neo-factors that play essential role in the process of placenta development. The structure of the PPI network is related to whether or not a given protein is essential [63]. We use CentiBiN to identify crucial elements of biological networks with multicentrality measures. Firstly, we ranked the top 20 key TFs by every centrality measure, and because different measures focus on different aspects of centrality, so the key TFs might be with different rank according to different centrality measures (see Table S2). Then, to select the top 20 key TFs from an integrated view, we ranked all TFs in Table 4 by the number of times they appeared there. Our analysis highlights a core set of TFs, which consists of multiple genes widely known to be important in placenta development, such as GATA3 [64], PPARG [58], ETS1 [65]. Furthermore, CREB1 recently have been identified critical for the regulation of furin expression during human TB syncytialization [66]. A total of 15/20 (75%) of the TFs have been reported (Table 4). Our analysis provides other potential TFs that play an important role during human TB differentiation process.

Table 4.

The top 20 key TFs corresponding to integrated analysis.

Rank	TF	Description	Number	Reference
1	CREB1	Cyclic AMP-responsive element-binding protein 1	9	[66]
2	ESR1	Estrogen receptor	9	[95]
3	FOS	Proto-oncogene c-Fos	9	[55,71]
4	JUN	TF AP-1	9	[71]
5	JUND	TF jun-D	9	[96]
6	MYOD1	Myoblast determination protein 1	9	[97]
7	MYOG	Myogenin	9	[97]
8	SMAD3	Mothers against decapentaplegic homolog	9
9	GATA1	Erythroid TF	8	[98]
10	GATA3	GATA binding protein 3	8	[64]
11	JUNB	TF jun-B	8	[99]
12	PPARA	Peroxisome proliferator-activated receptor alfa	8	[100]
13	PPARG	Peroxisome proliferator-activated receptor gamma	8	[58]
14	RUNX1	Runt-related TF 1	8	[101]
15	SMAD4	Mothers against decapentaplegic homolog	7
16	TCF3	TF E2-alpha	7
17	ETS1	Protein C-ets-1	6	[65]
18	GATA2	GATA-binding factor 2	6	[64]
19	ATF2	Cyclic AMP-dependent TF ATF-2	5
20	PAX2	Paired box gene 2	5

Discussion

Transcriptional regulation was first identified as a pivotal process that confers cellular identity and modulates the biological activities within a cell [67]. Cis-regulatory alterations accompanying different growth conditions or cell differentiation and cycling impact multiple regulators simultaneously and are difficult to understand.

Footprint data derived in vivo can address all factors simultaneously in their native state and detect regions of direct binding with nucleotide precision [68,69]. Our analysis shows that this technology could uncover key regulation sequences important to placenta development in human. By comparing the accessible chromatin landscape at HB to that at H1, we identified 17 472 TB-specific DHSs. Many of these DHSs are at the upstream of the genes that are associated with terms “blood vessel” and “trophectoderm.” By performing motif scanning analyses, we brought thousands of peaks down to 104 TFs. This set of DHSs is most significantly associated with the genes known to be important in processes of placenta development.

In order to identify the most likely TF family member binding to this set of DHSs, we analyzed public available RNA-seq and microarray data, and determined that BACH2, ESRRG, FOSB, GATA2, GATA3, GRHL1, PPARG, RUNX1, and TFAP2A are strong candidates for coregulating the placenta development genes. Notably, these TFs did not show upregulation in mesendoderm derived from BMP4+FGF2-treated hESCs [70].

We constructed and performed centrality analysis of TF network. AP-1 TFs FOS and JUN, which were previously identified as critical factor for cell invasion in EVT [71], were at the center of the network. This indicates that the network we constructed opens up multiple avenues of research that will help us understand in more detail the mechanisms regulating placenta development. Our analysis also indicates that BMP4-treated hESCs in the absence of FGF2 are unlikely to give rise to solely mesoderm derivatives as previously considered [72]. Hence, hESCs treated with BMP4 provide an invaluable model for studying transition to the EVT lineage.

Our study has a number of important limitations. For example, millions of regulatory sequences have been predicted in the human genome through epigenetic analysis, but few have been shown to regulate transcription in their native contexts [73]. Hence, a major drawback of this study is a lack of functional characterization the large number of predicted cis-regulatory elements with regard to their contributions to target gene expression. Furthermore, our analysis miss some important TFs. Transcription factors whose DNA-binding domains (DBDs) are highly similar were found to bind to motifs that were almost indistinguishable from one another [74]. This could explain that our analysis identifies ETS1, but not ETS2, which was considered key TF during the process of TB differentiation [75–79], for the reason that ETS1 and ETS2 have very similar DBDs [80]. We used JASPAR [26], TRANSFAC [27], PBM [28], and the hPDI databases [29] to characterize the DNA-binding specificities of TFs. Nevertheless, only the DNA-binding properties of less than half of all human TFs have been experimentally defined [81]. This incomplete TF motif catalog limits the ability to analyze the effects of genetic variation on TF-DNA binding, as a large fraction of human TFs cannot be taken into account in these studies. Recently developed technology such as selective microfluidics-based ligand enrichment followed by sequencing (SMiLE-seq) [82], which yields independent and quantitative DNA-binding data, would facilitate achieving a better understanding of genome function during the development of human placenta. Finally, our study only analyzes complex and heterogeneous mixture of cell types derived from BMP4-treated human ES cells [15,83]. Moreover, BMP4 in combination with two small compounds, the activin A signaling inhibitor A83-01 and the FGF2 signaling inhibitor PD173074, makes the TB-induced system more efficient and synchronous [15]. We expect that future epigenetic profile study based on this system combined single-cell technology [84] and would add more insight into TB development in human.

In this study, we annotate gene regulatory elements and investigate their function during the process of throphoblast development by RNA-seq and DNase-seq technology. With the development of high-throughput sequencing technology, biological big data analytics and integration strategies will play a fundamental role to better define cell types or states and the differences among cell types/states during the process of TB development. Liu et al. achieved epigenetic edit by fusion of Tet1 or Dnmt3a with a catalytically inactive Cas9 (dCas9) and altered gene expression in the neighboring loop by targeted de novo 415 methylation of a CTCF loop anchor site [86]. Further experimental work that combines recently developed clustered regularly interspaced short palindromic repeats (CRISPR) technology would facilitate the understanding on which inferred binding sites are functional and which genes they regulate [85].

In sum, transgenic experiments [87], high-throughput reporter assays [88,89], gene-knockout studies in mice [10], and more recently, CRISPR technology [90] in cultured cells could be employed to test predicted cis-regulatory elements activities in the future.

Supplementary data

Supplementary data are available at BIOLRE online.

Supplementary Figure S1. Characterization of DHS peaks in H1 cell lines.

Supplementary Figure S2. Identify peaks—regions with an enrichment of signal in H1. (A) Two-dimensional histogram of the “TB reads” (on X-axis) and “H1 reads” (on Y-axis), the “TB normalized” (on X-axis) and “H1 normalized” (on Y-axis) parameters found in the H1-specific DHS. A bar illustrating the relationship between count and coloring can be seen on the right side of the plot. The values on the X-axis and Y-axis are presented on a log scale with a base of 10. The counts within each bin are log transformed with a base of 10 before controlling the color. (B) Ratio metric heat map of the ratio between the “TB” and “H1” signals at the H1-specific DHS. A pseudo count of 1 was added to both signals. Regions were sorted according to decreasing pos reads. A bar showing the relationship between hue and ratio can be seen at the bottom of the plot. The density of the coloring depends on the highest intensity of the two samples at each position. (C) Track of the “H1” and “TB” superimposed signal (intensity on Y-axis) at 64402 regions from H1-specific DHS. Numbers correspond to chromosomal coordinates. Horizontal bar illustrates the size of the genomic region.

Supplementary Figure S3. Gene expression level of pluripotent-related genes and trophoblast genes in human H1, H9 ESC, BMP4-derived trophoblast [91], and preimplantation embryo [92] by microarray analysis.

Supplementary Figure S4. Gene expression level of pluripotent-related genes and trophoblast genes in mouse trophoblast stem cells (TSC), trophoblast giant cells (TGC) [94], and preimplantation embryo [93] by microarray analysis.

Supplementary Figure S5. Track of the “H1” and “TB” signal (intensity on Y-axis) at a single region corresponding to specific position. Numbers correspond to chromosomal coordinates. Horizontal bar illustrates the size of the genomic region.

Supplementary Figure S6. The interaction network of TFs enriched for H1-specific DHS.

Supplementary Figure S7. String high-resolution bitmap of the interaction network of TFs enriched for TB-specific DHS.

Acknowledgments

The authors would like to thank Jinjing Yang for critical reading and helpful discussion for the manuscript.

Author notes

†

Grant Support: YJL is supported by the postdoctoral scientific research start-up fund of Southeast University of China (1107030155) and National Basic Research Program of China (2012CB316500). HDL is supported by the National Natural Science Foundation of China (No. 31371339). XS is supported by Key Research and Development Program of Jiangsu province (BE2016002-3) and National Basic Research Program of China (2012CB316500).