FOXD3 Regulates Pluripotent Stem Cell Potential by Simultaneously Initiating and Repressing Enhancer Activity

Abstract

Early development is governed by the ability of pluripotent cells to retain the full range of developmental potential and respond accurately to developmental cues. This property is achieved in large part by the temporal and contextual regulation of gene expression by enhancers. Here, we evaluated regulation of enhancer activity during differentiation of embryonic stem to epiblast cells and uncovered the forkhead transcription factor FOXD3 as a major regulator of the developmental potential of both pluripotent states. FOXD3 bound to distinct sites in the two cell types priming enhancers through a dual-functional mechanism. It recruited the SWI/SNF chromatin remodeling complex ATPase BRG1 to promote nucleosome removal while concurrently inhibiting maximal activation of the same enhancers by recruiting histone deacetylases1/2. Thus, FOXD3 prepares cognate genes for future maximal expression by establishing and simultaneously repressing enhancer activity. Through switching of target sites, FOXD3 modulates the developmental potential of pluripotent cells as they differentiate.

Figure 1. Profiling of Changing Enhancers Identifies FOXD3 as a Key Player during the Transition from the ESC to EpiC State

(A) Schematic of transitions followed during early differentiation from ESCs (R) to EpiCs (dY). Colors represent expression from miR-290 locus (Red or R, mCherry), miR-302 locus (Green or G, eGFP), both (Yellow or Y), or neither (Black or B). The letter “d” prior to color represents cells undergoing differentiation from ESCs (R) to dR (16 hr off LIF and 2i) to EpiCs (dY) (52 hr off LIF and 2i). Subsequent differentiation yields non-homogenous transition of cells toward dG (green) and dB (black, no reporters), which are FACS sorted for purity. G is a derived self-renewing EpiSC state grown in FGF and Activin. It requires multiple passages in FGF and Activin prior to establishment, represented as a broken line. (B) Principal component analysis of microarray data of the populations in (D). Array data can be found at GEO: GSE54341 (Parchem et al., 2014). (C) Heatmap of six states of chromatin determined by ChromHMM using H3K4me1, H3K27ac, H3K27me3, and P300 ChIP-seq data in R, dR, dY, and G cells. States 1–3 represent the majority of the genome and show a general lack of these markers (percentages of genome in states 1 and 2 are shown in pie chart in Figure S4). State 3 shows increased H3K27me3 consistent with heterochromatin. State 4 similarly shows increased H3K27me3, but also detectable levels of H3K4me1. States 5 and 6 show low H3K27me3 and increased H3K4me1. State 5 shows low levels of H3K27ac, while state 6 shows elevated H3K27ac. Intensity of blue indicates average ChIP-seq signal for each modification in each state. (D) Genome-wide transitions from states 5 and 6 to all states 1–6. x axis shows the two starting states 5 and 6, and the bar graphs represent the distribution of the end states (1–6) for the matching regions following transition between ESC and EpiC. (E) Metagene analysis of H3K27ac levels centered around enhancers that transition from state 5 (“primed”) to state 6 (“active”) (left), or vice versa (right). Enhancers were defined by the identification of nucleosome-depleted regions in state 6 chromatin as determined by an adapted algorithm (Wamstad et al., 2012) (see Experimental Procedures). Metagene analysis for a 750 bp region surrounding the defined enhancer is shown, with a moving window average of 500 bp windows and 100 bp steps. (F) RSAT (regulatory sequence analysis tool) was used to identify enriched motifs within nucleosome-depleted regions of identified enhancers that transition between primed and active during ESC to EpiC differentiation (Thomas-Chollier et al., 2011). Shown is the most commonly identified motif found in all transitions between poised and active enhancers, which significantly overlaps with the JASPAR motif for FOXD3 (overlap underlined in red). Significance and enrichment of overlap with FOXD3 and other uncovered TF binding are shown in Table S2.

Figure 2. FOXD3 Is Mobilized during the Transition from the ESC to EpiC State, and Its Departure Is Associated with Increased Gene Expression

