Dynamically regulated transcription factors are encoded by highly unstable mRNAs in the Drosophila larval brain

  1. Ilan Davis1
  1. 1Department of Biochemistry, University of Oxford, Oxford OX1 3QU, United Kingdom
  2. 2Mathematical Institute, University of Oxford, Oxford OX1 3LB, United Kingdom
  1. Corresponding author: ilan.davis{at}bioch.ox.ac.uk
 
Next Section

Abstract

The level of each RNA species depends on the balance between its rates of production and decay. Although previous studies have measured RNA decay across the genome in tissue culture and single-celled organisms, few experiments have been performed in intact complex tissues and organs. It is therefore unclear whether the determinants of RNA decay found in cultured cells are preserved in an intact tissue, and whether they differ between neighboring cell types and are regulated during development. To address these questions, we measured RNA synthesis and decay rates genome wide via metabolic labeling of whole cultured Drosophila larval brains using 4-thiouridine. Our analysis revealed that decay rates span a range of more than 100-fold, and that RNA stability is linked to gene function, with mRNAs encoding transcription factors being much less stable than mRNAs involved in core metabolic functions. Surprisingly, among transcription factor mRNAs there was a clear demarcation between more widely used transcription factors and those that are expressed only transiently during development. mRNAs encoding transient transcription factors are among the least stable in the brain. These mRNAs are characterized by epigenetic silencing in most cell types, as shown by their enrichment with the histone modification H3K27me3. Our data suggest the presence of an mRNA destabilizing mechanism targeted to these transiently expressed transcription factors to allow their levels to be regulated rapidly with high precision. Our study also demonstrates a general method for measuring mRNA transcription and decay rates in intact organs or tissues, offering insights into the role of mRNA stability in the regulation of complex developmental programs.

Keywords

Previous SectionNext Section

INTRODUCTION

Development of a multicellular organism requires intricate control of protein production in response to specific signals. Much of this control is exercised at the level of transcription, with signaling cascades leading to binding of transcription factors to recognition sites in the genome. Transcription factor binding precipitates production of a target mRNA, and eventually, translation of protein required in a specific cell at a specific time. Another class of gene regulation, known as post-transcriptional regulation, operates on the mRNA after its production. mRNAs can be regulated post-transcriptionally by several means: translation, intracellular localization, and decay, among others. The stability of an mRNA determines whether it is available for translation into protein and whether it reaches the region of the cell where it is required. Regulation of mRNA stability is a key component of gene regulation because the inappropriate presence of even a single molecule of mRNA could have drastic consequences. According to some estimates, one molecule of mRNA directs the synthesis of 6000 molecules of protein on average (Schwanhäusser et al. 2011). In combination with translational control, rapid degradation of unneeded RNAs can prevent the accumulation of proteins in the incorrect cell type or at the incorrect developmental stage.

Multicellular development involves many rapid transitions in cell types and finely tuned spatial boundaries. Such cell fate changes are often accompanied by rapid changes in gene expression. For example, reiterated, segmental organization of the early Drosophila embryo depends on the striped expression of pair-rule segmentation genes in alternate segment-wide stripes (Chipman 2020). Expression of these transcripts is highly dynamic (Hafen et al. 1984; Lim et al. 2018), with the mRNAs of fushi tarazu (ftz) and other pair-rule transcripts having half-lives of only 7 min (Edgar et al. 1986, 1989). Instability of Ftz protein appears also to be critical for correct patterning (Kellerman et al. 1990).

Later in development, the embryo undergoes even more complicated patterning to create tissues and organs, the most intricate of which is probably the brain. In Drosophila, neural identities are determined by a defined series of temporal transcription factors (Grosskortenhaus et al. 2005), and similar transcriptional cascades have been described in mammalian neural stem cells (Beckervordersandforth et al. 2015). For several switches, a rapid transition in expression of successive factors depends on transcriptional repression of an early factor by its successor (Doe 2017). However, it seems likely that robust switching between stem cell identities also depends on rapid mRNA degradation or translational repression. Regulation of RNA stability is likewise important for neural differentiation, and specific RNA-binding proteins stabilize or destabilize their targets during fate transition (Moon et al. 2018; Samuels et al. 2020).

Current measurements of RNA stability have largely been derived from single-celled organisms or cultured mammalian cells, and so do not directly address in vivo mechanisms of cell-state changes in developing tissues. Very few studies have estimated RNA stability in multicellular organisms (Burow et al. 2015; Sami et al. 2022), and we are unaware of any data describing RNA stability across an entire brain. In this paper, we estimate synthesis and decay rates in the Drosophila larval brain genome wide by analyzing the metabolic labeling kinetics of more than 7000 RNAs. We show that decay rates of different RNA species vary by orders of magnitude, with rates differing greatly across different gene functions. In particular, we found that mRNAs that encode transcription factors are particularly unstable, especially if their expression is restricted to a limited number of cell types. Many of these cell type-specific transcription factor RNAs are derived from genes that are targets of histone H3K27 trimethylation, a chromatin modification that is associated with transcriptional repression. Our results suggest a link between transcriptional and post-transcriptional repression during brain development.

Previous SectionNext Section

RESULTS

To study RNA dynamics genome-wide, we first devised a protocol to address the technical challenges of estimating decay rates from small tissue samples (see Materials and Methods). We found that RNA transcription and decay rates can be estimated from these samples using metabolic labeling, the steady-state assumption, and a single time point design (de Pretis et al. 2015). In this method, synthesis rates are estimated directly from a pulse of 4-thiouridine (4sU) RNA labeling, and decay rates are inferred by comparing the synthesis rates with the measured total RNA levels from the same sample. The steady-state assumption requires that total RNA levels do not change over the course of the assay. To evaluate this assumption, we collected RNA from brains either directly after dissection or after 60 min of incubation in our culture conditions. We found that the global pattern of RNA levels was highly correlated between the two time points (r2 = 0.94), supporting our use of the steady-state model (Supplemental Fig. S1A). A similar result was found when we compared brains incubated for 60 min in the presence of 4sU with brains harvested directly after dissection (Supplemental Fig. S1B), suggesting that incubation in 500 µM 4sU does not cause large gene expression changes over this time frame. We also examined a measure of RNA dynamics, the ratio of degradation/processing rate that can be approximated by the ratio of premature/mature RNA in total RNA samples (Furlan et al. 2020). Likewise, we did not observe large changes in RNA dynamics after 60 min of brain incubation in the presence or absence of 4sU (Supplemental Fig. S1C,D).

We next generated genome-wide estimates of RNA synthesis and decay rates in the larval brain. We added 4sU to third-instar (L3) larval brains cultured briefly ex vivo (Hailstone et al. 2020), and then isolated RNA from the brains after 20 min of labeling time (Fig. 1A). This time frame is short enough to resolve transcriptional dynamics for many RNAs (Mukherjee et al. 2017), yet long enough to allow precision in the labeling duration and to recover a reasonable amount of labeled RNA. RNAs containing 4sU were purified and RNA sequencing libraries were constructed from both the labeled 4sU+ RNA and the total RNA. 4sU+ RNA was rich in unspliced intronic regions relative to total RNA, as expected, because 4sU+ RNAs are either recently transcribed or in the process of transcription (Fig. 1B). The data resulting from these sequencing experiments were then used to estimate RNA synthesis and decay rates (de Pretis et al. 2015) and to explore the roles of RNA dynamics in the brain.

FIGURE 1.
View larger version:
FIGURE 1.

Measurement of RNA transcription and decay rates in larval brains. (A) A protocol for measurement of RNA transcription and decay rates in larval brains. (B) Estimation of the percentage of spliced and unspliced transcripts in the total and 4-thiouridine enriched (4sU+) libraries. (C) Histogram of estimated synthesis and decay rates among all genes. The standard deviation (σ) is used to describe the variation between rates across all genes. The rates were calculated in units TPM/min (transcripts per million per min) and 1/min. For synthesis rates, σ = 0.79, and for decay, σ = 0.44, for log10-transformed rates. (D) Half-lives of all RNAs in the larval brain. Half-life estimates were capped at 1000 min, due to estimation uncertainty of very long half-lives (see Materials and Methods). (E) Box plot showing the synthesis and decay rates of coding mRNAs compared with noncoding RNAs. (*) P < 0.05, (***) P < 10 × 10−40. (F) Half-lives of noncoding RNAs. Named noncoding RNAs with long half-lives are indicated.

