Analysis of dynamic changes in retinoid-induced transcription and epigenetic profiles of murine Hox clusters in ES cells

Abstract

The clustered Hox genes, which are highly conserved across metazoans, encode homeodomain-containing transcription factors that provide a blueprint for segmental identity along the body axis. Recent studies have underscored that in addition to encoding Hox genes, the homeotic clusters contain key noncoding RNA genes that play a central role in development. In this study, we have taken advantage of genome-wide approaches to provide a detailed analysis of retinoic acid (RA)-induced transcriptional and epigenetic changes within the homeotic clusters of mouse embryonic stem cells. Although there is a general colinear response, our analyses suggest a lack of strict colinearity for several genes in the HoxA and HoxB clusters. We have identified transcribed novel noncoding RNAs (ncRNAs) and their cis-regulatory elements that function in response to RA and demonstrated that the expression of these ncRNAs from both strands represent some of the most rapidly induced transcripts in ES cells. Finally, we have provided dynamic analyses of chromatin modifications for the coding and noncoding genes expressed upon activation and suggest that active transcription can occur in the presence of chromatin modifications and machineries associated with repressed transcription state over the clusters. Overall, our data provide a resource for a better understanding of the dynamic nature of the coding and noncoding transcripts and their associated chromatin marks in the regulation of homeotic gene transcription during development.

Figure 1.

Analysis of changes in Hox and global gene expression during RA-induced differentiation of mouse ES cells. (A) Flow chart outlining overall experimental strategy. (B) Heatmap showing global changes in gene expression upon RA induction compared with uninduced ES cells as analyzed on Affymetrix Mouse Genome 430 2.0 arrays. The lower panel shows the distinct induction profiles of six clusters, identified by k-means clustering. Expression values are the average value from three independent biological replicates. The middle panel displays a heatmap of global changes in gene expression upon RA-induced differentiation; k-means clusters are indicated on the right. Only clusters with an absolute value of cluster mean >0.5 are shown. The upper panel shows changes in the expression profile for 15 of the most rapidly induced genes in cluster 1, which includes Hox genes and genes for the cofactors (Meis). (C) Heatmap of relative changes in Hox gene expression upon RA induction compared with uninduced ES cells as analyzed on Affymetrix Mouse Genome 430 2.0 arrays. Genes in the HoxA and HoxB clusters display a more rapid and robust induction than those of the HoxC and HoxD clusters. Several genes have multiple probes, and results are shown for each. (D) Temporal changes in Hox gene expression induced by RA quantitated by TLDA qPCR microfluidics cards. All data points are the average of three biological and two technical replicates. The y-axis in all clusters is shown on the same scale and illustrates relative levels of induction between Hox clusters. (E) RNA-seq analysis of Hox gene expression in RA treated ES cells compared to uninduced cells. Fold changes are shown as a heatmap.

Figure 2.

Single-molecule fluorescence in situ hybridization (FISH) analysis of gene expression in ES cells. (A,B) Quantification of single-color FISH data using either the Stellaris multiprobe approach (A) or the hybridization chain reaction (HCR) method (B) to measure expression of Cyp26a1. Each point represents total integrated intensity in arbitrary units per cell analyzed for uninduced colonies or colonies induced with RA for 4 or 24 h. The y-axis is on a log₁₀ scale. A 90% threshold of the uninduced distribution is used to estimate levels above which cells are “on” for expression of a gene. (C) Dual-color Stellaris FISH analysis using Cyp26a1-q570 and Hoxa1-q670 probes. (D) Dual-color Stellaris FISH analysis using Cyp26a1-q670 and Hoxb1-q570 probes. In C and D, total integrated intensity in arbitrary units is determined for each color in each cell. Lines on the 2D plot are generated corresponding to the 90% level of the uninduced cells for that color. The white box represents cells that are scored as off for both genes, the red box represents cells on for both genes, and the green and blues boxes represent cells that are positive for a single probe. For A–D, all experiments used a minimum of 169 total cells from at least nine different colonies for analysis. The fraction of positive cells for each probe and time point are indicated at the bottom of the respective panels.

Figure 3.

Global transcriptional activity in and around the four Hox clusters analyzed by Agilent tiling arrays. Heatmap showing global changes in gene expression over the time course of RA induction compared with uninduced ES cells as analyzed on custom Agilent 2x105K Hox tiling arrays. mRNA is labeled with Cy3. Heatmaps are generated on the Integrative Genomics Viewer (IGV) 1.5, where intensity of heat represents relative expression level. Upper and lower panels in each cluster represent sense and antisense strands, respectively. The x-axis denotes the chromosomal position and location of respective Hox genes and ncRNAs. The y-axis denotes the length of RA treatment in h.

Figure 4.