(A) Scatterplot of reads per million (RPM) at the union of all peaks identified by MACS (Zhang et al., 2008) in both ESCs and EpiCs. Peaks were re-characterized as ESC-specific (red), EpiC-specific (yellow), and common peaks (black) based on their distance from the x = y diagonal (threshold set at ± 0.1) to obtain high-confidence specific peaks (q < 0.01). ESC only peaks = 5,350. EpiC only peaks = 6,737. Common peaks = 1,650. Samples were prepared in triplicate. (B) mRNA expression (from Illumina microarray data) for genes closest to FOXD3 peaks in ESC (red line) and EpiC (yellow line), as well as all other genes on the array (blue line). Microarray experiments were performed in R (ESC), dR, dY (EpiC), dG, and dB cells. dG and dB (differentiated Black) represent further sequential differentiation following EpiCs during continued culture in the absence of LIF and 2i. For this experiment, all populations were collected by FACS. Data shown are all relative to all genes in the ESC state. (C) Fraction of genes with increased expression from ESCs (R) to EpiCs (dY), as well as from EpiCs (dY) to dG, for genes near FOXD3 peaks compared with all other genes. Genes associated with ESC FOXD3 peaks increase with transition to EpiC but then level off. In contrast, genes associated with EpiC FOXD3 peaks continue to go up after EpiC. Asterisk (*) denotes Z scores > 2, p < 0.05 relative to genes without FOXD3 peaks. (D) Heatmap showing hierarchical clustering (single linkage) of all genes shown in (F) and (G) for the transition from ESC (R) to EpiC (dY) to dB. Data are shown as a log2 transformed fold change relative to ESCs. (E) Top 12 GO terms identified by gene ontology analysis using DAVID for genes near FOXD3 peaks in ESCs. (F) Same as in (D) for genes near FOXD3 peaks in EpiCs.

Figure 3. FOXD3 Precedes Binding of Other TFs, Such as OCT4, at Enhancers

(A) Motifs enriched within a 200 bp window surrounding FOXD3 peaks in ESCs cells, shown as a percentage of peaks containing the motif. Shades of gray represent p values as shown (full data in Table S2). (B) Same as (A) but for FOXD3 peaks in EpiCs (full data in Table S2). (C) Average expression of TFs whose motifs were found to be neighboring FOXD3 in either ESCs or EpiCs. For motifs with multiple potential binding TFs, an average of all family members was used. The red line is for TFs to motifs associated with FOXD3 sites in ESCs and the yellow line is for TFs to motifs associated with Foxd3 sites in EpiCs. (D) Scatterplot of RPM in ESCs versus EpiCs at all OCT4 peaks. ESC-specific OCT4 peaks are in red, EpiC-specific peaks are in yellow, and common peaks are in black. Similar to FOXD3, OCT4 shows very distinct binding sites in ESCs and EpiCs, even though sites are distinct from FOXD3. Peaks were identified similarly to FOXD3 peaks in Figure 3A. (E) Scatterplots of RPM for peaks of OCT4 (y axis) versus RPM for peaks of FOXD3 (x axis) in ESCs. ESC-specific FOXD3 peaks are in red, EpiC-specific peaks are in yellow, and common peaks are in black. (F) Same as (E), but in EpiCs. (G) Metagene analysis summarizing OCT4 and FOXD3 binding for regions surrounding FOXD3 peaks in ESCs. Shown are the data for ESCs in a 4,000 bp region surrounding the FOXD3 peaks using a moving window average of 500 bp at 100 bp steps. (H) Same as (G), but in EpiCs. (I) Magnification of OCT4 binding data shown in (G) along with heatmap showing binding of OCT4 in dR and EpiCs, at sites bound by FOXD3 in ESCs.

Figure 4. FOXD3 Binds Nucleosome Occupied Sites and Establishes Primed Enhancers