We first used our genome-wide data set to assess the contribution of transcription and decay rates to the total RNA levels in the brain. Synthesis rates span a wider range than decay rates (σ = 0.79 vs. σ = 0.44 for the log-transformed rates) (Fig. 1C). We transformed the decay rates to half-lives and found that the median RNA half-life was 82 min (Fig. 1D), similar to recent estimates in mammalian cells (Schwanhäusser et al. 2011; Mukherjee et al. 2017). However, the distribution is highly skewed with many short half-lives and a long tail of more stable RNAs. Using our method, we were not able to resolve half-lives greater than approximately 1000 min (see Materials and Methods). Therefore, we grouped genes with half-lives greater than 1000 min into a single bin labeled ≥1000 min (Fig. 1D). Noncoding RNAs had much higher decay rates than coding transcripts, but only slightly lower synthesis rates (Fig. 1E). Of note, we found that many highly stable noncoding RNAs have known functions; for example, 7SL and RNase MRP have established metabolic roles, and roX1 and cherub have important developmental functions (Fig. 1F; Meller et al. 1997; Landskron et al. 2018). We then used this data set as a starting point to identify features linked to RNA stability in the developing brain.

We hypothesized that RNAs that are localized to neurites need to be stable in order to survive their long transit from the nucleus to the periphery. To address this hypothesis, we defined neurite-localized mRNAs based on homology with mammalian mRNAs that were found to be enriched in neurites across several studies (Kügelgen and Chekulaeva 2020) and compared their decay rates with those of other RNAs. We found that neurite-localized mRNAs have lower decay rates than other mRNAs in the genome (Supplemental Fig. S2). They also have higher synthesis rates and higher total RNA levels. We conclude that mRNAs that are present at the distal periphery of neuronal cytoplasmic projects are more highly synthesized, stable, and have higher total RNA levels. We propose that the low decay rates of these mRNAs helps them survive their long journey to the tips of neurites.

To investigate the relationship between the function of mRNAs and their stability, we grouped genes into functional categories using GO Slim annotations and determined the stability profiles of the gene groups (Fig. 2). We found the general principle that mRNA stability was intimately linked to gene function. Genes encoding transcription factors produced RNAs with much lower stability than average (median around 28th percentile). In contrast, genes encoding structural molecules, many of which are ribosomal protein genes, had very high stability (median around 79th percentile). In between these two outliers, gene groups had less severe bias, but many were nonetheless intriguing. RNAs encoding receptors had low median stability whereas RNAs encoding binders of receptors had a bimodal distribution, with some having relatively low and others relatively high stability (Fig. 2). There is only a weak global correlation between RNA decay rate and total RNA level (r2 = 0.04), and the correlation is absent when examining the genes encoding transcription factors alone (r2 = 0.00) (Supplemental Fig. S3A). Some highly expressed mRNAs are also unstable. For example, many unstable mRNAs encoding receptors are also highly expressed (Supplemental Fig. S3B). Together these results argue for tight regulation of RNA stability and coupling of RNA degradation to the functional roles of the mRNA.

FIGURE 2.
View larger version:
FIGURE 2.

Different functional classes of mRNAs have distinct decay rates. Stability percentile of gene ontology (GO) Slim molecular function categories for mRNAs. The GO Slim categories were taken from the Flybase functional ribbon annotations. The histograms are displayed as kernel density estimates.

Brain tissue consists of many different cell types. We tested whether RNAs that are highly regulated and expressed in only certain cell types would be more unstable than other RNAs that are more widely expressed in the brain. We found that in general, this is not true, but it is true for the subclass of RNAs that encode transcription factors. To classify RNAs as cell type-specific or nonspecific, we made use of a single-cell RNA-seq study of the larval brain (Corrales et al. 2022). We classified RNAs as cell type-specific (CTS) if they were both statistically enriched in a given cell type and at least twofold more abundant in that cell type relative to the rest of the brain. We note that this comparison is an approximation because our RNA stability data were collected from whole brains rather than specific cell types. RNAs enriched in certain cell types such as neuroepithelial cells were generally unstable (Fig. 3A). In contrast, RNAs enriched in hemocytes, glia, and several classes of mature neurons tended to be stable (Fig. 3A). We repeated this analysis with RNAs enriched in other cell types of the brain that had been further subdivided into more specific cell types using the same data set (Dillon et al. 2022). Analysis of these more specific cell types also revealed interesting differences between them. For example, RNAs enriched in type I and quiescent neuroblasts tended to be stable, whereas those enriched in type II neuroblasts tended to be unstable (Supplemental Fig. S4). When viewed as a group, CTS RNAs are more stable than average if RNAs encoding transcription factors are excluded (Fig. 3B).

FIGURE 3.
View larger version:
FIGURE 3.

Cell type-specific transcription factor RNAs are rapidly degraded. (A) The stabilities of RNAs (including both mRNAs and ncRNAs) that are enriched in specific cell types are shown. mRNAs encoding transcription factors (TFs) are indicated in blue. The asterisks at the right of the plot indicate that the stabilities for the indicated group of cell type-specific RNAs are different than other RNAs in the data set (Mann–Whitney U-test). The heat maps at the right display the fraction of each category that is composed of TFs and the significance of the enrichment. The enriched RNAs were extracted from a single-cell RNA-seq data set derived from larval brains (Corrales et al. 2022). (N) Neurons, (KCs) Kenyon cells, (NPCs) neural progenitor cells, (RGL) ring gland, (Hemos) hemocytes, (NE-OL) neuroepithelia/optic lobe, (NE Tch) neuroepithelia/trachea, (Immat) immature, (Chol) cholinergic, (Glut) glutamatergic, (Undet) undetermined, (Pept) peptidergic, (Octop) octopaminergic, (GABA) gabaergic, (Serot) serotoninergic, (Dopa) dopaminergic. (B) Stability of cell type-specific RNAs versus nonspecific RNAs for transcription factor genes (TFs) and other genes (non-TFs). The numbers in parentheses under the box plots indicate the number of genes in each category subdivided further into the ubiquitous and cell type-specific groups. (*) P < 0.05, (**) P < 10 × 10−10, (***) P < 10 × 10−40.

To explore this effect further, we analyzed the functions of RNAs in the top and bottom 10% of stability for CTS RNAs. Stable CTS RNAs were enriched for metabolic processes and regulation of transporters (Supplemental Fig. S5). In contrast, unstable CTS RNAs were most enriched for transcriptional regulators. Strikingly, CTS RNAs that encode transcription factors were more unstable than other transcription factor RNAs (16th percentile of stability vs. 28th percentile of stability) (Fig. 3B). Among the various cell types, we found that RNAs specific to immature neurons, neural progenitor cells, neuroepithelia, and cholinergic neurons were enriched with mRNAs encoding transcription factors (Fig. 3A).

To put our results into context, we examined the RNA stability and cell type enrichment for another class of regulators, the RNA-binding proteins (RBPs), and the subclass of these that is known to bind mRNAs, the mRBPs. RBPs and mRBPs were only enriched in RNAs specific to neural progenitor cells (Supplemental Fig. S6A). Cell type-specific mRBP mRNAs, but not cell type-specific RBP mRNAs, have lower stability than other cell type-specific mRNAs (Supplemental Fig. S6B). These analyses show that cell type specificity is linked to RNA stability, but only for certain classes of genes.

We next asked what features could explain the low stability of cell type-specific transcription factor RNAs. These RNAs had no significant enrichment for known RNA-binding protein motifs (see Materials and Methods; Krismer et al. 2020). However, we found that they are targets of epigenetic silencing. Silenced genes are marked by regions of histone H3K27 trimethylation, which often extend through both the gene body and its regulatory regions. These genes are often referred to as Polycomb group (PcG) target genes (Simon and Kingston 2013; Kassis et al. 2017). To identify potential PcG target genes, we examined H3K27 trimethylation sites in the fly genome collected from L3 larvae as part of the modENCODE project (Roy et al. 2010). A total of 339 genes for which we have determined RNA decay rates correspond to genes that have methylation sites within 1 kb of their transcription start site (upstream), within 1 kb of their transcription stop site (downstream) and within the gene body, the region corresponding to the transcript itself (including both introns and exons of the gene). Because these genes are predicted PcG targets, rather than experimentally confirmed targets, we refer to them as extended H3K27me3 genes. These genes have exceptionally low stability when compared with genes without H3K27 trimethylation, and other classes of genes with methylation covering only some of these regions also had lower stability than average (Fig. 4). This observation suggests a link between transcriptional silencing and high RNA decay rates.

FIGURE 4.
View larger version:
FIGURE 4.

