A spatiotemporally resolved atlas of mRNA decay in the C. elegans embryo reveals differential regulation of mRNA stability across stages and cell types

Abstract

During embryonic development, cells undergo dynamic changes in gene expression that are required for appropriate cell fate specification. Although both transcription and mRNA degradation contribute to gene expression dynamics, patterns of mRNA decay are less well-understood. Here we directly measured spatiotemporally resolved mRNA decay rates transcriptome-wide throughout C. elegans embryogenesis by transcription inhibition followed by bulk and single-cell RNA-sequencing. This allowed us to calculate mRNA half-lives within specific cell types and developmental stages and identify differentially regulated mRNA decay throughout embryonic development. We identified transcript features that are correlated with mRNA stability and found that mRNA decay rates are associated with distinct peaks in gene expression over time. Moreover, we provide evidence that, on average, mRNA is more stable in the germline compared to in the soma and in later embryonic stages compared to in earlier stages. This work suggests that differential mRNA decay across cell states and time helps to shape developmental gene expression, and it provides a valuable resource for studies of mRNA turnover regulatory mechanisms.

Figure 1.. Transcription inhibition by actinomycin D allows mRNA decay measurements in C. elegans embryonic cells.

(A) A schematic representation of the approach to measure mRNA half-lives in the C. elegans embryo using transcription inhibition and bulk RNA-sequencing. (B) Distribution of mRNA half-lives across three biological replicates that met a moderate filtering strategy. Median half-life was 54 minutes (black line). 72 genes with half-lives greater than 200 minutes not shown. (C) Select gene ontology categories enriched among the top 15% stable and unstable transcripts. Background set of genes used was all genes that met our moderate mRNA half-life filtering metric. (D) A schematic representation of a gene with highly persistent expression and a gene with highly transient expression and the predicted overall stability of their transcripts. (E) Left. Expression over time of highly persistent (rpl-35, pab-1) and highly transient (zip-7, hlh-14) genes from a whole-embryo RNA-sequencing time series (Hashimshony et al. 2015). Right. The measured mRNA decay of the corresponding genes from our bulk RNA-sequencing data, with each point representing normalized transcript abundance from one of three biological replicates. (F) Box plots showing the mRNA half-life distributions of genes characterized as highly transient, transient, persistent, or highly persistent transcriptome-wide. Numbers to the left of the box plots are median half-lives within each group. Numbers above the box plots are the number of genes with half-lives greater than 150 minutes within each group. P-values comparing median half-lives were calculated using the Wilcoxon rank sum test.

Figure 2.. Transcript stability correlates with specific sequence features.

The relationship between transcript stability and different sequence features was examined using the top 15% stable and unstable transcripts identified in our bulk RNA-sequencing data. The box plots compare different features between all, stable, and unstable transcripts. Numbers to the left of the box plots are median values within each group. P-values comparing median values were calculated using the Wilcoxon rank sum test. Outliers not shown. (A) Distribution of codon optimality scores, as measured using tRNA adaptation index (tAI) values for each gene. Higher tAI value corresponds to greater codon optimality. (B) Distribution of the number of introns averaged across all splice isoforms of a gene. (C) Distribution of 3’ UTR length averaged across all 3’ UTR isoforms of a gene. (D) Distribution of percent GC content in the coding sequence of genes. For genes with multiple coding sequence isoforms, the longest isoform was used. (E) Distribution of percent GC content in the 3’ UTR of genes. For genes with multiple 3’ UTR isoforms, the longest isoform was used. (F) Distribution of coding sequence length averaged across all splice isoforms of a gene. (G) Linear regression was used to identify the percent of variation in mRNA half-lives explained by individual sequence features. (H) Multiple linear regression was used to identify the percent of variation in mRNA half-lives explained by combinations of different sequence features. (I) The de novo motif-finding program MEME (Bailey et al. 2015) identified three motifs enriched in the 3’ UTRs of the top 15% stable transcripts compared to the 3’ UTRs of the top 15% unstable transcripts. For genes with multiple 3’ UTR isoforms, the longest isoform was used. The best match of each de novo motif to known motifs in mammals (Ray et al. 2013) is noted.

Figure 3.. High transcript accumulation is associated with increased mRNA stability.

