Global Genome Conformational Programming during Neuronal Development Is Associated with CTCF and Nuclear FGFR1-The Genome Archipelago Model

Abstract

During the development of mouse embryonic stem cells (ESC) to neuronal committed cells (NCC), coordinated changes in the expression of 2851 genes take place, mediated by the nuclear form of FGFR1. In this paper, widespread differences are demonstrated in the ESC and NCC inter- and intra-chromosomal interactions, chromatin looping, the formation of CTCF- and nFGFR1-linked Topologically Associating Domains (TADs) on a genome-wide scale and in exemplary HoxA-D loci. The analysis centered on HoxA cluster shows that blocking FGFR1 disrupts the loop formation. FGFR1 binding and genome locales are predictive of the genome interactions; likewise, chromatin interactions along with nFGFR1 binding are predictive of the genome function and correlate with genome regulatory attributes and gene expression. This study advances a topologically integrated genome archipelago model that undergoes structural transformations through the formation of nFGFR1-associated TADs. The makeover of the TAD islands serves to recruit distinct ontogenic programs during the development of the ESC to NCC.

Keywords: FGFR1; HoxA; chromatin structure; embryonic stem cells; neuronal committed cells.

Figure 1

Chromatin interactions occur across the genome in both ESCs and NCCs. (A) Full genome inter- and intra-chromosomal contact map. (B) Chromosome 6 intra-chromosomal contact map. (C) 49.5–54.0 mb chromosome 6 contact map. (D) Genome-wide data overviews of RNA-seq and nFGFR1 binding and interactions in ESCs and NCCs. (E) Data overviews at Chromosome 6 for RNA-seq, nFGFR1 binding, and interactions in ESCs and NCCs. Values are shown for ESCs in blue and NCCs in red.

Figure 2

Interactions correlate with gene coding, regulatory features, and nFGFR1 binding. Interactions, gene coding and regulatory features, and nFGFR1 are predictive of RNA expression. (A,B) PCA on ESC and NCC full genome measurements of correlations between interaction anchor strength, gene coding and regulatory features, and nFGFR1 binding. (C–F) Deep Neural Network machine learning using a two-window approach prediction model of interaction strength, gene coding and regulatory features, and nFGFR1 binding for prediction of RNA expression. (C,D) Two output categories (<1 or >=1 FPKM clusters) or (E,F) 3 categories (<1, 1–30, and >30 FPKM clusters) in ESC and NCC using window-ranges of 1/1–200/240 kb were used for prediction of RNA expression from interaction strength, gene coding and regulatory features, and nFGFR1. A two-window range (ranges in kb) neural network prediction model was applied in which a smaller window as a numerator contains the attributes used for prediction (FGFR1 binding, interaction strength, genome structural features, or combinations of them), and a larger (or equal to) window as a denominator contains region within which the RNA FPKM level is predicted. Analysis was completed using Keras, to build a deep neural network capable of providing high levels of accuracy for FPKM class prediction.

Figure 3

TADs occur in both ESC and NCC. ESC and NCC chr6 50–54 directionality indexes (left) and calculated TADs (right), with an example of directionality index (+) and TAD (*) reorganization during NCC differentiation.

Figure 4

Interaction strength correlates with gene coding, regulatory features, and nFGFR1 in ESC TADs. 200 TADs of 480 kb stacked as rows in descending order by interaction strength within each TAD. Heatmaps show the feature span (gene coding and regulatory features) or strength (interaction strength, RNA FPKM, nFGFR1 binding) throughout each TAD and a mean plot of the heatmap rows in 5 TAD bins. Pearson’s R values were calculated between interaction strength and each of the corresponding attributes in 5 TAD bins. Right panels to each heatmap show average interaction strength sorted feature enrichments across same sized TADs. Z-Score statistics indicate bins which are outside the mean of all bins p < 0.05; p < 0.01; p < 0.005.

Figure 5

Interaction strength correlates with gene coding, regulatory features, and nFGFR1 binding in both ESC and NCC TADs. (A–H) All TADs in ESC (NCC in Figure S6) split in half and aligned by their left and right borders (NCBI directionality). Heatmaps show features enrichments in each TAD. Top and right panels to each heatmap show average location (column means—5 kb bins) and size (row means—5 TAD bins) dependent feature enrichments across all TADs. Z-Score statistics indicate bins which are outside the mean of all bins p < 0.05; p < 0.01; p < 0.005.

Figure 6

