|Home | About | Journals | Submit | Contact Us | Français|
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.
Over the past decade, transcriptomic investigations into the of nature embryonic stem cell (ESC) differentiation have elucidated key biochemical features of stemness and differentiation. Increasingly, it has become apparent that understanding the dynamics and coordination of gene expression signatures over time during the key phases of differentiation is critical to adequate characterization of fundamental biological processes.
Mouse ESC differentiation is a highly complex cascade of gene expression changes that allow single pluripotent cells in culture to progress to an organoid composed of cells reflecting three germ lineages within only five days. The spontaneous differentiation of these cells in culture has provided key insights into the developmental processes underlying the generation of the primary germ cell layers1. Microarray and RNA sequencing have provided a means to characterize the molecular transitions in gene expression underlying ESC biology and more recently single cell transcriptomic studies have provided the first glimpses into the molecular history of these cells2. However, it is clear that much more of the transcriptional landscape of ESC remains to be elucidated3.
Access to new technologies, such as massively parallel sequencing (MPS), has led to a dramatic increase in our knowledge of the mammalian transcriptome. Early genomic tiling array analysis indicated that most of the genome was transcribed into RNA4. MPS of the transcriptome validated this observation and revealed that the majority of the mammalian genome is pervasively transcribed as interlaced and overlapping RNAs5, many of which lack protein-coding potential6. The large number of long-noncoding transcripts (lncRNA) has become the focus of significant interest due to their exquisite cell type specific expression7, potent biological function8, 9, and rapid transactivation of cellular processes10. However, in general, lncRNAs are lowly expressed and short lived11, possibly because, unlike mRNAs that require translation, are able to exert their function directly. These qualities obfuscate their identification and characterization with traditional approaches that are tuned to the properties of mRNAs12. Owing to the relative infancy of the field, the vast majority of noncoding transcripts are of unknown function13. Additionally, the expression patterns of these genes imply that their function is dependent on cellular context and likely regulatory8, thus the identification of these molecules and the context in which they act remains a research priority14.
Various expression profiling studies, using both microarrays and RNA-seq15–18, have been used to explore the molecular changes occurring during ES cell development, typically at 24-hourly or more. This potentially has lead to incomplete gene expression relationships through the phenomenon of temporal aggregation bias whereby each time point is assumed to represent all the signaling changes occurring in that time window19. In contrast to single cell based approaches- which provide insight into the state of individual cells - examinations of whole cell populations provides system-wide behavior and a practical means to explore gene expression dynamics across time. The combination of these techniques has recently shed light on the molecular framework of cellular differentiation20. Higher temporal resolution has also shown rapid induction (within two hours of retinoic acid stimulation) of lncRNAs associated with the HOX locus21. Furthermore, high temporal resolution has provided valuable insights into transcriptional annotation and regulation in drosophila22, 23, Xenopus24 and C.elegans25.
Here we show that additional temporal resolution of the global transcriptome in spontaneously differentiating mESC cells following LIF withdrawal enables the capture of the rapid and complex dynamic regulatory and noncoding changes occurring during ES development. We analyzed the transcriptome of differentiating mouse ESCs at six-hourly intervals over a five-day period, over which time the three primordial germ layers are specified. Using this fine-resolution temporal sampling approach, we identify significant transitions in the transcriptome and large-scale shifts in observable transcription factor activities that could not be observed at 24 hourly sampling periods. Moreover, we identify entirely novel coding and noncoding transcripts that are expressed only within specific sub-24-hour window. By leveraging the high sampling frequency of the data, we are able to both accurately recapitulate known regulatory cascades in ES development and predict and refine others. Finally, using correlative approaches, we can infer functions for uncharacterized lncRNAs and predict the regulatory centers across the genome that coordinate early development.
A median 42-million, paired-end 100-bp reads (Supplementary Fig. S1A) were mapped from stranded, poly-A derived cDNA libraries derived from biological duplicate, six-hourly time courses of mESC differentiation over five days where key differentiation programs occur (0–120hours, Fig. 1A). Transcript-level expression data was generated as previously described26, then normalized for library size and transformed for data visualization and differential gene expression analysis. Evaluation of 24 hourly time points indicated that our data was comparable to previously published data in a similar model27 (Supplementary Fig. S1B). An interactive gene expression portal was created to visualise this data (https://betsig.shinyapps.io/paper_plots).
To assess the reproducibility and provide confidence in the biological validity of the global transcriptome trends, a principle components analysis (PCA) was performed on the 2,000 most variable genes (Fig. 1B). This analysis indicated that biological replicates clustered closely, indicating that synchrony was retained, and that the major contributor to the determination of variance was explained by time. Deconvolution of the dimensions yielded time-dependent expression (in the first dimension) of genes enriched in focal adhesion/ECM interactions KEGG pathways. Interestingly, the second dimension deconvolution (PC2), in which undifferentiated ESCs resemble the more differentiated embryoblast, yielded genes enriched in MAPK-signaling and cancer pathways, implying that the process of differentiation involves a partial reacquisition of a mitogenic signaling. In the third component (PC3), in which the undifferentiated ES cell is separate, the axon-guidance pathway was enriched. We then evaluated expression patterns of genes associated with pluripotency, primitive streak formation and cell specialization (Fig. 1C). We observed that, although the gene expression patterns were broadly consistent with published studies (Supplementary Fig. S1B), there were changes in expression on less than 24 hourly timeframes that could not be attributed to stochastic expression changes (within the top 5% of deviation of all genes from expression values loess-smoothed over 24hours). Similarly, we observed that no single 24-hourly measure was representative of the average expression over that day (Mann-Whitney U p. adj. <<<0.0001) and that more than 1,000 genes displayed a more than 2-fold difference mostly in the first 24hours of differentiation (Supplementary Fig. 1C,D). Therefore we conclude that 24-houly expression profiles are unable to capture the intervening expression changes and that 6-hourly measurements reduce the phenomenon of temporal aggregation bias19, providing enhanced resolution of transcriptional changes in this system.
To evaluate characteristics of sub-24 hour gene expression patterns in in the transcriptome of developing ESCs, we observed that, compared to 24hour time points, 417 more genes had counts data considered sufficient for differential gene expression analysis (>1 CPM in at least two samples); this was associated with a relative increase in detected noncoding genes (13% (588 vs. 520); defined as antisense, lincRNA and processed transcript biotypes) over protein coding, (2% (13336 vs. 13036) despite being underrepresented in the total pool (chi-squared p value<0.001, Supplementary Table S1 and Supplementary Fig. S1E)). The additional time points allowed the assembly of 58% more novel multiexonic intergenic, antisense and intronic noncoding RNAs from the data - indicating that a substantial proportion of noncoding transcripts are present on timescales much shorter than 24hours. These results indicate that enhanced temporal resolution reduces the phenomenon of temporal aggregation bias and allows the observation of more distinct cell expression states than typical time-courses.
Increased sampling frequency can provide a powerful insight into understanding of the contribution of gene regulatory networks to cellular differentiation22. We utilized the DREM v2 analysis tool28 to evaluate transcription factor (TF) target gene expression patterns. Divergence of gene targets responsive to groups of TF at each time point, either 24-hourly or 6-hourly (Fig. 2A,B) was shown if the overall difference was significant at p<0.001. Compared to 24-hourly, the observed complexity was significantly higher, especially in the first 48hours. We observed that significant changes in gene regulation occurred continuously within the 24-hour windows. Most notably, first 24hours following depart from pluripotency resembles an ordered cascade of TF activity (Figs 2A and S2A) with large-scale changes in TF activity at 12, 18 and 24hours; of which little can be deduced measuring at just 24hours (Figs 2B and S2B). Focusing on the interplay between two key transcription factors (Otx2 and Pou5f1/Oct429, Fig. 2A), we observed a rapid rise in Otx2 activity in the first six hours and stable Pou5f1activity for the first 24hours (Red Box). Otx2 activity did not coincide with mRNA expression of the factor itself (Fig. 2A vs. vs.C),C), although previous studies have observed increased in Otx2 protein expression within 3-hours of differentiation29, however periodic drops in Pou5f1 mRNA expression appeared to coincide with decreases in POU5F1 target genes, we calculated the time taken for Pou5f1 expression to result in changes in highly positively correlated (r>0.8) target genes using a cross-correlation approach similar to ref. 30. We then evaluated how these “delays” enriched for certain Reactome pathways (Fig. 2D). We found rapid effects for targets enriched for “gene expression”- and a delayed effect on “cell cycle” pathways compared to a null distribution produced by 500 random “target” selections (grey). These were similarly observed in the DREM GO-term enrichment tool for Pou5F1 targets decreasing in expression at 42 (early- Transcription Factor Activity) and 54hours (late- Epithelial Proliferation; Fig. 2A, Blue Box & Supplementary Fig. S2C) and associated with the decrease in Pou5F1 expression (Fig. 2C, Blue Box). Importantly, Pou5F1 mRNA and protein expression are temporally correlated29. This result implies that TF-target genes may be activated in an ordered- time dependent fashion. To explore this more broadly, we evaluated other TF-target gene temporal dynamics for other TFs that exhibited strong positive or negative correlations between the TF and their target genes. We found evidence of highly structured TF-target expression patterns in time for negatively correlated Pou5f1 and Suz12 targets, as well as positively correlated Nanog, Myc, Sox2 and Suz12 targets (Supplementary Fig. S3).
These observations of precise temporal ordering of transcriptional events emphasize the importance of factoring time delays into understanding gene regulatory networks31. This highlights the capacity of increased temporal resolution to directly identify regulator-target gene interactions instead of relying on inference; which is common in cross-correlation approaches.
Having established that the increased temporal resolution markedly improves the molecular framework for evaluating the contribution of gene expression to ES differentiation, we next sought to identify gene expression signatures previously unable to be resolved using lower temporal resolution. For each 24-hour period, we identified genes that were differentially expressed between 0 and 6, 12 and 18hours but not between any 24-hourly measures (Supplementary Fig. S2D). We identified 1,135 genes with significant changes in gene expression that were unchanged between any 24-hourly comparison (adjusted p<0.0001). Of these, 354 were differentially expressed for more than half of the corresponding 24-hour window, mostly in the first and last 24-hour periods. These genes were described as short-lived RNAs (slRNAs). slRNA expression patterns over the first 24hours of differentiation were found to be positively correlated with the same time window of retinoic acid directed differentiation21 (Supplementary Fig. S2E) implying that these genes may form part of the early response to differentiation signals. K-means clustering and KEGG pathway analysis of the expression profiles of these genes (Fig. 2E) revealed enrichment in genes associated with the spliceosome (adjusted p<0.05) dramatically decreasing in expression over the first 24hours before returning slowly to baseline. To examine whether this impacted gene-splicing patterns, we employed a differential exon (DEX) analysis between consecutive six-hourly time points and counted the number of genes displaying DEX usage (adjusted p value<0.01 Fig. 2E). Consistent with previous studies, the alternate splicing was most highly associated with cell differentiation32 (Fig. 2E). Increased temporal resolution has elucidated that these changes happen very rapidly (majority of changes in the first six hours), and that slRNAs may be involved in suppressing the alternate splicing of genes and limiting transcriptional plasticity.
Some slRNAs appeared to have periodic expression profiles. We thus sought to uncover periodic expression patterns genome-wide, by applying a fast-Fourier transformation to our data (see Methods). Periodogram analysis was utilized to ascertain the dominant cycling period for each gene. We found 137 genes, which we termed cycling RNAs (cycRNAs), sharing the same dominant cycling period of less than 36hours in both biological replicate experiments (Supplementary Table S2). Supporting the efficacy of the approach, we found Clock, which encodes a key regulator of circadian rhythm in mammals, to have a period of 24.2hours. We identified 20 genes that displayed characteristics of both slRNAs and cycRNAs (Supplementary Fig. S2F), including Ewsr1 and Clk1, involved in gene splicing33, 34 as well as five uncharacterized lncRNAs. Given the highly specific expression patterns in this context, we propose these genes may similarly have roles in maintaining or establishing biological rhythms. Together these investigations show that the augmented temporal resolution approach provides access to gain insights from regulatory pathways by identifying transitions in expression that would otherwise have remained hidden.
Evaluating gene transcription at high temporal resolution in a highly dynamic process such as ES development, we anticipated that it might be feasible to dissect structural gene regulation within a given locus. To explore this possibility, we examined expression arising from transcripts that are oriented head-to-head as so-called bidirectional pairs35, 36. Interestingly, we observed that the antisense transcript for Evx1 (Fig. 1C) displayed a previously unobserved15 increase in expression in the first 24hours after departure from pluripotency that was reflected in its paired protein coding gene Evx1 (Supplementary Fig. S4A), highlighting the increased power of frequent sampling over time. In total, we identified 1,251 gene pairs with bidirectional transcriptional start sites (TSS) within 2,000bp and evaluated correlation coefficients across the time course, distance between TSS and median expression values. Consistent with other studies, we found expression correlation more positive for bidirectional gene pairs than random transcript pairs35 (Supplementary Fig. S4B). We were also able to show that the distance between TSS of highly correlated bidirectional gene promoters is typically less than 500bp (Fig. 3A), consistent with a common regulatory domain. Highly correlated or anti-correlated genes pairs displayed differences in total gene expression, particularly with discordant gene biotypes (Fig. 3B). We found that protein coding gene pairs were more likely to be of similar expression levels and positively correlated (Mann-Whitney p<0.05) than protein coding/noncoding pairs (Supplementary Fig. S4C). Applying a variant of the temporal offset analysis used to measure TF- gene target delays, we calculated the time taken and defined the apparent driver gene type for peak correlation in coding/noncoding bidirectional pairs (Supplementary Fig. S4D,E). This did not reveal a generalized bias in either time taken or particular “driving” gene type. However, this approach shows that the lncRNA Hotairm1, required for activation of Hoxa137, appears to have a six-hour delay between its expression changes and HoxA1. We present evidence of other examples of lncRNA-led expression of protein coding genes in small numbers of bidirectional pairs (Supplementary Fig. S5).
To investigate whether the strong correlative potential between gene pairs could facilitate the identification of regions of the genome that are coordinately regulated38, we scanned across the genome for regions containing five or more contiguous genes that were co-expressed (r>0.5). This revealed 59 regions with a mean size of 821kb -each containing 5–14 genes (mean of 6) genes. To examine the higher order chromatin architecture of these regions, we compared these regions to published data on topological associated domains (TADs) for mouse ESCs39. We found that the majority of the regions were each contained within a single TAD (Supplementary Fig. S4F), increasing the likelihood for a common regulatory architecture. Evaluation of gene-expression patterns across these regions revealed evidence of high co-expression at both the inter- and intra-chromosomal levels (Supplementary Fig. S4G). We assembled a map of regions of the mouse genome displaying high levels of clustered co-expression (Fig. 3C) by comparing the expression profiles of the regions. Two independent modules were identified with distinct decreasing (green)- and increasing (blue) expression patterns with differentiation. Given the independent location and expression patterns of these clusters, we suggest they may form core expression-factories of cellular differentiation. In support of this notion, this analysis identified the a region -associated with the “increasing module”- containing the imprinting locus of H19, Igf2, Tnn3 and Mrpl23 40 (Supplementary Fig. S4H); previously shown to be activated in concert during early stem cell differentiation41.
These investigations illustrate how analysis of high-resolution temporal transcriptomic data provides an independent and convenient approach (relying only RNA-Seq) to guide the partitioning of the genome into regulatory domains.
Having shown that rapid changes in lncRNAs are a key feature of ES differentiation, and that co-expression analysis is a powerful tool for understanding gene regulation with augmented temporal resolution, we sought to unravel the roles that lncRNAs might play in ESC differentiation.
Analysis of gene annotations yielded confident expression data for 588 lncRNA genes at six-hourly resolution (520 for 24-hourly, Supplementary Table S1). Indeed, added temporal resolution increased information of all noncoding transcript biotypes indicating that a proportion of these genes were only present for a short duration in this system. It is important to note that these do not represent the entirety of lncRNA expression in this process since the depth of sequencing was not standard for lncRNA coverage42 and the poly-A selection for cDNA would have missed enhancer RNAs, miRNA precursors and Nuclease P processed lncRNAs43. Clustering the observed lncRNA expression patterns with time-dependent protein coding gene expression showed that lncRNAs were enriched at lower expression levels and shared related expression profiles to protein coding genes (Fig. 4A). This relationship was further examined whereby K-means clustering of these expression profiles compared to clustering of a similar number of time-dependent protein coding genes (Figs 4B and S6A) revealed clusters of lncRNA genes resembling gene expression patterns associated with stemness (cluster a) primitive streak formation (cluster b) and WNT signaling (cluster c)15. As co-expression has been illustrated to provide valuable insight into lncRNA function18, the additional correlative strength afforded by this study is anticipated to more reliably guide the functional association of these lncRNAs with these processes. The dynamics of lncRNA expression observed here indicate that future studies using RNA capture sequencing or higher sequencing depth of Ribosome depleted RNA will provide more comprehensive insights into the role of lncRNA in the molecular events underlying cell differentiation.
As lncRNAs often exert their function through guiding or assembling transcriptional machinery, we sought to identify potential regulatory lncRNAs in this system. We selected 50 highly or variably expressed lncRNAs (Fig. 4A) and tested for evidence of gene regulatory behavior across the transcriptome. The temporal resolution allowed the use of time to resolve precedence, thus adding weight to a potential causal relationship. Since lncRNAs typically exert their function as a transcript, we set a maximum time offset of 18hours to avoid secondary (altered protein level) effects and examined patterns in the predicted gene targets of lncRNAs (r>0.8, divided by positive or negative associations). Reactome pathway analysis revealed that 11 of these lncRNAs (including well characterized lncRNAs, Supplementary Fig. S6B and C) were potentially involved in regulating networks of genes associated with key developmental processes (p.adj <0.05, Supplementary Fig. S6C). These analyses assigned target gene networks consistent with characterized lncRNA biological functions for Malat1 (oncogenic)44, Neat1 & Rian (association with gene repression)45 and Meg3 (tumour suppressor)46. Interestingly, these data suggest that the pro-tumorigenic function of Malat1 may be mediated through facilitating the increase of MAPK signaling molecules. Importantly, these data also provide testable evidence for seven previously uncharacterized lncRNAs role in ES development and describes a map of regulatory interactions potentially driven by lncRNAs (Fig. 4C) whereby lncRNAs expression may impact coding gene expression across the genome. The identification of lncRNAs with a predicted biological role is important for unraveling lncRNA function, providing candidate functional lncRNAs and providing a level of molecular detail that is currently lacking in many lncRNA studies.
Transcriptional regulation of key biological events is a key feature in understanding the complexity of cellular processes. Here we describe a detailed transcriptomic resource for research in cellular development, a framework for unraveling this detail and identifying new targets for analysis. We also present a comprehensively detailed survey of noncoding transcripts throughout early stem cell development. We have identified many previously uncharacterized noncoding RNAs with potentially pivotal roles in cellular differentiation. This will provide a valuable tool for researchers unraveling the transcriptional complexity of cellular differentiation.
The understanding of molecular events underlying the departure from pluripotency has been determined by the extant knowledge of how biological functions are exerted – often measured at 24 hourly or greater intervals. We hypothesized that interpretations of this model were missing detail in light of evidence indicating the unforeseen dynamics in RNA biology and regulation. By probing this detail with finer time distinctions, we show that gene expression profiles of well-characterized genes display significant variation of expression levels and that more detail can still be gleaned with increased sampling frequency (Supplementary Fig. S7A and B). Importantly, these variations are manifest in a significantly more complex gene regulatory framework. This is consistent with a reduction in temporal aggregation bias19 and highlights early array-based investigations in yeast demonstrating the importance of sufficient temporal resolution in understanding gene expression patterns47. As such, much detail is likely missing from other systems that involve a change in phenotype or cellular behavior. With large-scale transcriptomic analyses becoming increasingly accessible, it is opportune to revisit other well-studied transitions with the view of improving understanding and applicability of their results rather than relying on presuppositions about gene expression patterns48.
We have shown the benefit of frequent sampling over time in observing the transcription of genes that are observable only within sub-24 hourly windows. This approach highlights the importance of taking into account the presence of short-lived transcripts and shows that cells express more of the transcriptome in a time-dependent fashion. To this end, we have identified rapid changing and periodically expressed genes, which we term short-lived (slRNA) and cycling (cycRNA), that were unobservable outside this framework. That many slRNAs exhibited changes in expression over the first 24hours of differentiation is consistent with rapid initial cellular response to stimuli21, 49. Indeed, it is likely that significant gene expression changes- especially noncoding- occur on timeframes shorter than those presented that may not be amenable to optimal time point prediction strategies48. By probing deeper into time-dependent gene transcription-possibly by interpolating available datasets-47 it will be possible to uncover further complexity underlying cellular plasticity and gene regulation. These observations reinforce the concept that adequate temporal resolution is vital for describing biological transitions- for example in dissecting primary from follow on effects in gene knockdown studies – and that end-point analysis likely does not reflect the complex biology of phenotype changes.
Similarly, by using time to separate the order of gene transcription, we have been able to predict local gene regulation across the genome. We have been able to observe concerted gene expression (in trans) of hundreds of genes separated by large genome differences (in cis). Typical studies of this nature involve correlative analysis requiring large samples sizes and resources50. We have instead leveraged the time axis to achieve these as well as discriminate driver from passenger molecular events. This has allowed the estimation of the time delay for changes in expression of regulatory molecules to manifest in changes in their target gene transcription and we have been able to unravel a potentially complex network of gene profiles responding to lncRNA transcription. Finally, we have been able to use an integrated biological system to draw strong associations in trans relationships with bidirectional promoters. Typically these associations are observed by using thousands of gene expression profiles, yet here we have been able to do so with only 42 transcriptomes (duplicate time courses of 21 points each).
Biological duplicate, low passage number (P18) W9.5 ESCs were cultured and differentiated as described previously15, 51. Cultures were harvested every six hours from the induction of differentiation to 120hours post differentiation induction. Total RNA from cultures was purified using Trizol (Life Technologies) and DNase treatment was performed by RQ1 DNase (Promega) according to the manufacturer’s instructions. RNA integrity was measured on a Bioanalyzer RNA Nano chip (Agilent). RNA-Seq library preparation and sequencing of Poly-A-NGS libraries generated from 500ng total RNA using SureSelect Strand Specific RNA Library Preparation Kit (Agilent) were performed according to the manufacturer’s instructions at the same time to minimize batch effect. Paired-end libraries were sequenced to the first 100bp on a HiSeq 2500 (Illumina) on High Output Mode.
Library sequencing quality was determined using FastQC (Babraham Bioinformatics) and FastQ Screen (Babraham Bioinformatics). Illumina adaptor sequence and low quality read trimming (read pair removed if <20 base pairs) was performed using Trim Galore! (Babraham Bioinformatics: www.bioinformatics.babraham.ac.uk/). Tophat252 was used to align reads to the December 2011 release of the mouse reference genome (mm10) as outlined by Anders et al.26. Read counts data corresponding to GENCODE vM2 transcript annotations were generated using HTSeq53. de novo transcript assembly was performed on each merged BAM file using Cufflinks’ reference annotation based transcript (RABT) assembly54, using the Gencode vM2 transcriptome55 as a guide (options: -u -I 500000 -j 1.0 -F 0.005-trim-3-dropoff-frac 0.05 -g gencode.vM2.annotation.gtf–library-type fr-firststrand). Transcript assemblies were then merged using Cuffmerge56 using default parameters, and compared to the Gencode vM2 reference transcriptome using Cuffcompare56. Novel transcripts with a Cuffcompare class code of j, i, o, u or x were filtered using three steps to find novel lncRNAs. First, a Browser Extensible Data (BED) format file was generated using a python script (https://gist.github.com/davidliwei/1155568) and any single exon transcripts were removed. Second, the FASTA-formatted sequence for each transcript was obtained using BEDTools57, the nucleotide (nt) length and open reading frame (ORF) size found using Perl scripts, and those with a length less than 200 nt or a ORF size greater than 300 nt were removed. Lastly, transcript sequences were submitted to Coding Potential Calculator (CPC)58, and those with a coding potential of >0 were removed.
All analyses were performed in the R Statistical Environment59. Briefly, counts data were background corrected and normalized for library size using edgeR60, then transformed using voom61 for differential expression analysis using LIMMA62. Transcription Factor (TF) activity was inferred from gene expression data using DREM28 with a branching P-value of 0.001 based on curated TF-target gene lists associated with mouse ESC differentiation from ChEA63. TF-target gene was calculated by maximal Pearson’s correlation coefficient of >0.8 using a custom autocorrelation analysis and verified with the “ccf” function in R. Gene differential exon (DEX) usage was analyzed by DEXSeq64 on vM2 gene annotations using default settings and an adjusted p value cutoff of 0.001 for DEX between biological duplicates at each consecutive time-point. Genome position analyses were performed using genomic ranges65 based on vM2 annotations imported with ‘rtracklayer’66 and Pearson’s correlation coefficient of gene expression Bidirectional genes were defined as two genes with expression data on opposing strands with <2000 bp between the transcriptional start sites (TSS). Co-expressed gene clusters were defined as >5 contiguous genes with expression data displaying a Pearson’s Correlation Coefficient of >0.5 with neighbouring genes. Cluster co-expression data was visualized with corrplot67 and Cytoscape (v3.1.0)68, location of related clusters was visualized by Circos69. Gene expression periodicity was measured on 120 interpolated expression values70 for each replicate time series using GeneCycle71, candidate periodically expressed genes were identified as having the same calculated dominant cycling frequency between biological replicates. Time-dependent expression signatures were established using maSigPro72 with a replicate correlation coefficient cutoff of 0.8. Target genes of potential regulatory (top 50 most highly and/or variably expressed) lncRNAs were identified using the GeneReg package73 on 100 point-interpolated expression data based on fitted expression values between duplicates and setting a maximum time delay of 18hours and a global correlation coefficient of 0.9 and visualized using Cytoscape. Gene lists were functionally annotated with KEGG and Reactome pathways (adjusted p value<0.05) using the clusterProfiler and ReactomePA packages74.
Data has been deposited into GEO under accession number GSE75028.
The authors acknowledge Kenneth Sabir and Ruth Pidsley for reviewing the manuscript; the Garvan Foundation and the Peter Wills Bioinformatics Facility for providing facilities and Agilent Technologies for providing RNAseq kits. BG acknowledges constructive feedback from Aaron Statham, Mark Pinese, Nenad Bartonicek, Jesper Maag and Quek Xiucheng. B.G. is supported by Cancer Institute NSW Early Career Fellowship 13/ECF/1-45.
B.G. wrote the manuscript, assisted study conception, performed the analyses and library preparations (assisted by D.K.). M.D. conceived the study and assisted writing the manuscript. B.S. performed de-novo assembly and P.C. deconvolution, designed the web portal, assisted with figure composition and reviewed the manuscript. S.C., F.G. and D.K. performed lab work and reviewed the manuscript. A.P. provided design input, biological samples and facilities.
The authors declare that they have no competing interests.
Electronic supplementary material
Supplementary information accompanies this paper at doi:10.1038/s41598-017-06110-5
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.