(A) A schematic representation of genes that accumulate high, medium, or low transcript levels and the expected range of transcription and mRNA decay rates that could contribute to such accumulation. (B) Examples of genes that peak in expression ~200 minutes after the four-cell stage with high, medium, and low transcript accumulation, from left to right. Expression data taken from a whole-embryo RNA-sequencing time series (Hashimshony et al. 2015). (C) Examples of genes that peak in expression ~350 minutes after the four-cell stage with high, medium, and low transcript accumulation, from left to right. Expression data taken from a whole-embryo RNA-sequencing time series (Hashimshony et al. 2015). (D) Box plots showing the mRNA half-life distributions of genes that peak in expression ~200 minutes after the four-cell stage to low, medium, and high transcript levels. (E) Box plots showing the mRNA half-life distributions of genes that peak in expression ~350 minutes after the four-cell stage with low, medium, and high transcript levels. (F) The de novo motif-finding program MEME (Bailey et al. 2015) identified three motifs enriched in the 3’ UTRs of genes with high transcript accumulation compared to the 3’ UTRs of genes with low transcript accumulation ~200 minutes after the four-cell stage. For genes with multiple 3’ UTR isoforms, the longest isoform was used. The best match of each de novo motif to known motifs in mammals (Ray et al. 2013) is noted. (G) The de novo motif-finding program MEME (Bailey et al. 2015) identified three motifs enriched in the 3’ UTRs of genes with high transcript accumulation compared to the 3’ UTRs of genes with low transcript accumulation ~350 minutes after the four-cell stage. For genes with multiple 3’ UTR isoforms, the longest isoform was used. The best match of each de novo motif to known motifs in mammals (Ray et al. 2013) is noted. (H) The gene expression patterns for three putative RNA-binding protein genes in C. elegans from a whole-embryo RNA-sequencing time series (Hashimshony et al. 2015). The genes are homologs of mammalian RNA-binding protein genes whose corresponding proteins are known to bind motifs similar to those discovered in (F) and (G). Numbers to the left of the box plots are median half-lives within each group. Numbers above the box plots are the number of genes with half-lives greater than 125 minutes (D) or 175 minutes (E) within each group. P-values comparing median half-lives were calculated using the Wilcoxon rank sum test.

Figure 4.. Single-cell RNA-sequencing allows measurement of mRNA half-lives at high resolution throughout C. elegans embryogenesis.

(A) A schematic representation of the approach to measure mRNA half-lives in the C. elegans embryo using transcription inhibition and single-cell RNA-sequencing. (B) UMAP projection of the integrated dataset of three biological replicates. Cells are colored by the age of the embryo from which a cell was produced, estimated from correlations to a whole-embryo RNA-sequencing time series (Hashimshony et al. 2015). Trajectories corresponding to major cell types are labeled. (C) Scatter plot comparing the mRNA half-lives calculated in pseudobulk from the single-cell data and the half-lives calculated from the bulk-cell data on a log-log scale. Spearman correlation coefficient = 0.76. Dashed line is the x = y line. (D) Box plots showing the stage-specific mRNA half-life distributions of genes within Early-, Middle-, and Late-stage cells (50–200, 200–350, and 350+ minutes after the four-cell embryo stage, respectively). (E) Box plots showing the cell type-specific mRNA half-life distributions of genes within muscle, germline, epidermis, neuron, and pharynx. Numbers above the box plots are the number of genes with half-lives greater than 100 minutes (D) or 150 minutes (E) within each group. P-values comparing median half-lives were calculated using the Wilcoxon rank sum test.

Figure 5.. Differential mRNA decay occurs throughout different developmental stages of C. elegans embryogenesis.

(A) Scatter plots comparing mRNA half-lives specific to Early-, Middle-, and Late-stage cells in all pairwise comparisons. Each point represents a gene. Blue points correspond to the top 5% of genes with faster mRNA decay in the later stage compared to in the earlier stage. Pink points correspond to the top 5% of genes with slower mRNA decay in the later stage compared to in the earlier stage. Dashed line is the x = y line. (B) Select gene ontology categories enriched among the top 5% of genes with faster mRNA decay in the later stage of embryogenesis compared to the earlier stage of embryogenesis. Background set of genes used in each comparison was shared genes we were able to calculate stage-specific mRNA half-lives for between the relevant stages. (C, D, E, F) Left. Median scaled expression of gene subset using data from a whole-embryo RNA-sequencing time series (Hashimshony et al. 2015). Pink shading spans the Early stage, light purple shading spans the Middle stage, and dark purple shading spans the late stage. Right. Plot displaying the change in mRNA half-lives from earlier to later stage for gene subset.

Figure 6.. mRNA degradation may contribute to transcription factor dynamics at both the RNA and protein levels.