Characterization of Heater (HoxA EArly Transcribed REgion) transcripts. (A) Transcription from the Heater region analyzed by tilling arrays, RNA-seq, and ENCODE. The top of the panel shows rapid induction of Heater transcripts between 2 and 24 h of RA treatment detected by tiling array profiles. W and C represent Watson and Crick strands, respectively. In the middle of the panel is evidence of Heater expression in adult kidney and cultured cells based on ENCODE data. + and − indicate opposing strands. RNA-seq data in uninduced and 24-h RA-treated ES cells validate expression profiles observed in tiling arrays. At the bottom, multiple Heater transcripts, from Halr1 and Halr1os1, are shown schematically with the respective transcription start sites (TSS) and direction of transcription as indicated by arrows. The pink box denotes a region with a large number of unspliced transcripts generated from TSS4. (B) Heatmap of qPCR quantitation of Halr1, Halr1os1, Hotairm1, Hotairm2, and Hobbit1 transcripts in RA-induced ES cells. Levels of transcripts are compared against respective transcript levels in uninduced ES cells. Heater, Hotairm1, and Hotairm2 are induced at comparable levels, whereas Hobbit1 shows a higher level of induction in RA-induced differentiation. (C) Heatmap of qPCR quantitation of Halr1, Halr1os1, Hotairm1, Hotairm2, and Hobbit1 transcripts in developing mouse embryos. Levels of transcripts are compared against respective transcript levels in 10 dpc embryo. Heater and Hobbit1 are induced at comparable levels, whereas Hotairm1 and Hotairm2 show higher levels of induction during mouse embryonic development. Green represents down-regulation; red shows up-regulation. (D) Response of Halr1, Halr1os1, Hotairm1, Hotairm2, and Hobbit1 to RA in mouse embryos. The relative response to RA is calculated as fold change in transcript level after treating 10.0 dpc mouse embryos with RA compared to untreated mouse embryos at the same stage.

Figure 5.

ChIP-on-chip analysis of changes in the epigenetic state and retinoid receptor occupancy of HoxA and HoxB clusters during RA-induced differentiation. H3K27me3 is used as a repressive mark; H3K4me3 is used as a mark for an active chromatin state, and Pol II is used as a mark of active transcription. Tracks are configured by using windowing function as mean and smoothing windows as 0 pixels in the UCSC Genome Browser. All time points for a given antibody are normalized with the same y-axis, and the specific range of each y-axis is shown for respective uninduced samples. The schematic at the top indicates the relative positions of Hox genes, microRNAs, and noncoding transcripts. Within 4 h of RA treatment, there is a rapid gain of H3K4me3 and Pol II, whereas H3K27me3 is gradually lost over 36 h. Many major changes in Pol II occupancy and gain of H3K4me3 are related to Heater, Hotairm1, and Hobbit1 noncoding transcripts. Distinct dynamic changes in retinoic acid receptors and the NCOR corepressor is observed over HoxA and HoxB cluster.

Figure 6.

Analysis of changes in the epigenetic state and retinoid receptor occupancy in and around the Heater region during RA-induced differentiation. (A) Occupancy of retinoic acid receptors (RXRA, RARA, RARB, and RARG) and the NCOR corepressor in and around the Heater region. ChIP-on-chip analysis shows a large 2.5-kb region (H-AR1) bound by RAR/RXRs and NCOR upstream of Heater and a smaller domain (H-AR2) region downstream from uninduced and RA-treated ES cells. The schematic at the top shows the relative positions of these regions to Heater and the relative positions of predicted consensus Direct Repeat motifs recognized by retinoid receptors (see Supplemental Fig. S7). (B) Dynamic occupancy of Pol II and components of an elongation complex in the Heater region. H3K27me3 is used as a repressive mark, H3K4me3 is used as a mark for promoters and active chromatin state, and Pol II is used as a mark of active transcription. Along with Pol II, a bivalent mark formed by H3K4me3 and H3K27me3 is noticeable over TSS1-TSS4. There is a rapid recruitment of the transcription elongation factors (ELL2, AFF4, and CDK9) upon RA induction. Tracks were configured by using a windowing function as mean and smoothing windows as 10 pixels. At the bottom, multiple Heater transcripts are shown schematically with the respective transcription start sites (TSS) along with RNA-seq at 24 h. The direction of transcription is indicated by arrows. The pink box denotes a region with a large number of unspliced transcripts generated from TSS4. A putative TSS is also apparent in the H-AR2 region.

Figure 7.

Gradual and progressive loss of repressive marks from the four mouse Hox clusters during RA-induced differentiation. (A) ChIP-on-chip analysis shows a gradual loss of the H3K27me3 repressive mark and SUZ12 occupancy over each Hox cluster upon RA treatment. Although many of the Hox genes are expressed early in the time course, over a whole cluster, H3K27me3 is greatly reduced by 36 h of RA treatment. The gradual loss of repressive marks is observed from anterior to posterior genes in a Hox cluster over a time course correlating with colinearity. (B) Kinetics of reduction of SUZ12 occupancy over TSS of HoxA and HoxB genes in paralogy groups 1 and 9 during RA-induced differentiation of ES cells. Anterior Hox genes rapidly lose SUZ12 over their TSS, as illustrated by changes for Hoxa1 and Hoxb1 at 2 h of RA treatment; whereas posterior genes show little change over their TSS in this time frame. The differences in the kinetics of loss of SUZ12 between genes correlates with their respective time of activation. A 500-bp region around the TSS is shown, and a 50-bp region around TSS is marked by a light blue band. The y-axis shows relative levels of occupancy of SUZ12.