Unstable transcription factor RNAs are targets of H3K27 trimethylation. (A) The RNA stability percentile of genes with H3K27me3 binding sites is shown. The asterisks at the right of the plot indicate that the stabilities for the indicated group of cell type-specific RNAs are different than other RNAs in the data set (Mann–Whitney U-test). The heat maps at the right display the fraction of each category that is composed of transcription factors (TFs) and the significance of the enrichment. The H3K27me3 binding sites were derived from ChIP-seq of whole L3 larvae that were collected as part of the modENCODE project (Celniker et al. 2009). Genes were classed as having a binding site within 1 kb of their transcription start site (upstream [A]), within 1 kb of their transcription termination site (downstream [C]), or within the gene body itself, including both introns and exons (gene body [B]). We assigned genes to a methylation category based on the combination of regions that had a methylation binding site. Each gene is assigned to only one category. Numbers in parentheses indicate the number of genes in each category. We excluded the upstream–downstream (AC) category because of scarcity; only 10 genes fell into this category. (*) P < 0.05, (**) P < 10 × 10−10, (***) P < 10 × 10−40.

We wondered whether the relationships between RNA decay rates and cell type specificity, H3K27 trimethylation, and transcription factor status could be disentangled from each other. These gene groups are highly overlapping, with RNAs encoding transcription factors making up 49% of extended H3K27me3 genes and cell type-specific transcription factors making up 13% of extended H3K27me3 genes (Fig. 5A). These overlaps are highly significant. Transcription factor genes are enriched for extended H3K27me3 (odds ratio 18.3, CI = 14.4–23.2, P < 0.001), and CTS transcription factor RNAs are also enriched for extended H3K27me3 (odds ratio 17.8, CI = 11.8–26.9, P < 0.001) (Fig. 5A). The overlap among these gene groups complicates any simple interpretations that do not consider the relationship between them.

FIGURE 5.
View larger version:
FIGURE 5.

Relationships between RNA decay and gene attributes. The relationship between RNA decay rates and gene membership in the following categories is shown: transcription factor (TF), extended H3K27me3 genes (me3), and cell type-specific genes (CTS). H3K27me3 target genes were defined as those that have H3K27me3 binding sites upstream of, downstream from, and covering the gene body, as specified in Figure 4. (A) The overlap between each gene group is shown. CTS TFs are cell type-specific transcription factor genes. The P-value for the overlap enrichment is shown. Note that the overlap between CTS TF and TF and the overlap between TF and TF is not shown because these categories have 100% overlap by definition. The numbers in parentheses indicate the number of genes in each group. (B) The log transformed and winsorized decay rates are shown. (C) The equation for the linear regression model is shown. The model predicts the log10 decay rate using gene membership in each of the three groups (TF, me3, and CTS) as features. The model was built either with or without interactions between the features. β0–β7 indicate the coefficients, with only β0–β3 determined for the model without interactions. The model without interactions can be derived from this equation simply by setting the interaction coefficients (β4–β7) to 0. (D) A multiple linear regression model was built to predict RNA decay rate as a function of TF, CTS, and H3K27me3 status. The blue diamonds indicate the values determined using all data, whereas the gray box plots show the range of values obtained from subsampling the data set with bootstrapping. Bootstrapping was performed by subsampling each nonoverlapping gene group with a unique combination of factors using a sample size equivalent to the smallest group size (n = 39). Ten-thousand bootstrap samples were taken. (Note that the eight nonoverlapping gene groups are shown in F.) (E) P-values obtained from the linear regression model in D. The P-values are the −log10 value of the original P-values. The P-values derived from the model using the entire data set are written above each coefficient in blue. The P-values corresponding to β0 using the entire data set are 0; thus, the −log10 value is undefined but could be interpreted as an infinitely large number and is indicated by “inf.” (F) The decay rates of all eight possible nonoverlapping groups derived from three possible explanatory features is shown. Each group is indicated by a vector (TF, me3, or CTS), where 1 = present and 0 = not present. (G) For each of the three possible explanatory features (TF, me3, and CTS), the mean change between a group that is positive for the feature is compared with the group that is negative for the feature but is identical with respect to the other features. The overlaid numbers indicate the number of genes in each group. An asterisk marks the comparisons that show the effect of adding either TF or me3 to a group that already has the other feature (i.e., adding TF to a group that is already positive for me3 or vice versa).

To determine which features(s) have the largest correlation with RNA decay rates, we performed multiple linear regression to model the decay rates as a function of gene membership in each of the three groups (TF = encodes a transcription factor, me3 = extended H3K27 trimethylation, and CTS = cell type specificity). We found that these features differ in their effects, and that interactions between the features change the interpretation of the effects. Before regression, we log transformed the decay rates and winsorized the top and bottom 1% of the data to reduce the effect of outliers (Fig. 5B). Next, we modeled the decay rates as a function of RNA membership in each gene group, either with or without interactions between the features (Fig. 5C–E). The obtained β coefficients indicate the strength and directionality of the effect, whereas the P-values indicate the probability that the coefficient is zero. The coefficients β0, β1, β2, and β3 are present in both models, but they can differ. For example, β1 from the noninteraction model indicates the effect of TF without considering that the TF category overlaps the CTS category, but the model with interactions takes this overlap into account, with the coefficient β5 describing the effect of TF and CTS combined. The model with interactions shows that TF has an effect both independently of other categories and in combination with CTS (both β1 and β5 are positive) (Fig. 5D, blue diamonds).

In order to focus our analysis on genes that belong to our groups of interest, we rebuilt the model with 10,000 bootstrapped data sets sampled from our data (Fig. 5D,E). This procedure will give each of our gene sets of interests the same importance in the model, rather than importance proportional to their frequency in the original data set. For the bootstrapped models, we randomly sampled each independent gene group (groups shown in Fig. 5F) at the same depth (39 genes per group per iteration, which corresponds to the minimum size of all gene groups). The bootstrapped models largely agree with the models built using the whole data set, and the range of values obtained in the bootstrapped models reinforce confidence in the effects of some features. The coefficients obtained from bootstrapping differ slightly from the coefficients derived from the whole data set in the model without interactions (Fig. 5D), showing, for example, that the effect of β1 (TF) is higher when the analysis is focused on our gene groups of interest than when using the whole data set. In contrast, the coefficients are nearly identical between the bootstrapped models and the model, considering the whole data set in the case of the model with interactions; this is likely because many of the independent gene groups that are used to determine the interaction coefficients are small and highly or completely sampled in the bootstrapping procedure. We found that β2 (me3) and β5 (TF and CTS) are positive in nearly all the bootstrapped models, reinforcing our confidence that both H3K27 trimethylation and encoding a transcription factor in combination with cell type specificity are associated with higher RNA decay rates (Fig. 5D, right). β1 is positive and β3 is negative in most bootstrapped models, showing that encoding a transcription factor and cell type specificity are associated with high and low RNA decay rates, respectively (Fig. 5D, right). Overall, our linear models suggest that both H3K27 trimethylation and encoding a TF are predictors of high RNA decay rates, as is the combination of the features cell type-specific and encoding a transcription factor.

To address the same question in a different way, we examined the decay rates of each independent gene group in which the combination of features is constant (Fig. 5F). We found that extended H3K27 trimethylation is associated with a similar increase in decay rates regardless of the other features it is combined with, whereas the other groups show more variable behavior. Each group is denoted by a vector that specifies its combination of features, with zero indicating that the feature is not present and one indicating that the feature is present (see the legend for Fig. 5). For example, (1,0,0) identifies genes that encode transcription factors but do not have any other features. This gene set does not overlap with the set denoted by (1,1,0), which specifies genes that both encode TFs and have extended H3K27 trimethylation. Direct comparison of each independent group gives similar insights with those obtained via linear regression. The CTS feature alone is not a consistent indicator of high RNA decay rate and in some comparisons is associated with a lower decay rate (Fig. 5G). In contrast, extended H3K27 trimethylation is a consistent indicator of high decay rate across all comparisons, whereas the effect of encoding a transcription factor is strongest in combination with the CTS feature (Fig. 5G).