(A) Box plots showing the mRNA half-life distributions of transcription factor genes and all other genes in Early-, Middle-, and Late-stage cells. (B) Grid displaying the percentage of transcription factors with transient or persistent mRNA expression and transient or persistent protein expression. The probability that RNA and protein dynamics are independent of one another was calculated using Fisher’s exact test. (C) Box plots showing the single-cell RNA-sequencing pseudobulk mRNA half-life distributions of transcription factors with transient mRNA and protein expression, transient mRNA and persistent protein expression, persistent mRNA and transient protein expression, and persistent mRNA and protein expression. (D, E) Sublineages with coloring representing reporter GFP expression from a single-cell transcription factor protein expression atlas of the C. elegans embryo (Ma et al. 2021). Protein and mRNA dynamics are characterized below, along with pseudobulk measured mRNA half-life. Numbers to the left of the box plots are median half-lives within each group. Numbers above the box plots in (A) are the number of genes with half-lives greater than 100 minutes within each group. P-values comparing median half-lives were calculated using the Wilcoxon rank sum test.

Figure 7.. mRNA stability is correlated with cell type-specific functions.

(A) Box plots showing the mRNA half-life distributions of cell type-specific genes within muscle, germline, epidermis, neuron, and pharynx cells. (B) Box plots showing the mRNA half-life distributions of broadly expressed genes within muscle, germline, epidermis, neuron, and pharynx. (C) Box plots showing the muscle-specific mRNA half-life distributions of mesoderm/muscle fate-specifying transcription factor genes, muscle structure genes, and all other genes. (D, E) Scatter plots of the normalized transcript abundance of the mesoderm/muscle fate-specifying transcription factor genes ceh-34, tbx-2, ceh-20, eya-1, hlh-1 and the muscle structure genes deb-1, unc-78, unc-112, sgn-1, tmd-2 throughout a 40 minute transcription inhibition time course in muscle cells. Each point represents normalized transcript abundance from one of three biological replicates. (F) Box plots showing the neuron-specific mRNA half-life distributions of neural fate-specifying transcription factor genes, synapse-/dendrite-/axon-associated genes, and all other genes. (G, H) Scatter plots of the normalized transcript abundance of the neural fate-specifying transcription factor genes ceh-5, egl-46, unc-86, ceh-43, hlh-2 and the synapse-/dendrite-/axon-associated genes lnp-1, unc-37, mig-2, ced-10, unc-53 throughout a 40 minute transcription inhibition time course in neuronal cells. Each point represents normalized transcript abundance from one of three biological replicates. Numbers to the left of the box plots are median half-lives within each group. Numbers above box plots are the number of genes with half-lives greater than 150 minutes within each group. P-value comparing median half-lives was calculated using the Wilcoxon rank sum test.

Figure 8.. Differential mRNA decay occurs across germline and somatic cells in the C. elegans embryo.

(A) Box plots showing the cell type-specific mRNA half-life distributions of genes across the germline and soma. (B) Scatter plot comparing soma- and germline-specific mRNA half-lives to one another. Each point represents a gene. Pink points correspond to the top 10% of genes with longer mRNA decay in the germline compared to in the soma. Blue points correspond to the top 10% of genes with faster mRNA decay in the germline compared to in the soma. Dashed line is the x = y line. (C) A cartoon of box plots representing the expected mRNA half-life distribution of genes with less, constant, or greater expression over time if mRNA decay primarily contributes to gene expression. (D) A cartoon of box plots representing the expected mRNA half-life distribution of genes with less, constant, or greater expression over time if transcription primarily contributes to gene expression. (E) Box plots showing the germline-specific mRNA half-life distributions for genes that decrease or increase in expression over time in the germline from a fold-change of 1–2 or greater than 2. (F) Box plots showing the soma-specific mRNA half-life distributions for genes that decrease or increase in expression over time in the soma from a fold-change of 1–2 or greater than 2. (G) Box plots showing the germline-specific mRNA half-life distributions of mRNA 3’ UTR-binding genes, other mRNA-binding genes, and genes not annotated as RNA-binding. (H) Box plots showing the soma-specific mRNA half-life distributions of mRNA 3’ UTR-binding genes, other mRNA-binding genes, and genes not annotated as RNA-binding. (I) Plot displaying the soma- and germline-line specific mRNA half-lives for mRNA-binding genes whose decay is more rapid in the germline than in the soma. Numbers to the left of the box plots are median half-lives within each group. Numbers above the box plots are the number of genes with half-lives greater than 150 minutes (A, G, H) or 325 minutes (E, F) within each group. P-value comparing median half-lives was calculated using the Wilcoxon rank sum test.