Owens ND, Blitz IL, Lane MA, Patrushev I, Overton JD, Gilchrist MJ, Cho KW, Khokha MK.

R01 GM099149 NIGMS NIH HHS , P40 OD010997 NIH HHS , R01 HD073179 NICHD NIH HHS , 5R01GM099149 NIGMS NIH HHS , Medical Research Council , Cancer Research UK

	Data generation, mRNA content of the embryo and gene expression dynamics, see also Fig S1, S2 A – Experimental design and sample collection B – Left: eef1a1o and ERCC-00171 abundance in relative normalization (TPM), Clutch A polyA+. RNA standard decreases with time. Right: absolute normalization of left. Red line indicates ERCC-00171 transcripts added based on the manufacturer’s datasheet. The RNA standard is constant at the correct abundance. C – Total mRNA (PolyA+ and rdRNA) in the embryo in nanograms with time. Gray region and line respectively mark Gaussian process 95% confidence interval (CI) and median of Clutch A & B PolyA+ data. D – Dynamics of gene expression of PolyA+ RNA in transcripts per embryo. Circles mark Clutch A PolyA+, triangles mark Clutch B PolyA+. The horizontal axis marks time (bottom) and NF developmental stage (top). Shaded regions mark Gaussian process 95% CIs of the data.
	Clutch A vs Clutch B PolyA+ Differential Expression A – Clutch A vs Clutch B histogram of differential expression scores as measured by the Bayesian Information Criterion (BIC), larger scores indicate greater differential expression (Supplemental Experimental Procedures). Three regions marked: 1) no differential expression BIC ≤ 0, 2) differential expression small effect 0 < BIC ≤ 60, 3) differential expression large effect BIC > 60. See B for explanation of thresholds. B – Differential expression effect size as measured by the log10 mean overlap between Clutch A and Clutch B Gaussian process models (Supplemental Experimental Procedures). Mean overlap decreases with increasing BIC. BIC = 0, 60 and mean overlaps of 60%, 10% and 2.5% marked. At BIC = 60 approximately all genes have less than 10% overlap, and all genes with less than 2.5% overlap have BIC > 60. C – Differential expression examples with decreasing BIC (top right). Genes on the boundary for strong differential expression have highly correlated expression profiles in A and B. Vertical lines mark convergences (red) or divergences (black). D – Left: Examples of convergence (dhx32) and divergence (Xetro.K05169). Right: Histogram of convergence and divergence times.
	Isoform Differential Expression in Clutch A, see also Fig S5 A – Isoform switching event examples: maternal to maternal; maternal to zygotic; and zygotic to zygotic switching events (left). Histogram of isoform switching events by category and time (right). B – bicd2 isoforms switch in PolyA+ data (circles), and do not switch in rdRNA data (squares) indicating differential polyadenylation. PolyA+ and rdRNA abundances agree within uncertainty of the absolute normalization (Supplemental Experimental Procedures). C – Normalized read depth on bicd2 locus in Clutch A polyA+ between 0–4.5 hpf. Note temporal dynamics shared by the final three bicd2(1) exons. Heatmap and depth of pile ups corrected for changing total mRNA. D – Spatial expression of three genes with isoforms with differential dynamics. All isoforms show different temporal and spatial expression domains. Grey lines and points mark total expression for each gene (sum of the red and blue isoforms).
	Temporal Synexpression in Clutch A PolyA+, see also Fig S6 A – Temporal map of the transcriptome: an enumeration of all gene expression dynamics in the embryo. Heat maps display all genes normalized by maximum expression and ordered by similarity. Inset (right) expands on genes transiently expressed with early onset. Vertical bars mark Gene Ontology (GO) enrichment (colors correspond to GO terms), numbers appended to GO bars indicate the number of genes of given category. All reported GO blocks are enriched with p < 2×10−4 (Fisher Exact) for entire block. S1, S2 and V1, V2 mark the location of somite and vision temporal synexpression genes respectively in B. “Transcription Factors *” labels an annotation of transcription factors separate to GO (Supplemental Experimental Methods). B – Somite (left) and Vision (right) synexpression groups.
	Transcript copy numbers and kinetics per embryo in Clutch A A – Transcript copy number histogram in polyA+ sequencing at 0.0 (egg), 4.5 hpf (early stage 9), 10 hpf (stage 12.5, end of gastrulation), 66 hpf (tadpole stage 42). Detection limits give number of transcripts required to generate a single read averaged over all transcript lengths (Fig S2E). Transcript copy number distributions transition smoothly between these timepoints (data not shown). B – Histogram of maximum rates of increase in transcripts/hour between any consecutive measurement points calculated from Gaussian process medians for each gene. Blue – PolyA+ 0–23.5 hpf; Red – rdRNA 0–23.5 hpf; Green – PolyA+ 24–66 hpf. Vertical lines mark the gene with max. rate of increase for each category. C – As B but kb/hour. D – Maximum rate of increase in transcripts/hour (left) and kb/hour (right) for 0–3.5 hpf and 4–23.5 hpf time intervals. PolyA+ and rdRNA distributions are discrepant for 0–3.5 hpf reflecting polyadenylation of maternal RNAs. Distributions are identical over 4–23.5 hpf.
	Accumulation kinetics of pri-mir427 per allele in Clutch A rdRNA, see also Fig S7 A – pri-mir427 expression over first 12 hpf in rdRNA data and in situ hybridizations. pri-mir427 show ubiquitous expression matching RNA-seq temporal data. Examples of most abundant transcripts shown for comparison. Inset: first detection of pri-mir427 above detection limit at 2.0 hpf (8–16 cell transition). B – pri-mir427 abundance in kb per cell. Line marks Gaussian process median, error bar is 95% CI (Supplemental Experimental Procedures). Peak per cell abundance occurs after 11th division at 4.4 hpf. C – pri-mir427 accumulation rate in kb/min/allele with time. Line is median of differential of Gaussian Process, error bar gives 95% CI (Supplemental Experimental Procedures). Peak accumulation rate is achieved at 4.2 hpf just after 10th cell division. D – Organization of pri-mir427 locus on X. tropicalis Chromosome 3 (X. tropicalis v8 genome, Supplemental Experimental Procedures). Clusters of mir427 hairpins appear in tandem with a repeated sequence arranged symmetrically on opposite strands around a gap in the genome assembly.
	Characteristic Timescales of Gene Expression in Clutch A PolyA+ A – Distribution of max. rate of increase of all transcripts (gray) (as Fig 5), transcription factors (blue), and ribosomal proteins rps* and rpl* (red). B – Theoretical covariance (Matérn covariance function, Supplemental Experimental Procedures) between two measurements by time interval for example genes in C. C – Histogram of characteristic timescale over all genes in in Clutch A polyA+. Top row annotates example genes (timescale inset), fast/short timescales are on the left slow/long timescales on the right, Gaussian process 95% CI (gray) and median (blue) marked. Bars below give GO enrichment (calculated as Fig 2, Fig S6A) of histogram, numbers indicate total genes with given category. All enrichments have p < 2.6×10−5 (Table S4). “Transcription Factors *” labels an annotation of transcription factors separate to GO (Supplemental Experimental Methods).
	Fig S1 – RNA-seq data quality, related to Fig 1, Supplemental Experimental Procedures (A) Spearman correlation between all pairs of all samples for relative normalized TPM abundances. Neighboring samples are as correlated as replicates. (B) Scatter plots of all TPM abundances for all samples. Color indicates density of points. Horizontal and vertical strips indicate measurements only present in one sample. (C) RNA standard performance. Scatter plots compare RNA standard TPM abundances between all pairs of datasets. Color indicates spike-in species. Red boxes mark ArrayControl spikes present only in Clutch B PolyA+. Although, we find small numbers of reads that align uniquely to ArrayControl sequences when ArrayControl spikes are not present in the sample. We note that this is not true for ERCC spikes, we never find ERCC reads when ERCC spikes are not present (data not shown). ERCC PolyA+ and rdRNA comparisons lie off the diagonal indicating poor PolyA+ performance. ERCC-00116 has the worst PolyA+ performance and is not used in absolute normalization. (D) Identification of poorly performing libraries. Normalized leave-one-out score (Supplemental Experimental Procedures), higher scores indicate a poorly performing sample. Samples adjacent to a bad sample can also receive high scores. Red and blue samples in Clutch A PolyA+ are selected for repeat library construction and sequencing (E). Red samples excluded from all analysis. (E) Technical replicates in Clutch A PolyA+ of poorly performing samples reveal that the poor performance is highly reproducible indicating that performance issues are not due to library construction and sequencing. Red boxes marked those discarded, black boxes marked those retained. (F) Projections onto the first two principle components for Clutch A PolyA+ (left), Clutch B PolyA+ (center), and Clutch A rdRNA (right). Note strong temporal correlation between samples, major changes in the transcriptome are visible.
	Fig S2 – Statistics of Absolute Normalization, related to Fig 1, Supplemental Experimental Procedures (A) Absolute normalization of RNA standards. Left: relative normalization in TPM on linear and log scale for the three datasets. The trend shared by all spikes is clearly visible and differs between PolyA+ and rdRNA samples. Right: absolute normalization in transcripts per embryo on log and linear scale. Trend in spikes is removed, the sample noise of the spikes is retained and does not contaminate the absolute normalization of native transcripts. (B) Absolute normalization factors. Circles mark per sample factors, lines mark Gaussian Process smoothed factors. Left: Estimates of 𝛽𝑗 (Supplemental Experimental Procedures), for three datasets. Center-left: per dataset correction factors e−𝛽𝑗, note discrepancy in magnitude but not in trend between datasets. Center-right: correction for PolyA+ bias, and averaging Clutch A and Clutch B spike levels results in excellent agreement between the correction factors. Right: Smoothed correction factors calculated with non-stationary Gaussian Process (Supplemental Experimental Procedures). (C) Absolute normalization consistency and accuracy. Scatter plots of ERCC standards spiked vs calculated for Clutch A PolyA+, Clutch A rdRNA, Clutch B PolyA+ and Clutch A/Clutch B PolyA+ combined. Points are averaged over all samples in time to remove spikein sample noise. ArrayControl spikes are given with fold change errors for Clutch B PolyA+, fold changes range from 1.11-1.25. ArrayControl spikes EC02 and EC12 have are spiked in with the same copy number and are marked with the same vertical line.Combined Clutch A/Clutch B gives 95% confidence interval for a linear model with Gaussan noise, with slope unity and with standard deviation 𝜎𝑠 = 0.25. Note that residuals do not depend on expression level. (D) Absolute normalization uncertainty model. Left: model as described in C (far-right), with 25%, 50%, 75%, 95% confidence intervals for 106 transcripts. Right: Propagation of uncertainty model to absolute normalization of eef1a1o Clutch A polyA+, again 25%, 50% 75% and 95% of confidence intervals for the true number of eef1a1o transcripts/embryo marked. (E) RNA-seq detection limits – Transcript abundance required to produce a single read with time, detection limits increase as RNA content in the embryo increases. Top row: Detection limits in kb, points give the per sample detection limits, lines and shaded area give Gaussian Process median and 95% confidence interval. Bottom: Detection limit in transcripts/embryo, calculated by averaging kb Gaussian Process distribution over all transcript lengths. (F) The depletion of non-ribosomal RNAs in rdRNA sequencing. Comparisons between Clutch A PolyA+ and rdRNA sequencing (top). Examples of extreme depletion: mixer and foxh1.2.
	Fig S3 – Spatial expression of selected genes by in situ, related to Fig 1 All panels give Clutch A PolyA+ (circles), Clutch B PolyA+ (triangles) and Gaussian Process median line and 95% confidence intervals. Numbers on in situ images indicate the Gaussian process median number of transcripts for the given stage. (A) Genes involved in muscle induction. (B) Genes with onset midway through the timecourse. (C) Genes with onset late in the timecourse.
	Fig S4 – Genes with constant expression – loading control candidates, related to Table S3 and Supplemental Experimental Procedures (A) Histogram of fold change envelope between maximum and minimum Gaussian Process median abundances in transcripts per embryo for all genes for Clutch A PolyA+ and rdRNA. Candidates for loading control with a fold change < 2 are highlighted. Histogram displays up to 8-fold envelope. (B) Traditional loading controls gapdh, odc1, eef1a1 have large fold changes between maximum and minimum expression. Fold change envelopes indicated. (C) Candidates for loading controls with maximum/minimum fold changes < 2. Full list in Table S3.
	Fig S5 – Differential Isoform Dynamics in Clutch A polyA+, related to Fig 3. (A) Explanation and detection of differential isoform dynamics. Left: isoform expression in transcripts/embryo. Center: Joint model of all isoforms assumed to have the same expression on z-normalized scale (see Supplemental Experimental Methods). Right: Models of all isoforms considered separately. If Bayesian Information Criterion (BIC) favors the joint model, the gene has differential abundance, whereas if the BIC favors the individual models a gene has differential dynamics. eef1a1 has differential abundance, nek2 and lima1 have differential dynamics. (B) Differential dynamics of groups of isoforms, BIC graph (left) between all lima1 isoforms. lima2(1) and lima2(2) have differential abundance, and they have differential dynamics with lima2(3) defining two groups of isoforms (right).
	Fig S6 – Temporal synexpression and gene ontology in Clutch A PolyA+, related to Fig 4. (A) Sliding window gene ontology enrichment. Bottom: Heatmap of all genes (Fig 4) arranged by similarity. Middle: The location of all transcripts factors, blue line marks the number of transcription factors per 500 gene window. Red line marks the expected number of transcripts for any 500 gene window. Top: -log(p) chi-squared test for the association between the window position and the transcription factor density. A region that exceeds –log(p) > 5 is considered for annotation (Supplemental Experimental Procedures). (B) Somite Synexpression. Examples from the somite temporal synexpresssion group. Ranks give similarity to ckb between 34-66 hpf (Supplemental Experimental Procedures). Image source given at bottom.
	Fig S7 – Early transcription events, related to Fig 6. (A) Stage 9 RNA pol II ChIP - Pile ups on mir427, eef1a1o, eef1a1, krt5.7 and cirbp loci of RNA pol II ChIP-seq reads (van Heeringen et al., 2014). See Supplemental Experimental Procedures for calculation of mir427 pile up. (B) Accumulation rates in kb/min/allele averaged over all cells in the embryo – for nodal3, nodal5, nodal6, sia1, sia2, klf17, admp and ventx3.1. Left column: kb/embryo. Center column: Detection of activation of each gene during the cleavage stages. Note square root vertical scale. Vertical gray dashed lines indicate cell divisions. Horizontal dashed line gives average detection limit of 970kb for the cleavage stages. Right column: Accumulation rate in kb/min/allele, rates account for exponential increase in cell numbers and are correct if each transcript were transcribed in all cells (Supplemental Experimental Procedures).
	ace2 (angiotensin I converting enzyme 2) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 43, lateral view, anterior right, dorsal up.
	btbd3 (BTB (POZ) domain containing 3) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 33, lateral view, anterior right, dorsal up.
	cpb1 (carboxypeptidase B1) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 41, lateral view, anterior right, dorsal up.
	ctrb1 (chymotrypsinogen B1) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 41, lateral view, anterior right, dorsal up.
	ctrl (chymotrypsin-like) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 42, lateral view, anterior right, dorsal up.
	g6pc2 (glucose-6-phosphatase, catalytic, 2) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 29, lateral view, anterior right, dorsal up.
	hbe1 (hemoglobin subunit epsilon 1) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 29, lateral view, anterior right, dorsal up.
	mep1a (meprin A subunit alpha) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 43, lateral view, anterior right, dorsal up.
	mep1b (meprin A subunit beta) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 43, lateral view, anterior right, dorsal up.
	thdl17(thyroid hormone down-regulated protein (gene 17) ) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 43, lateral view, anterior right, dorsal up.
	thdl18 ( thyroid hormone down-regulated protein (gene 18) ) gene expression in Xenopus tropicalis embryo, assayed via in situ hybridization, NF stage 40, lateral view, anterior right, dorsal up.