(A) Metagene analysis of H3K4me1 MNase-ChIP-seq data for ESCs (red line) and EpiCs (yellow line) in a 4,000 bp window surrounding FOXD3-bound sites in ESCs. (B) Same as (A) but for EpiC-bound sites. Asterisk (*) denotes paired Student's t test, p < 0.05. (C) Metagene analysis of H3K4me1 MNase-ChIP-seq data from WT and Foxd3 KO ESCs using a 4,000 bp window surrounding FOXD3-bound sites in ESCs. Data are shown for WT (black line) and Foxd3 KO (green line) cells. (D) Same as (C) but in EpiCs. In addition, data are shown for KO of Foxd3 16 hr prior to differentiation (orange line; KO-l stands for KO long). (E) Metagene analysis of H3K4me1 MNase-ChIP-seq data from WT and Foxd3 KO ESCs using a 4,000 bp window surrounding FOXD3-bound sites in EpiCs. Data are shown for WT (black line) and Foxd3 KO (green line) cells. (F) Same as (E) but in EpiCs. In addition, data are shown for KO of Foxd3 16 hr prior to differentiation (orange line; eKO stands for early KO). (G) MNase tiling qPCR at an ESC enhancer (see Supplemental Experimental Procedures for primers) in WT ESCs, WT EpiCs, and early Foxd3 KO EpiCs. The graph represents an average of n = 2–3 experiments for each primer pair. Unpaired Student's t tests were conducted for WT and Foxd3 KO EpiCs, and the heatmap for the p values for this test is shown below the graph. ANOVA for a three-way comparison of WT ESCs, WT EpiCs, and Foxd3 KO EpiCs was also performed (p = 0.97964, not shown on graph). (H) Same as (F) but for an EpiC enhancer, with t test p values shown as a heatmap. ANOVA p = 1.111e–07, not shown on graph.

Figure 5. FOXD3 Maintains Lower Levels of Histone Acetylation at Enhancers

(A) Metagene analysis of H3K27ac ChIP-seq data from WT and Foxd3 KO ESCs using a 4,000 bp window surrounding FOXD3-bound sites in EpiCs. Data are shown for WT (black line) and Foxd3 KO (green line) cells. Asterisk (*) denotes paired Student's t test, p < 0.05. (B) Data in (A) shown as heatmaps. Brackets represent genes with significant peaks of H3K27ac (signed-rank test, p < 0.05). (C) Same as (A) but for H3K27ac ChIP-seq at EpiC sites. Asterisk (*) denotes paired Student's t test, p < 0.05. (D) Data in (C) shown as heatmaps. Brackets represent genes with significant peaks of H3K27ac (signed-rank test, p < 0.05). (E) Scatterplot showing average gene expression on the x axis and fold change between WT and Foxd3 KO cells on the y axis. Data are shown on the left for ESCs and the right for EpiCs. Points in red are significantly upregulated (adj. p < 0.05) and points in green are significantly downregulated (adj. p < 0.05). (F) Venn diagram showing overlap of upregulated and downregulated genes upon Foxd3 KO in ESCs and EpiCs (adj. p < 0.05). (G) Top ten GO terms identified by gene ontology analysis using DAVID for genes affected by Foxd3 KO in ESCs (top) and EpiCs (bottom). (H) Microarray expression data for WT and Foxd3 KO ESCs. Plotted are all genes that are significantly changed upon Foxd3 KO (“all”) versus genes near FOXD3 binding sites. (Asterisk denotes p < 0.05 and Z > 2 in comparison between two populations). All n = 2142 and FOXD3 n = 511. Inset: expression of FOXD3-bound and Foxd3 KO-affected genes (n = 511) during differentiation. (I) Same as (G), but for EpiCs. All n = 1,800 and FOXD3 n = 412. Inset: expression of FOXD3-bound and Foxd3 KO-affected genes (n = 412) during differentiation. ChIP-seq samples were performed in duplicate.

Figure 6. FOXD3 Recruits the Chromatin Remodeling Swi/Snf Complex and the Histone Deacetylases HDAC1/2 to Target Enhancers