Because of the strong and consistent correlation of H3K27 trimethylation with high RNA decay rates, we wondered whether there might be a functional link between PcG targeting and RNA decay. Such a link would ensure that transiently expressed RNAs could not wreak havoc in another cell type as a consequence of leaky transcription and that they could be cleared effectively before transitioning to another cell state. The most obvious explanation for the link is that convergent evolutionary selection processes have led to both PcG-driven transcriptional silencing and low RNA stability for these genes. However, the most parsimonious explanation for the link would be that the two are functionally related; that is, whatever processes serve to transcriptionally silence the gene also cause the locus-derived RNA to be rapidly degraded. To distinguish between these possibilities, we examined a data set that measured both nascent and total RNA levels after knockdown of components of the Polycomb repressive complex 1 (PRC1) in a Drosophila neural-derived cell line BG3 (Pherson et al. 2017). As expected, knockdown of the Ph subunit of PRC1 increased the levels of many RNAs at the nascent RNA level. These RNAs also increased in total RNA levels (Supplemental Fig. S7A–C). Nascent RNA sequencing captures the active transcriptional response to a change. It is expected that a change in transcriptional regulation of a gene would change its nascent RNA levels, but that the total RNA levels would mirror the changes in nascent RNA with a slight time delay depending on the half-life of the RNA in question. If a perturbation acts by altering RNA stability, then we would expect it to change the ratio between the nascent and the total RNA levels. We found that the ratio between the total and nascent RNA levels was slightly increased after Ph knockdown for genes that were defined as probable targets of PcG-driven silencing (PcG domain genes) (Supplemental Fig. S7D). This experiment suggests that Ph depletion primarily derepresses the transcription of PcG domain genes, but that the stability of the resulting RNAs remains similar after they are transcriptionally derepressed. Therefore, we failed to find evidence of a direct causative link between PcG targeting and RNA decay.

Finally, we examined the RNA stabilities for two groups of genes that are known to be spatiotemporally regulated in either the Drosophila embryo or larva. One group of genes is involved in neural fate determination, via specifying the developmental timing window in which a neuron is born. Another group of genes regulates anterior–posterior patterning in the Drosophila embryo. These groups of genes had low RNA stability, in agreement with our finding that rapid RNA decay is needed for RNAs that are required transiently in development (Fig. 6). We note that five of these genes (hb, Kr, tll, and slp1/2) are shared between the groups and have functions in both processes. The majority of these RNAs encode transcription factors (Fig. 6). Therefore, we suggest a model wherein mRNAs encoding transiently expressed and dynamically regulated transcription factors are selected to be unstable so that rapid developmental transitions can occur robustly.

FIGURE 6.
View larger version:
FIGURE 6.

RNAs that are highly regulated during development are rapidly degraded. The stability percentile of RNAs (both mRNAs and ncRNAs) that are annotated as being involved in neural fate determination or early patterning of the Drosophila embryo are shown. For neural fate determination, the RNAs are a combination of temporal transcription factor genes and other genes that undergo temporal regulation in neural precursor cells (Doe 2017). For pattern formation, the genes were taken from a curated list hosted on the Interactive Fly website and encompassing gene sets known as the gap genes and the pair-rule genes (Brody 1999). Most of these genes are transcription factors (TFs); the non-TFs are indicated. Genes that are in our group of extended H3K27me3 genes (me3) are also indicated.

Previous SectionNext Section

DISCUSSION

Here, we report a systematic analysis of RNA transcription and decay rates in the developing Drosophila brain. We accomplished this feat by adapting 4sU metabolic labeling and sequencing to small tissue samples. Although we collected RNA from 150 larval brains per replicate, we have shown that the experiment can be performed with as little as 5 µg of total RNA, which can be collected from as few as 25 brains (see Materials and Methods). Our study serves as proof of principle that RNA labeling experiments and kinetic data can be collected from complex tissue types, which are often more biologically relevant than monolayer cells in culture. There are many potential uses for our data set. For example, our analysis of decay rates identified many lncRNAs with higher than average RNA stability. Among these lncRNAs was the RNA cherub (half-life = 248 min), which was recently found to have an important function in brain tumor formation (Landskron et al. 2018). It will be interesting to examine the functions of other high stability lncRNAs, which are likely to have important roles.

Using these data, we showed that transcription rates span a larger range than decay rates, but decay rates are nonetheless biologically relevant. As shown in previous studies of RNA decay from tissue culture cells and single-celled organisms, RNA decay is tightly coupled to the function of the encoded protein (Miller et al. 2011; Schwanhäusser et al. 2011; Neymotin et al. 2014; Mukherjee et al. 2017). Interestingly, not all high expression RNAs are stable. For example, some RNAs encoding receptors are highly abundant but also unstable. Such dynamics might allow cells to respond to changes in conditions rapidly, even for abundant RNAs.

By comparing our data with single-cell RNA-seq data sets (Corrales et al. 2022; Dillon et al. 2022), we found that some cell types are defined by particularly unstable RNAs. Although it is true that our RNA stability measurements collected from whole brains do not directly correspond to the individual cell types profiled in these single-cell data sets, our measurements are likely to be similar to the actual cell type-specific values for RNAs that are highly enriched in a specific cell type. In other cases, our RNA stability measurements derived from whole brains should reflect the median stability across these different cell types, and this median measurement reveals differences between RNAs derived from genes of different functional classes.

Perhaps the most striking finding was the high turnover rate of RNAs encoding transcription factors. Although some past studies of RNA kinetics in tissue culture cells have observed lower stability for RNAs encoding transcription factors, our study takes this observation a step further. We found that the low stability status of these RNAs is linked to their tendency to be transiently expressed and developmentally regulated. These characteristics of the unstable transcription factor RNAs are demonstrated via their enrichment in certain transient cell populations, such as among immature neurons. Likewise, their high levels of H3K27 trimethylation in data derived from whole larvae indicates that these genes tend to be kept in a transcriptionally silent state in most cell types.

There are two possible models to explain the relationship between transcriptional regulation and the high decay rates of these RNAs. The first possibility is that the transient TF RNAs contain cis-acting elements, such as RNA-binding protein motifs. On the whole, we have not found evidence for universal RBP targeting of cell type-specific TF mRNAs, but it is possible that evolutionary pressures have selected sequence characteristics that destabilize these RNAs on a case-by-case basis and/or that destabilization depends on complex combinatorial action of multiple RNA-binding proteins that cannot be easily detected.

A second possibility is that transcriptional regulation of these RNAs is directly linked to their high degradation rates. Using data from another study (Pherson et al. 2017), we found that knockdown of PRC1 subunit Ph generally increases both the nascent and total RNA levels of PcG target RNAs, which argues against a direct effect of the PRC1 complex on RNA decay. However, this observation does not rule out the possibility that PcG targeting could be linked to RNA degradation of the same gene product through an unknown mechanism. This RNA destabilizing mechanism might remain active even after the gene is derepressed via experimental manipulation or in vivo during development. There is a growing body of evidence that factors that regulate the RNA life cycle can be deposited on the RNA during the process of transcription, a process that has been shown to depend on the promoter (Bregman et al. 2011; Dori-Bachash et al. 2012; Lavoie et al. 2012; Trcek et al. 2013; Catala and Abou Elela 2019). Likewise, it seems possible that factors that regulate transcriptional silencing could also recruit RNA destabilizing factors, which could be deposited on an RNA during transcription. An exciting candidate for a PRC linked destabilizing factor is the rixosome, a complex involved in rRNA processing that was recently shown to interact with PRC complexes (Zhou et al. 2022). Other candidates include the LSM2–8 complex, which was shown to repress the levels of mature PcG target RNAs in cooperation with Xrn2, presumably in the nucleus (Mattout et al. 2020). Future experiments involving careful manipulation of these complexes and measurements of the effects on both nuclear and cytoplasmic RNA decay rates will be required to clarify the roles of these pathways in PcG target RNA degradation.

Previous SectionNext Section

MATERIALS AND METHODS

Choice of labeling approach and protocol optimization

We considered different approaches for estimating RNA decay rates from our small tissue samples. We decided to use a pulse labeling approach with 4-thiouridine (4sU) and to extract RNA kinetic parameters from a single time point (Rabani et al. 2014; de Pretis et al. 2015). Classically, RNA decay has been measured with a pulse-chase approach. However, this method is difficult to perform with tissue samples. A pulse period and subsequent chase period in culture will expose the tissue to more time and potential damage from the ex vivo conditions. Furthermore, many chase time points are needed to capture the half-lives of RNAs at different ends of the half-life spectrum. In addition to the extra resources needed to process these numerous samples and their replicates, taking time points from tissue involves destructive sampling; that is they are collected from different sets of animals. Calculations that involve comparison between time points will compound error for measurement of RNA half-lives from noisier sources such as tissue. In the single time point method, both labeled RNA and total RNA are quantified from the same time point, and the resulting data are used to estimate RNA synthesis and decay rates using a steady-state model (Rabani et al. 2014; de Pretis et al. 2015).

Drosophila culture and RNA extraction