Aanes, Zebrafish mRNA sequencing deciphers novelties in transcriptome dynamics during maternal to zygotic transition. 2011, Pubmed

Aanes, Zebrafish mRNA sequencing deciphers novelties in transcriptome dynamics during maternal to zygotic transition. 2011, Pubmed
Äijö, Methods for time series analysis of RNA-seq data with application to human Th17 cell differentiation. 2014, Pubmed
Ardehali, Tracking rates of transcription and splicing in vivo. 2009, Pubmed
Baker, The External RNA Controls Consortium: a progress report. 2005, Pubmed
Bowes, Xenbase: gene expression and improved integration. 2010, Pubmed , Xenbase
Brennecke, Accounting for technical noise in single-cell RNA-seq experiments. 2013, Pubmed
Collart, High-resolution analysis of gene activity during the Xenopus mid-blastula transition. 2014, Pubmed , Xenbase
Crowder, A comparative analysis of spindle morphometrics across metazoans. 2015, Pubmed , Xenbase
Curran, Circadian genes are expressed during early development in Xenopus laevis. 2008, Pubmed , Xenbase
Curran, Circadian genes, xBmal1 and xNocturnin, modulate the timing and differentiation of somites in Xenopus laevis. 2014, Pubmed , Xenbase
Fakhro, Rare copy number variations in congenital heart disease patients identify unique genes in left-right patterning. 2011, Pubmed , Xenbase
Fang, Transcriptome analysis of early organogenesis in human embryos. 2010, Pubmed
Giraldez, Zebrafish MiR-430 promotes deadenylation and clearance of maternal mRNAs. 2006, Pubmed
Giraldez, MicroRNAs regulate brain morphogenesis in zebrafish. 2005, Pubmed
Heinonen, Detecting time periods of differential gene expression using Gaussian processes: an application to endothelial cells exposed to radiotherapy dose fraction. 2015, Pubmed
Hellsten, The genome of the Western clawed frog Xenopus tropicalis. 2010, Pubmed , Xenbase
Honkela, Model-based method for transcription factor target identification with limited data. 2010, Pubmed
Islam, Quantitative single-cell RNA-seq with unique molecular identifiers. 2014, Pubmed
Junker, Genome-wide RNA Tomography in the zebrafish embryo. 2014, Pubmed
Kalinka, Gene expression divergence recapitulates the developmental hourglass model. 2010, Pubmed
Karlebach, Modelling and analysis of gene regulatory networks. 2008, Pubmed
Karpinka, Xenbase, the Xenopus model organism database; new virtualized system, data types and genomes. 2015, Pubmed , Xenbase
Kimelman, The events of the midblastula transition in Xenopus are regulated by changes in the cell cycle. 1987, Pubmed , Xenbase
Klein, Genetic and genomic tools for Xenopus research: The NIH Xenopus initiative. 2002, Pubmed , Xenbase
Langfelder, WGCNA: an R package for weighted correlation network analysis. 2008, Pubmed
Lund, Deadenylation of maternal mRNAs mediated by miR-427 in Xenopus laevis embryos. 2009, Pubmed , Xenbase
Marguerat, Quantitative analysis of fission yeast transcriptomes and proteomes in proliferating and quiescent cells. 2012, Pubmed
Miller, Visualization of RNA synthesis on chromosomes. 1972, Pubmed
Morin, Sequencing and analysis of 10,967 full-length cDNA clones from Xenopus laevis and Xenopus tropicalis reveals post-tetraploidization transcriptome remodeling. 2006, Pubmed , Xenbase
Newport, A major developmental transition in early Xenopus embryos: II. Control of the onset of transcription. 1982, Pubmed , Xenbase
Niehrs, Synexpression groups in eukaryotes. 1999, Pubmed , Xenbase
Paranjpe, A genome-wide survey of maternal and embryonic transcripts during Xenopus tropicalis development. 2013, Pubmed , Xenbase
Pourquié, Vertebrate somitogenesis: a novel paradigm for animal segmentation? 2003, Pubmed
Qing, mRNA enrichment protocols determine the quantification characteristics of external RNA spike-in controls in RNA-Seq studies. 2013, Pubmed
Risso, Normalization of RNA-seq data using factor analysis of control genes or samples. 2014, Pubmed
Sagata, A study on the steady-state population of poly(A)+RNA during early development of Xenopus laevis. 1980, Pubmed , Xenbase
Satija, Spatial reconstruction of single-cell gene expression data. 2015, Pubmed
SEQC/MAQC-III Consortium, A comprehensive assessment of RNA-seq accuracy, reproducibility and information content by the Sequencing Quality Control Consortium. 2014, Pubmed
Skirkanich, An essential role for transcription before the MBT in Xenopus laevis. 2011, Pubmed , Xenbase
Spencer, A spatial and temporal map of C. elegans gene expression. 2011, Pubmed
Stegle, A robust Bayesian two-sample test for detecting intervals of differential gene expression in microarray time series. 2010, Pubmed
Stegle, Computational and analytical challenges in single-cell transcriptomics. 2015, Pubmed
Tan, RNA sequencing reveals a diverse and dynamic repertoire of the Xenopus tropicalis transcriptome over development. 2013, Pubmed , Xenbase
Tebbenkamp, The developmental transcriptome of the human brain: implications for neurodevelopmental disorders. 2014, Pubmed
Tomancak, Systematic determination of patterns of gene expression during Drosophila embryogenesis. 2002, Pubmed
Trapnell, Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. 2010, Pubmed
Trapnell, TopHat: discovering splice junctions with RNA-Seq. 2009, Pubmed
Tu, Quantitative developmental transcriptomes of the sea urchin Strongylocentrotus purpuratus. 2014, Pubmed
van Heeringen, Principles of nucleation of H3K27 methylation during embryonic development. 2014, Pubmed , Xenbase
Vesterlund, The zebrafish transcriptome during early development. 2011, Pubmed
wa Maina, Inference of RNA polymerase II transcription dynamics from chromatin immunoprecipitation time course data. 2014, Pubmed
Wang, A genome-wide study of gene activity reveals developmental signaling pathways in the preimplantation mouse embryo. 2004, Pubmed
Yanai, Mapping gene expression in two Xenopus species: evolutionary constraints and developmental flexibility. 2011, Pubmed , Xenbase
Yang, Beta-catenin/Tcf-regulated transcription prior to the midblastula transition. 2002, Pubmed , Xenbase
Zaidi, De novo mutations in histone-modifying genes in congenital heart disease. 2013, Pubmed