(A) Immunoprecipitation using M2 anti-FLAG resin (FLAG-IP) on extract from WT (C, Control) and FOXD3-3X-FLAG targeted (F, 3XFLAG FOXD3) cell lines (1, 2). Western blots for FLAG and BRG1 are shown. Inset: 10% of input. Representative blot for n = 3 is shown. (B) Overlap of mRNA expression changes in Brg1 KO and Foxd3 KO relative to WT mouse ESCs. The proportion of genes affected by Brg1 KO is shown as a total of all genes, genes altered by Foxd3 loss, and genes both bound by FOXD3 and altered by Foxd3 loss. The latter two categories are significantly enriched for genes whose expression is changed by Brg1 deletion. Results of a Chi-square test are shown. *χ² = 125, p < 2.2e–16. **χ² = 346, p < 2.2e–16. (C) Metagene analysis (left) of BRG1 ChIP-seq data from WT and Foxd3 KO ESCs using a 4,000 bp window surrounding FOXD3-bound sites in ESCs. WT (black line) and Foxd3 KO (green line) cells are shown. Data are visualized in heatmaps on the right. (D) Same as (C), but in EpiCs. (E) Same as (C), but for FOXD3 EpiC bound sites. (F) Same as (D), but for FOXD3 EpiC bound sites. Asterisk (*) denotes paired Student's t test, p < 0.05. (G) Immunoprecipitation using M2 anti-FLAG resin (F) on extracts from WT (Control) and FOXD3-3X-FLAG targeted (3XFLAG FOXD3) cell lines. Western blots for FLAG, HDAC1, HDAC2, and HDAC3 are shown. Inset: 10% input. Representative blots for n = 3 are shown. (H) Overlap of mRNA expression changes in HDAC1/2 and Foxd3 KO relative to WT mouse ESCs. The proportion of genes affected by HDAC1/2 KO is shown as a total of all genes, genes altered by Foxd3 loss, and genes both bound by FOXD3 and altered by Foxd3 loss. The latter two categories are significantly enriched for genes changed by HDAC1/2 deletion. Results of a Chi-square test are shown. *χ² = 49, p = 2.5e–12. **χ² = 127, p < 2.2e–16. (I) Metagene analysis (left) of HDAC1 ChIP-seq data from WT and Foxd3 KO ESCs using a 4,000 bp window surrounding FOXD3-bound sites in ESCs. Data are shown for WT (black line) and Foxd3 KO (green line) cells. Asterisk (*) denotes paired Student's t test, p < 0.05. Data are visualized in heatmaps on the right. Top bracket represents genes that have significant peaks of HDAC1 in ESCs that are lost in EpiCs (signed-rank test, p < 0.05), and bottom bracket represents genes that have significant peaks of HDAC1 in ESCs that are maintained in EpiCs (signed-rank test, p < 0.05). (J) Same as (I) but in EpiCs (same brackets on heatmaps as C). Expr, log2 fold change in expression of associated genes from ESC to EpiC, showing a greater fold increase in expression of genes that lose HDAC1. The groups are significantly different, as per an unpaired Student's t test (p < 0.05). (K) Same as (I) but for FOXD3 EpiC bound sites. (L) Same as (I) but for FOXD3 EpiC bound sites. Asterisk (*) denotes paired Student's t test, p < 0.05. ChIP-seq samples were performed in duplicate.

Figure 7. FOXD3 Recruits HDAC1/2 and the SWI/SNF Complex to Sites to Both Initiate and Attenuate Enhancers

(A) Top: immunoprecipitation using IgG, anti-HDAC1, or anti-BRG1 antibodies on extracts from ESCs. Western blots for HDAC1 and BRG1 are shown. Bottom: same as top, but with western blot for FLAG (i.e., FOXD3). Inset: 10% input. Representative blot for n = 3 is shown. (B) Sequential ChIP for HDAC1 and BRG1 at four ESC FOXD3-bound sites and two EpiC FOXD3-bound sites in ESCs. Shown are all combinations of sequential ChIP with IgG, HDAC1, and BRG1 antibodies. In all cases, the results are expressed relative to IgG in second IP control. Background is shown as IgG used in first IP, which is used for all significance calculations. Error bars = SD (n = 3–4). Significance calculations were performed using a pairwise t test. p < 0.05. (C) Model of FOXD3 function at developmental enhancers. FOXD3 primes enhancers for future gene expression. Upon FOXD3 departure, these genes can be directly activated or repressed by alternative TFs.