Drosophila melanogaster wild-type strain Oregon-R flies were raised on standard cornmeal-agar medium at 25°C. For each replicate 150 wandering third instar larvae were collected and their brains were dissected in brain culture media (BCM: 80% [v/v] Schneider's Drosophila medium [Thermo Fisher 21720024], 20% [v/v] heat-inactivated fetal bovine serum [FBS; Sigma F4135], 10 µg/mL human insulin [Sigma I9278]). The BCM used here is based on a previous recipe used for longterm culturing of larval brains (Hailstone et al. 2020) but modified to exclude whole larval extract, which contains RNases and nonbrain RNAs that could contaminate our sequencing libraries. Brains were collected in small batches of 15–20 in order to minimize total dissection and preparation time before labeling, which was kept to <20 min. Brains were transferred to 1.7-mL epitubes and the media was replaced with BCM containing 500 µM 4sU (Carbosynth NT06186) for 20 min. For the incubation condition testing experiments shown in Supplemental Figure S1, 50 brains were collected per replicate and the brains were incubated for 60 min. Our current version of the labeling protocol is able to recover adequate 4sU-labeled RNA from as few as 25 brains, corresponding to ∼5 µg of total RNA as input for the 4sU RNA purification. In addition to the methods described below, we now recommend using SPRI bead-based purification instead of chloroform extraction and using buffers that further increase 4sU-labeled RNA purification specificity (see the Supplemental Protocol provided at doi: 10.5281/zenodo.7699832).

Labeling was performed with tubes laid on their side with gentle shaking in the dark (80 rpm on an orbital shaker; New Brunswick Innova 2000). At the end of the labeling period, the brains were quickly washed with Schneiders media + 500 µM 4sU. All incubation and wash steps used 200 µL of the specified media. The wash with Schneider's medium was performed to remove traces of FBS, which we found to have high levels of a thiol-reactive molecule that could possibly interfere with our 4sU RNA purification approach. The samples were then flash-frozen in liquid nitrogen. Next, RNA was extracted with the hot acid phenol method. Brains were homogenized on ice in SEA buffer (50 mM NaOAc at pH 5.2, 10 mM EDTA, 0.3% Sarkosyl) for 20 sec using a motorized pestle. After homogenization, 1% SDS was added, and the extraction was continued as previously described (Thompson et al. 2020). Previous experiments using TRIzol for extraction led to sporadic failure of 4sU purification, which we believe was due to reducing agent carryover from TRIzol to the sample.

Purification of thiolated RNAs and sequencing library construction

RNA was biotinylated in a 360-µL reaction containing 10 mM HEPES (pH 7.5), 1 mM EDTA, and 24 µg of brain RNA. After adding other components, MTSEA biotin-XX (MTS-biotin; Biotium 90066) was added to a final concentration of 0.04 mg/mL as suggested (Duffy et al. 2015). The stock concentration of MTS-biotin was dissolved in anhydrous DMSO (ThermoFisher D12345) to a concentration of 0.2 mg/mL and stored at −80°C before use. The biotinylation reactions were purified with two rounds of extraction with an equal volume of chloroform and cleaned up with isopropanol precipitation. Pellets were resuspended in 50 µL of H2O and subjected to purification of on µMACs streptavidin columns as previously described, although reaction volumes were scaled down twofold (i.e., 50 µL of RNA added to 50 µL of beads) (Rädle et al. 2013). Pull-down RNA was isopropanol-precipitated and resuspended in 10 µL of H2O. All precipitations were done with 150 mM NaOAc (pH 5.2), 50% isopropanol, and 50 µg/mL GlycoBlue (ThermFisher AM9516) and then washed in cold 75% ethanol.

Ribosomal RNA depletion was then performed using biotinylated RNAs that are antisense to the 5S, 5.8S, 18S, and 28S ribosomal RNAs, as previously described (Stewart et al. 2010) with adaptations. We used in vitro transcribed antisense rRNA against the 5S, 5.8S, 18S, and 28S rRNAs (see below). For total RNA samples, 500 ng of RNA was mixed with asrRNA in a 50-µL hybridization reaction. All recovered pull-down RNA was mixed with asrRNA in a 25-µL hybridization reaction. For each 100 ng of input RNA, a mix containing 0.3 pmol of each asrRNA was added per each 100 ng of input RNA, which is approximately twofold molar excess over the rRNA in the sample. The input hybridization reaction was purified with 250 µL of Dynabeads MyOne Streptavidin C1 (Thermo Fisher, 65001) in 100 µL of binding volume. The pull-down hybridization reaction was purified with 25 µL of beads in 25-µL volume. Beads were washed in Dynabeads solution A prior to use. The binding reactions were incubated for 10 min at 22°C in the thermomixer at 600 rpm. The supernatant (rRNA-depleted) was then cleaned up and treated with DNase using the RNAqueous-Micro total RNA isolation kit (Thermo Fisher AM1931).

The RNA was eluted in 10 µL from the kit, and 6 µL was input into the library construction with the Lexogen Sense total RNA-seq library preparation kit for Illumina (Lexogen 009). Libraries were PCR-amplified using primers with unique 6mer i7 indices with 19 cycles (input libraries) or 18 cycles (pull-down libraries). Libraries were pooled and sequenced on an Illumina HiSeq 2500 and sequenced as two 125-bp paired-end reads with two sets of nine cycles of additional index read sequencing. Samples for the incubation condition testing experiments were processed similarly, except that ribosomal RNA subtraction was not performed and 500 ng of total RNA was used to construct poly(A)-primed libraries using the QuantSeq 3′ mRNA-seq library preparation kit FWD for Illumina (Lexogen 015) and amplified for 13–15 PCR cycles. Quant-seq libraries were loaded onto an Illumina NextSeq 500/550 high-output v2 cartridge and subjected to single-end sequencing on a NextSeq 500 instrument for 85 cycles with six cycles of index sequencing.

In vitro transcription and spike-ins

Drosophila rRNAs corresponding to the 5S, 5.8S, 18S, and 28S were cloned in the antisense orientation behind the T7 promoter in the ERCC plasmid backbone and linearized with BamHI. In vitro transcription was performed with the T7 Megascript kit (Thermo Fisher AM1334) with 40% of the UTP substituted with Biotin-16-UTP (Tebubio N-5005). Reactions were treated with Turbo DNase (Thermo Fisher AM2238), cleaned up with the RNA Clean & Concentrator-25 kit (Zymo Research R1017), and eluted with water.

We added two types of synthetic RNA spike-ins to use for quality control. Nonthiolated spike-in RNAs were from the Lexogen SIRV E2 mix (Lexogen 025.03), and thiolated spike-in RNAs (mix m2) were derived from a subset of the ERCC plasmids donated by NIST including ERCC-00057, ERCC-00073, ERCC-00077, ERCC-00104, ERCC-00112, ERCC-00142, ERCC-00150, ERCC-00162, ERCC-00165, and ERCC-00168. The thiolated spike-ins were made using T7 RNA polymerase (NEB M0251) on BamHI-digested templates. The transcription reactions contained 0.5 mM each NTP and additional 4sUTP (TriLink N-1025) added in proportion to the number of Us in each transcript: >350 Us, 0.125 mM; 190–350 Us, 0.25 mM; 100–190 Us, 0.5 mM; and <100 Us, 1 mM. The thiolated spike-ins were quantified with the Qubit RNA BR assay (Thermo Fisher Q10210) and mixed at equal final concentrations. The nonthiolated spike-ins were added to both the input and pull-down samples before ribosomal RNA removal. The thiolated spike-ins were added to the pull-down samples before biotinylation and pull-down, whereas they were added to the input samples before ribosomal RNA removal. For the 20-min labeling experiments, 2.25 ng of thiolated spike-ins was added to each input sample and 0.25 ng was added to each pull-down sample. An amount of 4.8 ng of Lexogen E2 mix was added to each input sample, and 6.1 ng was added to each pull-down sample. For the 60-min incubation test samples, 0.15 ng of thiolated spike-ins and 0.31 ng of Lexogen E2 mix were added to each total RNA sample.

RNA-seq data processing

Adapter removal and/or quality trimming were performed following recommendations from Lexogen. For SENSE RNA-seq libraries, the first 9 nt of the R1 read, the first 6 nt of the R2 read, and adapter sequences were removed with Cutadapt (Martin 2011) using the command “-m 20 -u 9 -U 6.” For the QuantSeq libraries, BBDuk (https://sourceforge.net/projects/bbmap/) was used with the parameters “k = 13 ktrim = r mink = 5 qtrim = r trimq = 10 minlength = 20.” Raw sequence data from BG3 cells that were generated using the Ion total RNA-seq kit v2 (Thermo Fisher 4475936) (Pherson et al. 2017) were downloaded from the SRA and used without further processing because the data had already been trimmed. Drosophila melanogaster genomic and transcript sequences were downloaded from Ensembl (release 99, assembly BDGP6.28). Kb-python (Melsted et al. 2021) was used to build a Kallisto (Bray et al. 2016) index with the coding, noncoding, and spike-in transcripts, as well intronic sequences with 30 nt of flanking sequence on either side. Kallisto was used to estimate counts mapping to each intron and each transcript using the options “–rf-stranded” for the Sense RNA-seq libraries and the options “-l 220 -s 130 –single –fr-stranded –single-overhang” and “-l 100 -s 100 –single –fr-stranded –single-overhang” for the QuantSeq and Ion total RNA-seq libraries, respectively. Transcript per million (TPM) values were then recalculated for both intronic and exonic regions for each gene after removing reads that were assigned to ribosomal RNA or spike-in RNAs. INSPEcT (de Pretis et al. 2015) was then used to estimate RNA synthesis, processing and decay rates from the 4sU-labeled and total RNA data. Snakemake (Koster and Rahmann 2012) and Conda (Grüning et al. 2018) were used to run the RNA-seq pipeline and manage software dependencies. Downstream analysis and figure construction were completed using open-source scientific computing software, including Numpy, Scipy, Pandas, and Matplotlib, run in Jupyter notebooks (Hunter 2007; McKinney 2010; van der Walt et al. 2011; Kluyver et al. 2016; Virtanen et al. 2020). Further details of the RNA-seq pipeline and the generated figures are available from Zenodo, doi: 10.5281/zenodo.7699832.

RNA-seq data analysis

Before downstream analyses, the data sets were filtered to exclude genes that were not adequately expressed. For each experiment, we chose a cutoff of 10 counts in at least 2/3 of the libraries for at least one condition. Because small RNAs are not recovered quantitatively by general RNA-seq protocols, in our analysis we focus on mRNAs and noncoding RNAs >200 nt (this includes both lncRNAs and other noncoding RNAs such as the signal recognition particle 7SL RNA). For analysis of decay and synthesis rate trends in specific gene groups versus other genes, the Mann–Whitney U-test was used. To calculate the error in the synthesis and decay measurements reported in Supplemental Table S1, we used the variance of the synthesis and total RNA levels reported by INSPEcT (de Pretis et al. 2015). We converted these variance measures to coefficient of variation (CV) values and also calculated CV values for degradation rates based on the formula used by INSPEcT for decay rate calculation at a steady-state [decay rate = synthesis rate/(total RNA − pre-mRNA)]. We observed that some genes had very long calculated half-lifes (>1000 min), but that the variation between half-lives from these genes had large interreplicate variation due to estimate uncertainty. We chose 1000 min as an upper limit based on the expected purification specificity of 4sU+ RNA relative to total RNA in the 4sU+ sample. Based on an empirically determined estimate from spike-in RNAs, we determined that some unlabeled RNAs will be present in the 4sU+ RNA sample. We chose a cutoff of a 1:100 ratio of 4sU+ to mature RNA in the total RNA library. This ratio corresponds to a calculated half-life of 1386 min, which we rounded to 1000. For scatterplots where experiments are compared, a pseudocount corresponding to the minimum value across experiments (i.e., min TPM) was added to all values before log transformation to allow plotting.

Comparison of RNA stability data with other data sets, gene lists, and gene features

Genes encoding transcription factors were obtained from the Flybase Gene Groups page. RNA-binding proteins were taken as genes mapping to GO term (RNA binding; GO:0003723) and mRNA-binding proteins were taken as genes mapping to GO term (mRNA binding; GO:0003729). Genes enriched in specific cell types in the L3 brain were taken from a single-cell RNA-seq data set and were defined as genes statistically enriched in a given cell type (P < 0.05) and with a fold-change enrichment of at least twofold (Corrales et al. 2022). PcG domain genes were taken from a study that identified them based on their high enrichment levels in both Ph Chip-seq and H3K27me3 ChIP-chip (Pherson et al. 2017). Neurite-localized RNAs were extracted from a meta-analysis and correspond to genes found to have significant neurite enrichment (P < 0.1) in at least three studies (Kügelgen and Chekulaeva 2020). DIOPT 9.0 (Hu et al. 2011) was used to extract Drosophila homologs of gene sets from other species. The DIOPT query was run using https://www.flyrnai.org/cgi-bin/DRSC_orthologs.pl and filtered to remove homologs found in fewer than eight sources (DIOPT score > 7). The search for enriched RBP motifs was done with Transite using the TSMA mode against their entire database of RBP motifs (Krismer et al. 2020). GO analysis was performed with clusterProfiler (Wu et al. 2021) and the returned categories were simplified with the semantic similarity score cutoff set to 0.5. For modeling decay rates as a function of gene features, decay rates were first log10 transformed and winsorized at the top and bottom 1%. The rates were fit to a linear model that considered the effect of each feature alone and in combination with each other (see Supplemental Fig. S6C). The model was run on either the entire data set or subsampled 10,000 times with bootstrapping. Bootstrapping was performed by subsampling each nonoverlapping gene group with a unique combination of features using a sample size equivalent to the smallest group size (n = 39). Regression was performed with the OLS function in the Statsmodels package (Seabold and Perktold 2010).

Previous SectionNext Section

DATA DEPOSITION

RNA sequencing data, including processed files summarizing RNA kinetic rates and expression levels, have been deposited in the Gene Expression Omnibus under accession GSE219202.

Previous SectionNext Section

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.

Previous SectionNext Section

ACKNOWLEDGMENTS

We thank the US National Institute of Standards and Technology (NIST) for the generous donation of the DNA plasmids used for making thiolated spike-in RNAs. We also thank Lexogen for providing the SIRV spike-ins used in this study and for performing part of the sequencing as part of the Lexogen Research Award. Sequencing of some samples was performed by the NGS Facility at the Vienna Biocenter Core Facilities (VBCF), member of Vienna Biocenter (VBC), Austria. Sequencing of other samples was performed in the Oxford Department of Zoology with technical assistance from Amanda Williams. This work was funded by a Wellcome Investigator Award 209412/Z/17/Z and a Leverhulme Trust Research Project Grant RPG-2017-393 to I.D. M.K.T. was also funded by the European Union's Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement 750928.

Author contributions: M.K.T. and I.D. designed experiments and wrote the manuscript. M.K.T. performed experiments and data analysis. A.C. contributed to data analysis via mathematical modeling.

Previous SectionNext Section

Footnotes

  • Received December 8, 2022.
  • Accepted March 14, 2023.

This article, published in RNA, is available under a Creative Commons License (Attribution 4.0 International), as described at http://creativecommons.org/licenses/by/4.0/.

Previous Section
 

REFERENCES

  1. ↵
    Beckervordersandforth R, Zhang C-L, Lie DC. 2015. Transcription-factor-dependent control of adult hippocampal neurogenesis. Cold Spring Harb Perspect Biol 7: a018879. doi:10.1101/cshperspect.a018879
    Abstract/FREE Full Text
  2. ↵
    Bray NL, Pimentel H, Melsted P, Pachter L. 2016. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol 34: 525–527. doi:10.1038/nbt.3519
    CrossRefMedlineGoogle Scholar
  3. ↵
    Bregman A, Avraham-Kelbert M, Barkai O, Duek L, Guterman A, Choder M. 2011. Promoter elements regulate cytoplasmic mRNA decay. Cell 147: 1473–1483. doi:10.1016/j.cell.2011.12.005
    CrossRefMedlineWeb of ScienceGoogle Scholar
  4. ↵
    Brody T. 1999. The interactive fly: gene networks, development and the internet. Trends Genet 15: 333–334. doi:10.1016/S0168-9525(99)01775-8
    CrossRefMedlineWeb of ScienceGoogle Scholar
  5. ↵
    Burow DA, Umeh-Garcia MC, True MB, Bakhaj CD, Ardell DH, Cleary MD. 2015. Dynamic regulation of mRNA decay during neural development. Neural Dev 10: 11. doi:10.1186/s13064-015-0038-6
    CrossRefGoogle Scholar
  6. ↵
    Catala M, Abou Elela S. 2019. Promoter-dependent nuclear RNA degradation ensures cell cycle-specific gene expression. Commun Biol 2: 211. doi:10.1038/s42003-019-0441-3
    CrossRefGoogle Scholar
  7. ↵
    Celniker SE, Dillon LAL, Gerstein MB, Gunsalus KC, Henikoff S, Karpen GH, Kellis M, Lai EC, Lieb JD, MacAlpine DM, et al. 2009. Unlocking the secrets of the genome. Nature 459: 927–930. doi:10.1038/459927a
    CrossRefMedlineWeb of ScienceGoogle Scholar
  8. ↵
    Chipman AD. 2020. The evolution of the gene regulatory networks patterning the Drosophila blastoderm. Curr Top Dev Biol 139: 297–324. doi:10.1016/bs.ctdb.2020.02.004
    CrossRefGoogle Scholar
  9. ↵
    Corrales M, Cocanougher BT, Kohn AB, Wittenbach JD, Long XS, Lemire A, Cardona A, Singer RH, Moroz LL, Zlatic M. 2022. A single-cell transcriptomic atlas of complete insect nervous systems across multiple life stages. Neural Dev 17: 8. doi:10.1186/s13064-022-00164-6
    CrossRefGoogle Scholar
  10. ↵
    de Pretis S, Kress T, Morelli MJ, Melloni GEM, Riva L, Amati B, Pelizzola M. 2015. INSPEcT: a computational tool to infer mRNA synthesis, processing and degradation dynamics from RNA- and 4sU-seq time course experiments. Bioinformatics 31: 2829–2835. doi:10.1093/bioinformatics/btv288
    CrossRefMedlineGoogle Scholar
  11. ↵
    Dillon N, Cocanougher B, Sood C, Yuan X, Kohn AB, Moroz LL, Siegrist SE, Zlatic M, Doe CQ. 2022. Single cell RNA-seq analysis reveals temporally-regulated and quiescence-regulated gene expression in Drosophila larval neuroblasts. Neural Dev 17: 7. doi:10.1186/s13064-022-00163-7
    CrossRefGoogle Scholar
  12. ↵
    Doe CQ. 2017. Temporal patterning in the Drosophila CNS. Annu Rev Cell Dev Biol 33: 219–240. doi:10.1146/annurev-cellbio-111315-125210
    CrossRefMedlineGoogle Scholar
  13. ↵
    Dori-Bachash M, Shalem O, Manor YS, Pilpel Y, Tirosh I. 2012. Widespread promoter-mediated coordination of transcription and mRNA degradation. Genome Biol 13: R114. doi:10.1186/gb-2012-13-12-r114
    CrossRefMedlineGoogle Scholar
  14. ↵
    Duffy EE, Rutenberg-Schoenberg M, Stark CD, Kitchen RR, Gerstein MB, Simon MD. 2015. Tracking distinct RNA populations using efficient and reversible covalent chemistry. Mol Cell 59: 858–866. doi:10.1016/j.molcel.2015.07.023
    CrossRefMedlineGoogle Scholar
  15. ↵
    Edgar BA, Weir MP, Schubiger G, Kornberg T. 1986. Repression and turnover pattern fushi tarazu RNA in the early Drosophila embryo. Cell 47: 747–754. doi:10.1016/0092-8674(86)90517-9
    CrossRefMedlineWeb of ScienceGoogle Scholar
  16. ↵
    Edgar BA, Odell GM, Schubiger G. 1989. A genetic switch, based on negative regulation, sharpens stripes in Drosophila embryos. Dev Genet 10: 124–142. doi:10.1002/dvg.1020100303
    CrossRefMedlineGoogle Scholar
  17. ↵
    Furlan M, Galeota E, del Gaudio N, Dassi E, Caselle M, de Pretis S, Pelizzola M. 2020. Genome-wide dynamics of RNA synthesis, processing, and degradation without RNA metabolic labeling. Genome Res 30: 1492–1507. doi:10.1101/gr.260984.120
    Abstract/FREE Full Text
  18. ↵
    Grosskortenhaus R, Pearson BJ, Marusich A, Doe CQ. 2005. Regulation of temporal identity transitions in Drosophila neuroblasts. Dev Cell 8: 193–202. doi:10.1016/j.devcel.2004.11.019
    CrossRefMedlineWeb of ScienceGoogle Scholar
  19. ↵
    Grüning B, Dale R, Sjödin A, Chapman BA, Rowe J, Tomkins-Tinch CH, Valieris R, Köster J. 2018. Bioconda: sustainable and comprehensive software distribution for the life sciences. Nat Methods 15: 475–476. doi:10.1038/s41592-018-0046-7
    CrossRefMedlineGoogle Scholar
  20. ↵
    Hafen E, Kuroiwa A, Gehring WJ. 1984. Spatial distribution of transcripts from the segmentation gene fushi tarazu during Drosophila embryonic development. Cell 37: 833–841. doi:10.1016/0092-8674(84)90418-5
    CrossRefMedlineWeb of ScienceGoogle Scholar
  21. ↵
    Hailstone M, Waithe D, Samuels TJ, Yang L, Costello I, Arava Y, Robertson E, Parton RM, Davis I. 2020. CytoCensus, mapping cell identity and division in tissues and organs using machine learning. Elife 9:e51085. doi:10.7554/eLife.51085
    CrossRefGoogle Scholar
  22. ↵
    Hu Y, Flockhart I, Vinayagam A, Bergwitz C, Berger B, Perrimon N, Mohr SE. 2011. An integrative approach to ortholog prediction for disease-focused and other functional studies. BMC Bioinformatics 12: 357. doi:10.1186/1471-2105-12-357
    CrossRefMedlineGoogle Scholar
  23. ↵
    Hunter JD. 2007. Matplotlib: a 2D graphics environment. Comput Sci Eng 9: 90–95. doi:10.1109/MCSE.2007.55
    CrossRefMedlineGoogle Scholar
  24. ↵
    Kassis JA, Kennison JA, Tamkun JW. 2017. Polycomb and trithorax group genes in Drosophila. Genetics 206: 1699–1725. doi:10.1534/genetics.115.185116
    Abstract/FREE Full Text
  25. ↵
    Kellerman KA, Mattson DM, Duncan I. 1990. Mutations affecting the stability of the fushi tarazu protein of Drosophila. Genes Dev 4: 1936–1950. doi:10.1101/gad.4.11.1936
    Abstract/FREE Full Text
  26. ↵
    Kluyver T, Ragan-Kelley B, Pérez F, Granger B, Bussonnier M, Frederic J, Kelley K, Hamrick J, Grout J, Corlay S, et al. 2016. Jupyter Notebooks—a publishing format for reproducible computational workflows. In Positioning and power in academic publishing: players, agents and agendas (ed. Loizides F and Schmidt B), pp. 87–90. IOS Press.
    Google Scholar
  27. ↵
    Koster J, Rahmann S. 2012. Snakemake—a scalable bioinformatics workflow engine. Bioinformatics 28: 2520–2522. doi:10.1093/bioinformatics/bts480
    CrossRefMedlineWeb of ScienceGoogle Scholar
  28. ↵
    Krismer K, Bird MA, Varmeh S, Handly ED, Gattinger A, Bernwinkler T, Anderson DA, Heinzel A, Joughin BA, Kong YW, et al. 2020. Transite: a computational motif-based analysis platform that identifies RNA-binding proteins modulating changes in gene expression. Cell Rep 32: 108064. doi:10.1016/j.celrep.2020.108064
    CrossRefGoogle Scholar
  29. ↵
    Kügelgen N, Chekulaeva M. 2020. Conservation of a core neurite transcriptome across neuronal types and species. Wiley Interdiscip Rev RNA 11: e1590. doi:10.1002/wrna.1590
    CrossRefGoogle Scholar
  30. ↵
    Landskron L, Steinmann V, Bonnay F, Burkard TR, Steinmann J, Reichardt I, Harzer H, Laurenson AS, Reichert H, Knoblich JA. 2018. The asymmetrically segregating lncRNA cherub is required for transforming stem cells into malignant cells. Elife 7: e31347. doi:10.7554/eLife.31347
    CrossRefGoogle Scholar
  31. ↵
    Lavoie M, Ge D, Abou Elela S. 2012. Regulation of conditional gene expression by coupled transcription repression and RNA degradation. Nucleic Acids Res 40: 871–883. doi:10.1093/nar/gkr759
    CrossRefMedlineGoogle Scholar
  32. ↵
    Lim B, Fukaya T, Heist T, Levine M. 2018. Temporal dynamics of pair-rule stripes in living Drosophila embryos. Proc Natl Acad Sci 115: 8376–8381. doi:10.1073/pnas.1810430115
    Abstract/FREE Full Text
  33. ↵
    Martin M. 2011. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J 17: 10. doi:10.14806/ej.17.1.200
    CrossRefMedlineGoogle Scholar
  34. ↵
    Mattout A, Gaidatzis D, Padeken J, Schmid CD, Aeschimann F, Kalck V, Gasser SM. 2020. LSM2-8 and XRN-2 contribute to the silencing of H3K27me3-marked genes through targeted RNA decay. Nat Cell Biol 22: 579–590. http://dx.doi.org/10.1038/s41556-020-0504-1.
    Google Scholar
  35. ↵
    McKinney W. 2010. Data structures for statistical computing in Python. In Proceedings of the 9th Python in Science Conference. Vol. 445, pp. 56–61. doi:10.25080/Majora-92bf1922-00a
    CrossRefGoogle Scholar
  36. ↵
    Meller VH, Wu KH, Roman G, Kuroda MI, Davis RL. 1997. roX1 RNA paints the X chromosome of male Drosophila and is regulated by the dosage compensation system. Cell 88: 445–457. doi:10.1016/S0092-8674(00)81885-1
    CrossRefMedlineWeb of ScienceGoogle Scholar
  37. ↵
    Melsted P, Booeshaghi AS, Liu L, Gao F, Lu L, Min KHJ, da Veiga Beltrame E, Hjörleifsson KE, Gehring J, Pachter L. 2021. Modular, efficient and constant-memory single-cell RNA-seq preprocessing. Nat Biotechnol 39: 813–818. doi:10.1038/s41587-021-00870-2
    CrossRefGoogle Scholar
  38. ↵
    Miller C, Schwalb B, Maier K, Schulz D, Dümcke S, Zacher B, Mayer A, Sydow J, Marcinowski L, Dölken L, et al. 2011. Dynamic transcriptome analysis measures rates of mRNA synthesis and decay in yeast. Mol Syst Biol 7: 458. doi:10.1038/msb.2010.112
    Abstract/FREE Full Text
  39. ↵
    Moon B-S, Bai J, Cai M, Liu C, Shi J, Lu W. 2018. Kruppel-like factor 4-dependent Staufen1-mediated mRNA decay regulates cortical neurogenesis. Nat Commun 9: 401. doi:10.1038/s41467-017-02720-9
    CrossRefMedlineGoogle Scholar
  40. ↵
    Mukherjee N, Calviello L, Hirsekorn A, de Pretis S, Pelizzola M, Ohler U. 2017. Integrative classification of human coding and noncoding genes through RNA metabolism profiles. Nat Struct Mol Biol 24: 86–96. doi:10.1038/nsmb.3325
    CrossRefMedlineGoogle Scholar
  41. ↵
    Neymotin B, Athanasiadou R, Gresham D. 2014. Determination of in vivo RNA kinetics using RATE-seq. RNA 20: 1645–1652. doi:10.1261/rna.045104.114
    Abstract/FREE Full Text
  42. ↵
    Pherson M, Misulovin Z, Gause M, Mihindukulasuriya K, Swain A, Dorsett D. 2017. Polycomb repressive complex 1 modifies transcription of active genes. Sci Adv 3: e1700944. doi:10.1126/sciadv.1700944
    FREE Full Text
  43. ↵
    Rabani M, Raychowdhury R, Jovanovic M, Rooney M, Stumpo DJ, Pauli A, Hacohen N, Schier AF, Blackshear PJ, Friedman N, et al. 2014. High-resolution sequencing and modeling identifies distinct dynamic RNA regulatory strategies. Cell 159: 1698–1710. doi:10.1016/j.cell.2014.11.015
    CrossRefMedlineWeb of ScienceGoogle Scholar
  44. ↵
    Rädle B, Rutkowski AJ, Ruzsics Z, Friedel CC, Koszinowski UH, Dölken L. 2013. Metabolic labeling of newly transcribed RNA for high resolution gene expression profiling of RNA synthesis, processing and decay in cell culture. J Vis Exp e50195.
    Google Scholar
  45. ↵
    Roy S, Ernst J, Kharchenko PV, Kheradpour P, Negre N, Eaton ML, Landolin JM, Bristow CA, Ma L, Lin MF, et al. 2010. Identification of functional elements and regulatory circuits by Drosophila modENCODE. Science 330: 1787–1797. doi:10.1126/science.1198374
    Abstract/FREE Full Text
  46. ↵
    Sami JD, Spitale RC, Cleary MD. 2022. mRNAs encoding neurodevelopmental regulators have equal N6-methyladenosine stoichiometry in Drosophila neuroblasts and neurons. Neural Dev 17: 9. doi:10.1186/s13064-022-00166-4
    CrossRefGoogle Scholar
  47. ↵
    Samuels TJ, Järvelin AI, Ish-Horowicz D, Davis I. 2020. Imp/IGF2BP levels modulate individual neural stem cell growth and division through myc mRNA stability. Elife 9: e51529. doi:10.7554/eLife.51529
    CrossRefMedlineGoogle Scholar
  48. ↵
    Schwanhäusser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, Chen W, Selbach M. 2011. Global quantification of mammalian gene expression control. Nature 473: 337–342. doi:10.1038/nature10098
    CrossRefMedlineWeb of ScienceGoogle Scholar
  49. ↵
    Seabold S, Perktold J. 2010. statsmodels: econometric and statistical modeling with python. Proceedings of the 9th python in science conference. http://statsmodels.sourceforge.net/
    Google Scholar
  50. ↵
    Simon JA, Kingston RE. 2013. Occupying chromatin: polycomb mechanisms for getting to genomic targets, stopping transcriptional traffic, and staying put. Mol Cell 49: 808–824. doi:10.1016/j.molcel.2013.02.013
    CrossRefMedlineWeb of ScienceGoogle Scholar
  51. ↵
    Stewart FJ, Ottesen EA, DeLong EF. 2010. Development and quantitative analyses of a universal rRNA-subtraction protocol for microbial metatranscriptomics. ISME J 4: 896–907. doi:10.1038/ismej.2010.18
    CrossRefMedlineWeb of ScienceGoogle Scholar
  52. ↵
    Thompson MK, Kiourlappou M, Davis I. 2020. Ribo-Pop: simple, cost-effective, and widely applicable ribosomal RNA depletion. RNA 26: 1731–1742. doi:10.1261/rna.076562.120
    Abstract/FREE Full Text
  53. ↵
    Trcek T, Sato H, Singer RH, Maquat LE. 2013. Temporal and spatial characterization of nonsense-mediated mRNA decay. Genes Dev 27: 541–551. doi:10.1101/gad.209635.112
    Abstract/FREE Full Text
  54. ↵
    van der Walt S, Colbert SC, Varoquaux G. 2011. The NumPy array: a structure for efficient numerical computation. Comput Sci Eng 13: 22–30. doi:10.1109/MCSE.2011.37
    CrossRefGoogle Scholar
  55. ↵
    Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, Burovski E, Peterson P, Weckesser W, Bright J, et al. 2020. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods 17: 261–272. doi:10.1038/s41592-019-0686-2
    CrossRefMedlineGoogle Scholar
  56. ↵
    Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, et al. 2021. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation 2: 100141. doi:10.1016/j.xinn.2021.100141
    CrossRefMedlineGoogle Scholar
  57. ↵
    Zhou H, Stein CB, Shafiq TA, Shipkovenska G, Kalocsay M, Paulo JA, Zhang J, Luo Z, Gygi SP, Adelman K, et al. 2022. Rixosomal RNA degradation contributes to silencing of Polycomb target genes. Nature 604: 167–174. doi:10.1038/s41586-022-04598-0
    CrossRefGoogle Scholar

Articles citing this article

  • Imp and Syp in vivo temporal RNA interactomes uncover networks of temporal regulators of Drosophila brain development bioRxiv July 2, 2024 0: 2024.06.30.600407v1-2024.06.30.600407
| Table of Contents
OPEN ACCESS ARTICLE

This Article

  1. Published in Advance April 11, 2023, doi: 10.1261/rna.079552.122 RNA 29: 1020-1032 © 2023 Thompson et al.; Published by Cold Spring Harbor Laboratory Press for the RNA Society

Article Category

ORCID

  1. Profile for Thompson, M. K.http://orcid.org/0000-0002-4947-6048
  2. Profile for Ceccarelli, A.http://orcid.org/0000-0002-9598-8845
  3. Profile for Ish-Horowicz, D.http://orcid.org/0000-0001-5684-7129
  4. Profile for Davis, I.http://orcid.org/0000-0002-5385-3053

Share

Navigate This Article

Current Issue

  1. July 2026, 32 (7)
  1. Current Issue
  1. Alert me to new issues of RNA