High resolution temporal transcriptomics of mouse embryoid body development reveals complex expression dynamics of coding and noncoding loci

Abstract

Cellular responses to stimuli are rapid and continuous and yet the vast majority of investigations of transcriptional responses during developmental transitions typically use long interval time courses; limiting the available interpretive power. Moreover, such experiments typically focus on protein-coding transcripts, ignoring the important impact of long noncoding RNAs. We therefore evaluated coding and noncoding expression dynamics at unprecedented temporal resolution (6-hourly) in differentiating mouse embryonic stem cells and report new insight into molecular processes and genome organization. We present a highly resolved differentiation cascade that exhibits coding and noncoding transcriptional alterations, transcription factor network interactions and alternative splicing events, little of which can be resolved by long-interval developmental time-courses. We describe novel short lived and cycling patterns of gene expression and dissect temporally ordered gene expression changes in response to transcription factors. We elucidate patterns in gene co-expression across the genome, describe asynchronous transcription at bidirectional promoters and functionally annotate known and novel regulatory lncRNAs. These findings highlight the complex and dynamic molecular events underlying mammalian differentiation that can only be observed though a temporally resolved time course.

Figure 1

Global and gene-specific evaluation of augmented temporal resolution in mES differentiation. (A) Schematic of mouse embryonic stem cell (ESC) differentiation into embryoid bodies (EB) over the time course evaluated here. (B) Analysis of the top three principle components (PCs) based on the 2,000 most variable genes from biological duplicate-6 hourly transcriptomes and KEGG pathway enrichment for 500 genes contributing most to each of the top three PCs. (C) Expression profiles of genes associated with pluripotency, primitive streak formation and cell specialization.

Figure 2

Insights into regulatory and gene expression kinetics. (A and B) Observable regulatory network dynamics at 24- and 6-hourly measures with Otx2 and Pou5f1 target containing profiles annotated and in bold, See Supplementary Fig. S2 for full figure. Transcriptomes at 24- (top) and 6-hourly (bottom) were subjected to DREM analysis of mouse TF/target gene interactions. Lines represent the median fold change (relative to time 0) of grouped TF target genes- representing activity of the TF itself, line colors are assigned by branch and are not comparable between panels. A p-value cutoff of 0.001 was applied to calculating divergent TF activity (splits). (C) Expression of the key transcription factors Pou5F1 (Oct4) and Otx2. Red and blue boxes correspond to the time points highlighted in part A. (D) Distribution of the number of genes and the time delay required to meet a maximum correlation (>0.8) between gene targets of Pou5f1 and Pou5f1 itself compared to 95% quantiles of 500 random gene selections. (E) Two k-means clusters of short-lived RNA (slRNA) genes displaying differential expression without changes at 24-hourly time points (adj. p < 0.0001).

Figure 3

Analysis of gene co-expression patterns using augmented temporal resolution. (A) Smoothed scatter plot showing the correlation coefficient across the time course vs. distance between transcriptional start sites (TSS) of bidirectional gene pairs. Blue indicates no gene pairs; yellow and red indicate increasing numbers of pairs sharing similar properties. (B) Expression patterns of example bidirectional genes of the same or different gene biotype. Spearman’s correlation coefficient is reported for each pair. (C) Genomic location (circos) and expression pattern (line plot) of two independent co-expressed groups of 5 or more contiguous genes sharing correlated expression (r > 0.5).

Figure 4

Augmented temporal resolution of ncRNA expression in cellular differentiation. (A) Hierarchical clustering of lncRNAs (dark blue) with time-dependent protein coding genes (light blue) by their expression patterns over time. Dendrogram was manually colored to reflect gene expression levels of the top-level clusters. (B) K-means clustered expression profiles of protein coding genes compared to the same number of lncRNA gene expression clusters. Common profiles are marked with arrows. (C) Expression profiles of four lncRNAs predicted to have regulatory roles in ES development as well as the genome location & pathways enriched in their gene targets.  Malat1 and IRX3os display a positive association with their targets, whereas 1700057H21Rik and 1700042O10Rik have a putative repressive impact.