Downregulated TADs express stage specific genes and are reorganized during ESC to NCC differentiation. At interacting (q < 0.001) upregulated and downregulated anchor-anchor midpoints: (A) Differential gene expression. (B) TAD overlap. (C) Directionality Index. (D) Gene Ontology Category Enrichment. Two-way ANOVA Tukey Method: ESC versus NCC for each location shown, p < 1⁻²; p < 1⁻⁴; p < 1⁻⁶; p < 1⁻⁸; p < 1⁻¹⁰. At TAD borders aligned left and right based on regulated gene directionality (adjusted p < 0.05): (E) TAD reorganization. (F) Differential gene expression. (G) Differential interaction anchor strength. (H) Differential nFGFR1 binding. (I) Differential gene coding and regulator feature enrichment. Z-Score statistics indicate bins which are outside the mean of all bins p < 0.05; p < 0.01; p < 0.005. (J) Differential chromatin looping. (K) Differential nFGFR1 looping. (L) Differential CTCF looping. Paired T-Test Bonferroni adjusted p < 10⁻⁵; p < 10⁻¹⁵; p < 10⁻²⁵.

Figure 7

Upregulated TADs express stage specific genes and are reorganized during ESC to NCC differentiation. At interacting (q < 0.001) upregulated and downregulated anchor-anchor midpoints: (A) Differential gene expression. (B) TAD overlap. (C) Directionality Index. (D) Gene Ontology Category Enrichment. Two-way ANOVA Tukey Method: ESC versus NCC for each location shown, p < 1⁻²; p < 1⁻⁴; p < 1⁻⁶; p < 1⁻⁸; p < 1⁻¹⁰. At TAD borders aligned left and right based on regulated gene directionality (adjusted p < 0.05): (E) TAD reorganization. (F) Differential gene expression. (G) Differential interaction anchor strength. (H) Differential nFGFR1 binding. (I) Differential gene coding and regulator feature enrichment. Z-Score statistics indicate bins which are outside the mean of all bins p < 0.05; p < 0.01; p < 0.005. (J) Differential chromatin looping. (K) Differential nFGFR1 looping. (L) Differential CTCF looping. Paired T-Test Bonferroni adjusted p < 10⁻⁵; p < 10⁻¹⁵; p < 10⁻²⁵.

Figure 8

HoxA and HoxB clusters connect interchromosomally and reorganize during NCC differentiation. (A,B) Interchromosomal contact map showing locations enriched for interactions between Chr6 and Chr11 (1–10 examples numbered) with the HoxA-B interaction block circled. (C,D) Interchromosomal contact map showing enhanced view of HoxA-B interaction block with (+) indicating the intersection between the HoxA and HoxB cluster midpoints.

Figure 9

HoxA and HoxB cluster intrachromosomal interactions reorganize during ESC to NCC differentiation in parallel with gene expression and nFGFR1 binding changes. Inhibition of nFGFR1 by PD173064 in the HoxA cluster disrupts interactions at early- and mid-cluster anchor points, lowers CTCF binding, and alters gene expression. (A,B) Differential looping. Paired T-Test Bonferroni adjusted p < 10⁻⁵; p < 10⁻¹⁵; p < 10⁻²⁵. (C,D) Differential interaction anchor strength. (E,F) nFGFR1 binding. (G,H) Differential gene expression. (I) Interactions of HoxA1 with downstream HoxA2-A13. (J) nFGFR1 binding. (K) CTCF binding. (L) Gene expression. (M,P) PD fold change of HoxA1:HoxA2-A13 interactions. (N,Q) PD fold change of nFGFR1 binding. (O,R) PD fold change of CTCF binding. (S) PD fold change of gene expression. (T) ESC to NCC fold change with and without PD. Two-way ANOVA Fisher’s LSD Test: p < 0.05 —p < 0.01 =, p < 0.005 ≡, p < 0.001 ≡.

Figure 10

ESC nFGFR1 binding colocalizes with CTCF, NCC nFGFR1 binding targets new locations distinct from CTCF. (A) Immunocytochemistry of ESC and NCC stained with CTCF ab (red) and nFGFR1 ab (green). (B) Enhanced view of single ESC and NCC nuclei. (C) Correlation analysis of CTCF and nFGFR1 colocalization using Fiji ImageJ Coloc 2. (D) Genome Archipelago Model-hubs of transcriptional activity occur throughout the genome as TAD and multi-TAD islands. TADs containing intra-TAD genes looped together form islands alone or with other TADs which are in the vicinity of each other throughout the nucleus. TAD islands integrate multiple co-regulated genes under control of NCC− and NCC+ specific proteins, dependent on the protein binding motifs which are present. The TAD islands and looping components change during ESC to NCC differentiation. Gene co-expression relationships are conserved evolutionarily even as genes move around the genome because maintaining the same protein binding motifs allows the genes to form islands with diverse members of the same co-regulated gene group despite their distinct genomic locations. TAD island formations are controlled by CTCF and pluripotency transcription factors in ESC, and by FGFR1 and neuronal transcription factors in NCC.