Skip to main content

Direct measurement of pervasive weak repression by microRNAs and their role at the network level

Abstract

Background

A gene regulatory network (GRN) comprises many weak links that are often regulated by microRNAs. Since miRNAs rarely repress their target genes by more than 30%, doubts have been expressed about the biological relevance of such weak effects. These doubts raise the possibility of under-estimation as miRNA repression is usually estimated indirectly from equilibrium expression levels.

Results

To measure miRNA repression directly, we inhibited transcript synthesis in Drosophila larvae and collected time-course data on mRNA abundance, the decline of which reflects transcript degradation. The rate of target degradation in the absence of miR310s, a moderately expressed miRNA family, was found to decrease by 5 to 15%. A conventional analysis that does not remove transcript synthesis yields an estimate of 6.5%, within the range of the new estimates. These data permit further examinations of the repression mechanisms by miRNAs including seed matching types, APA (alternative polyadenylation) sites, effects of other highly-expressed miRNAs and the length of 3’UTR. Our direct measurements suggest the latter two factors have a measurable effect on decay rate.

Conclusion

The direct measurement confirms pervasive weak repression by miRNAs, supporting the conclusions based on indirect assays. The confirmation suggests that this weak repression may indeed be miRNAs’ main function. In this context, we discuss the recent proposal that weak repression is “cumulatively powerful” in stabilizing GRNs.

Background

In any gene regulatory network (GRN), there is a small percentage of gene-to-gene interactions whereby the regulatory gene (usually a transcription factor, TF) up- or down-regulates its target genes by several fold, often measured using ChIP-seq -type analyses [1, 2]. The phenotypic consequences and functional significance of such strong interactions have been extensively documented [3, 4]. In contrast, a different group of regulatory molecules, namely the microRNAs (miRNAs), repress hundreds of targets per miRNA weakly and broadly [5,6,7] in the GRNs of metazoans. Each miRNA binds mainly to the 3′ untranslated regions (3’UTR) of target transcripts that have a sequence matching the 7–8 bp seed region of the miRNA. The binding induces transcript degradation and/or translational inhibition [5, 8,9,10].

The broad and weak repression by miRNAs has been a central conundrum in the control of GRNs. A view that has been increasingly supported is the “few targets” hypothesis [11,12,13,14]. In this view, only the most strongly repressed targets are biologically meaningful. In some cases, a single gene (such as the let-7 miRNA) is believed to repress a single target (lin-41) to yield a dramatic phenotypic consequence [12]. In contrast, Hunter et al. (2013), analyzing the same let-7 miRNA and the same phenotype in the same species, reach the opposite conclusion [15]. They suggest that more than 20 target genes are involved in controlling the vulva phenotype. At the center of the debate is the two aspects of the same issue: First, is there a single regulator-target interaction that is strong enough to exert a phenotypic effect? Second, if not, then what might be the functions of all those weak interactions?

Liufu et al. (2017) have recently shown that each regulatory molecule (miR310s in this case) governs the same phenotype through multiple targets incoherently [16]. For example, miR310s enhances egg hatchability by repressing the Mad target gene but reduces hatchability via three other targets, E2f2, EcR, and Mef2. They conclude that the many weak repressions function to stabilize the phenotype. If the task is to correct any deviations from the norm, weak effects should suffice. Chen et al. (2017; in review but posted on BioRxiv) introduces the May-Wigner theory used in studying food web stability to the analysis of GRNs [17]. They conclude that the weak miRNA regulatory effects are cumulative. Repression each of 50 targets by 1% has a greater effect on system stability than repression of one target gene by 50% in their simulations.

At the heart of the debate, theoretical arguments, and empirical tests is the assertion that miRNAs indeed weakly repress many targets. However, this conclusion is based on indirect measurements. When an miRNA is deleted, it is generally observed that the expression of predicted target genes increases [5]. However, the increase is usually very modest. It is often less than 20% for even highly-expressed miRNAs [18, 19], and even smaller for the less expressed molecules [20]. Such weak effects have led to two different views. In one view, the pervasive weak expression is biologically real as miRNAs exert their influence through the entire RNA:RNA network. Weak effects function cumulatively at the systems level [21,22,23,24].

In contrast, as one of the conventional explanations, the weak repression could be the result of biased measurements [25,26,27,28,29,30]. It is sometimes suggested that miRNAs’ effects should be measured at the level of proteins, not mRNAs [31,32,33]. Several studies address this issue directly [34,35,36,37] and report that mRNA and protein measurements are reasonably well correlated. Given this correlation, it seems unlikely that weak repressions of, say, 10% at the mRNA level would be translated into strong repressions of, say, 50% at the protein level.

In this study, we further explore the possibility that miRNA effects are under-estimated. Previous estimates are affected by both mRNA synthesis and degradation. Therefore, the rate of degradation could potentially be severely under-estimated if the two processes are inter-dependent as several studies have concluded [38,39,40]. Here, we measure miRNA target degradation directly by turning off transcription using Actinomycin D (ActD), a widely-used chemical for transcription inhibition [41]. We took care to apply the drug only for a short duration and monitored larval viability to reduce the likelihood of side effects.

The miR310 cluster (miR310s) of D. melanogaster is chosen for this study because its function, especially in the context of evolution in Drosophila, has been analyzed in some detail [42,43,44,45,46]. The mid-level expression level of this cluster [47, 48] also makes it a suitable candidate for addressing the estimation of miRNA repression levels.

Results

I. Theory

Let x i (t) denote the mRNA concentration of a miRNA target gene at time t. When the system is at an equilibrium, \( \frac{{\mathrm{dx}}_{\mathrm{i}}}{\mathrm{dt}}=0 \) (also dxj/dt = 0 for all other genes j). Near equilibrium, we approximate transcript change by a linear system. For gene i, the system is described by the following equation (Eq.)

$$ {\mathrm{dx}}_{\mathrm{i}}/\mathrm{dt}={\mathrm{B}}_{\mathrm{i}}\hbox{--} {\mathrm{D}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{i}} $$
(1)

where Bi (= bi + Si) is the synthesis rate and Di (= di + mi) is the degradation rate.

Here, bi is the presumably constant basal transcription rate and

$$ \mathrm{Si}={\sum}_{\mathrm{j}=1,\mathrm{j}\ne \mathrm{i}}^{\mathrm{N}} aijxj $$
(2)

is the aggregate effect of other genes on gene i with a ij being the regulation strength of gene i by gene j. Further, di is the degradation rate of transcript i in the absence of a matching miRNA and mi is the added degradation due to miRNA action. At equilibrium, Xi = Bi/Di where the capital letter indicates the equilibrium value. The equations can also be applied to the genetic background where the regulatory miRNA is deleted. In that background, Xi*, Bi* and Di* are substituted into Eqs. (1–2) and Xi* = Bi*/di (note that Di* = di).

To quantify the repression effect of miRNAs, one would compare the two equilibrium values, Xi (= Bi/Di) vs. Xi* (= Bi*/di). It is usually assumed that Bi/Bi* ~ 1. Hence, if Xi*/ Xi = 1.1, the interpretation is that Di/di = (di + mi)/di = 1.1 and mi is said to increase the degradation by about 10%. In reality, the ratio of Bi/Bi* is nearly impossible to quantify because Bi and Bi* both contain the summation term of Eq. (2), where xj’s and xj*‘s are very likely different. An incorrect assumption of Bi/Bi* ~ 1 could thus lead to potentially serious under-estimation of miRNAs repression activity. For example, if Xi*/ Xi = 1.1 and Bi/Bi* = 2, then (di + mi)/di = 2.2. In that case, the miRNA actually increases degradation by 120% instead of 10%.

Given that the confounding factor of Bi/Bi* cannot be estimated, we apply an Actinomycin D (ActD) solution to inhibit transcription, making Bi = Bi* = 0. Eq. (1) is reduced to

$$ {\mathrm{dx}}_{\mathrm{t}}/\mathrm{dt}=\hbox{-} \left(\mathrm{d}+\mathrm{m}\right)\ {\mathrm{x}}_{\mathrm{t}} $$
(3)
$$ {{\mathrm{dx}}_{\mathrm{t}}}^{\ast }/\mathrm{dt}=\hbox{-} \mathrm{d}\ {{\mathrm{x}}_{\mathrm{t}}}^{\ast } $$
(4)

in the wild-type (WT) and knockout (KO) background, respectively. As we only deal with gene i but at different time points, the subscript i is replaced by t. xt is measured by RNA-Seq at different time points after ActD is administered as described below.

II. Direct vs. indirect measurement of transcript decay

Transcription shutdown in vivo in Drosophila larvae

ActD is the standard reagent to repress both rRNA and mRNA synthesis [49, 50]. We chose the L3 stage during larval development for transcription repression because larvae appear to be able to withstand ActD treatment for 8 h without dying. When L3 larvae are immersed in an 80 μg/ml of ActD solution for 8 h, all larvae remain viable and do not appear to show obvious abnormalities. To ensure that a sufficient concentration of ActD is administered without causing undue damage to the cells, a preliminary experiment was done using a gradient of ActD concentration. As ActD affects rRNA synthesis at a much lower concentration than it affects mRNA synthesis, we expect the proper ActD concentration sufficient to disrupt mRNA production should lead to difficulties in pupation and eclosion. We chose 80 μg/ml for which the eclosion rate is still appreciable, more than 10% of the normal rate. In either the WT or miR310s-KO background, xt is measured by RNA-Seq at 4 time points: 0, 2, 4, and 8 h after ActD is administered at the pre-determined concentration. A high between-replicate correlation was observed (Additional file 1: Figure S1).

On the basis of Eqs (3) and (4), we obtain

$$ \ln\ \left({\mathrm{x}}_{\mathrm{t}}\right)=\ln\ \left(\mathrm{B}/\mathrm{D}\right)\hbox{--} \left(\mathrm{d}+\mathrm{m}\right)\ \mathrm{t} $$
(5)
$$ \ln\ \left({{\mathrm{x}}_{\mathrm{t}}}^{\ast}\right)=\ln\ \left({\mathrm{B}}^{\ast }/\mathrm{d}\right)\hbox{--} \left(\mathrm{d}\right)\ \mathrm{t}. $$
(6)

When ln (xt) or ln (xt*) is plotted against t (see Methods), the regression slope is d + m in the WT or d in the knock-out background. Plots of three representative transcripts including two miR310s non-target and one target gene are given in Fig. 1a-c. Note that the intercept in each panel is not informative because the intercept, unlike the slope, depends on the calibration of xt’s (see legends for detail).

Fig. 1
figure1

Degradation plots of (a-b) two untargeted transcripts, (c) one miR310s targeted transcript and (d) all transcripts. The expression level at 0 h is adjusted to 1. D is estimated by fitting the exponential model with log transformed normalized expression values. The slope corresponds to the decay rate

The distribution of the decay rate D (= d + m) in the WT background is given in Fig. 1d. In computing the distribution, we restrict the analysis to 4500 most highly-expressed transcripts to control for the uncertainty in estimating D for weakly expressed genes. The distribution has a long tail on the right indicating that a fraction of transcripts do decay with extreme rapidity. Because the corresponding distribution of D* (=d) in the knockout background exhibits the same trend (Fig. 1d), the analysis suggests that the rapid decay is not a consequence of miRNA action. Furthermore, there exists a small fraction of transcripts that decay at an exceptionally low rate. This can be seen if we plot the distribution of half-life of all transcripts (t1/2 = ln2/D) as shown in Additional file 2: Figure S2. The tail to the right indicates exceptionally slow decay. The degree of kurtosis seen in the combined Fig. 1d and Additional file 2: Figure S2 could make the comparison of two distributions (in the WT vs. KO background) more difficult to interpret because the outliers might exert strong influence on the mean and median.

Transcripts with high and low decay rate

We identify 200 most unstable and stable transcripts in the WT line and perform a Gene Ontology (GO) term enrichment analysis. The genes coding for unstable transcripts are disproportionally involved in development of the larval cuticle (Additional file 3: Table S1). This is consistent with the observations that numerous proteins degrade rapidly during larval development after chitin synthesis is inhibited [51]. By contrast, the stable transcripts tend to be implicated in translation (Additional file 4: Table S2), consistent with previous results [52,53,54,55].

We hypothesize that the fast decay of unstable transcripts is induced by AU-rich elements (ARE), GU-rich elements (GRE) and U-rich elements (URE), which are targeted by the RNA binding protein (RBP) [56,57,58,59,60]. Consistently, decay rates of mRNAs encoding these elements are statistically marginally higher than the remaining mRNAs in the WT line (Wilcoxon rank sum test, P = 0.09, Additional file 5: Figure S3). This pattern is much more pronounced in the KO line (P = 8.8 × 10− 15). The difference between two genetic backgrounds could be due to the fact that miR310s target three RBPs including Rb97D, CG4119, and CG4896. When miR310s are knocked out, RBPs get upregulated (Additional file 6: Table S3) and transcripts encoding elements like ARE are degraded even faster.

The decay rate of targets in the WT vs. KO background

To quantify the contribution of miRNAs to transcript degradation (i.e., mi’s), we identify 292 miR310s targets (see Methods). We measure the repression effect of miR310s by the equilibrium value - Xi (= Bi/Di) vs Xi* (= Bi*/di) (Fig. 2a). Indeed, the distribution of Xi* is shifted to the right suggesting higher expression of the target transcripts when the miR310s cluster is removed. The median value shows a 6.5% (=16.3/15.3–1) increase, confirming earlier suggestions of weak repression by miR310s.

Fig. 2
figure2

Distribution of the equilibrium expression and degradation change of 3’UTR targets between miR310s WT and KO. Targets’ degradation rates were normalized by the mode of the background decay rate. The median of equilibrium expression value is 15.3 for WT and 16.3 for KO. The degradation rate median value is 1.16 for WT and 1.11 for KO. (a) The cumulative distribution of equilibrium expression change. (b) The cumulative distribution of degradation change. (c) The contour distribution of miR310s targets’ degradation between WT and KO. Here, we performed a 2D kernel density estimation with an axis-aligned bivariate normal kernel, evaluated on a square grid and displaying the results with contours. The line indicates the contour estimated by the density of points. The red dot marks the densest point (1.04, 0.91). Strongly modulated target genes tend to be in the “high” or “long” group, which account for > 80% of the outliers

We then examine cumulative distributions of Di (= mi + di, in the WT background) vs. di (in the KO background, Fig. 2b). The distribution of Di is shifted to the right, as expected. However, the shift is even smaller than in Fig. 2a. In Fig. 2b, the median values are 1.16 and 1.11, respectively. The increase contributed by mi is only 4.50% (=1.16/1.11–1). Therefore, if we consider the median, the direct measurement of miRNA repression is smaller than the indirect measurement of miRNA effect based on equilibrium expression.

The analysis of the median in Fig. 2 corroborates the general view of weak repression by miRNAs. Nevertheless, the long tail of the decay distribution (Fig. 1d) makes it difficult to assess the population effect robustly. While the median is a conventional choice, we also considered the mode of the distribution as an additional inference statistic. Fig. 2c shows that the distribution of di and mi + di values where the contour plot reveals the densest concentration of points centers on (1.04, 0.91). This analysis is feasible because the degradation measurements among the 292 target transcripts are clustered relatively tightly and within the same order of magnitude, and the distribution is unimodal. Thus, using the distribution mode as the statistic, the overall degradation increase contributed by miR310s may be as high as 14.3% (1.04/0.91–1). To make sure our conclusion does not depend on specific bioinformatics parameters, we re-defines miR310s targets using D. yakuba, a more diverged outgroup species. Our observations remain essentially unchanged (Additional file 7: Figure S4).

We next checked if mi310s sites in coding sequences (in addition to the already considered 3’-UTRs) affect our conclusions. Such sites have been found to be functional before [61]. Analogous to the procedure for 3’-UTR targets, we extracted targets conserved in CDSs across D. melanogaster and D. simulans and examine both equilibrium and degradation change. The median equilibrium value of WT and KO is 13.0 and 14.6, respectively (Fig. 3a). That is, the increase is 12.3% (14.6/13.0–1). Correspondingly, the median degradation decrease is 11.4% (1.17/1.05–1) (Fig. 3b). However, the densest point in the contour plot shows that the degradation change is roughly − 2.0% (0.99/1.01–1) (Fig. 3c). Thus, miR310s changes degradation of CDS targets is about − 2.0% ~ 11.4%, which appears weaker than transcripts with binding sites in 3’ UTRs (4.5%~ 14.3%).

Fig. 3
figure3

Distribution of the CDS target equilibrium expression and degradation change between miR310s WT and KO. The median degradation rate is 1.17 for WT and 1.05 for KO. (a) The cumulative distribution of equilibrium expression change. (b) The cumulative distribution of degradation change. (c) The contour distribution of miR310s targets’ degradation between WT and KO. The red dot on the contour plot marks the densest point (0.99, 1.01). Strongly modulated target genes tend to be in the “high” or “long” group, which account for > 80% of the outliers

In summary, the degradation effect contributed by miR310s as measured the direct assay on transcripts with 3’-UTR binding sites is between 5 and 15% whereas the effect estimated from equilibrium expression is 6.5%. In the previous section, it was suggested that the two estimates might differ by as much as 10-fold (10 to 120%) due to the confounding factor of new transcript synthesis. Such a possibility is ruled out by this study. Given the weaker effect of targets in CDS, this possibility is even smaller for these targets.

III. Degradation factors other than miR310s

Transcripts targeted by miR310s at their 3’ UTRs differ in the number of target sites, matching types (i.e., 8mer-1a, 7mer-1a and 7mer-m8), alternative polyadenylation (APA), whether they are also co-targeted by other miRNAs, and the length of UTRs. Since all of these factors may affect degradation [62,63,64,65] we performed a series of binary analyses to test their potential effects.

First of all, since 90% of the targeted transcripts harbor only one miRNA site, we group them by matching types. Equilibrium expression gradually increases in the KO line according to the following order: 8mer-1a (15.1%), 7mer-1a (13.4%), and 7mer-m8 (11.9%) (Fig. 4a). Accordingly, degradation decreases by 9.4, 5.8, and 4.8%, respectively. Although the difference between groups is not significant, the trend is consistent with previous reports where 8mer-1a sites exert stronger repression intensity than 7mer-1a and 7mer-m8 [66, 67]. (We note that Figs. 2 and 4 use different statistics due to the nature of data partitioning. Figure 2 shows the difference in the median values of the control vs. KO lines whereas Fig. 4 shows the median value of the control-KO differences. The different treatments lead to the same conclusion.)

Fig. 4
figure4

Distribution of miR310s contribution to the equilibrium and degradation of different types of miR310s 3’UTR targets. We divide the targets by (a) seed match types (7mer-1a, 7mer-m8, and 8mer-1a), (b) potential alternative polyadenylation location, (c) whether targeted transcripts are also targeted by other highly-expressed miRNAs or not, and (d) UTR length. The left panels indicate equilibrium change and the right panels show degradation change. P values are from the Wilcoxon rank sum test

To determine whether a target can be subject to APA, we looked for the presence of the canonical polyadenylation site AATAAA or the non-canonical site ATTAAA. We classify targets into two groups: those with the APA site prior to the miR310s targeting site (the “APA” group) and those without an APA site (the “noAPA” group). Since APA results in a deletion of the miRNAs targeting site [63], we expect that the former group will show a smaller increase in expression. We observe the expected pattern: a 10.2% in the APA group and a 15.9% elevation in the non-APA group. Accordingly, the degradation decrease is 4.7% in APA and 6.0% in the non-APA group (Fig. 4b), although the between-group difference is not statistically significant.

Third, the miR310s cluster is only moderately expressed during the third instar larvae and one transcript could be targeted by multiple miRNAs. Based on the expression profile, we identify the top 20 most highly-expressed miRNAs. We then classify miR310s targets into two groups: “high” (transcripts co-targeted by highly-expressed miRNAs) and “low” (transcripts not targeted by highly-expressed miRNAs). Although the equilibrium changes are not significant between the “high” and “low” groups (12.0% vs. 10.1%, P > 0.05, Fig. 4c), the degradation changes are (− 2.7% vs. 15.2%, P < 0.05). Such a strong contrast suggests that during degradation, highly-expressed miRNAs exert even stronger inhibition effects on targets once miR310s are knocked out. Several possible mechanisms may account for this compensatory effect. One of them is the RNA:RNA cross talks via the sponge mechanism [68, 69]. Alternatively, the absence of one miRNA may facilitate the binding by other miRNAs. This alternative explanation is consistent with the actions of highly expressed miRNAs [45, 62, 70, 71].

Finally, we classify the targets into two groups: the subset of targets with the longest 3’ UTR among all the alternative transcripts (“long” subgroup) and another with relatively short UTR (“short” subgroup). The equilibrium changes are again not significant between these “long” and “short” groups (14.3% vs. 12.7%, P > 0.05, Fig. 4d), but the degradation rate changes are (2.0% vs. 12.8%, P < 0.05). Such a pattern suggests that the stability of an alternative transcript with longer UTR may be regulated by multiple mechanisms, with miR310s regulation representing only one layer.

In summary, the equilibrium changes are always insignificant in our tests, but the degradation changes show a significant change in two cases. Moreover, except for the “high” and “long” groups, degradation rate change is consistent with the equilibrium value change. Considering more than 50% of transcripts overlap between these two groups, the expression level may be affected by multiple factors besides degradation.

Discussion

The pervasive weak repression of target transcripts by miRNAs has been the central conundrum about miRNA function in animals. In particular, the repressions by miRNAs are generally smaller than the level of natural variation in transcript abundance. This weak effect has led to a published view (e.g., Pinzon et al.) that the bulk of repressions by miRNAs may not be biologically relevant. Nevertheless, this view itself is also in dispute [72].

There are many hypotheses. The simplest one posits that most of these repressions are merely noise [11]. An alternative hypothesis proposes that each miRNA governs phenotypes via multiple targets in coordination [16]. Neither hypothesis is supported by experimental evidence.

An alternative solution is to question the essence of the conundrum – that target repression is generally weak. As presented in the Introduction, the claim that weak repression at the mRNA level does not reflect the strength of repression at the protein level is not compatible with available results. In fact, there are few convincing reports that support miRNAs’ roles in reducing protein output strongly [34,35,36,37].

In this study, we propose an additional possibility that has attracted relatively little attention in the literature [52]. This new hypothesis is based on the transcription-degradation coupling. For example, closely related species of yeast tend to show a strong correlation between the rate of transcription and degradation. A gene product that is more strongly degraded tends to be more highly transcribed in the same species [38]. Several of the possible molecular mechanisms of the coupling have been suggested [38,39,40, 73,74,75]. If some of these mechanisms are incorporated into the miRNA control circuitry, the indirect measurement of the equilibrium level of gene expression would under-estimate miRNA repression effect. In theory, the indirect and direct measurements could differ by orders of magnitude. This coupling effect is expressed in the S term of Eq. (2).

By turning off transcription in Drosophila larvae, we remove this confounding effect from estimation. Most importantly, the direct and indirect measurements of miRNA repression overlap. Although the direct measurement of transcript decay yields more variable values (both technically and biologically), it is clear that the accepted view of pervasive weak repression by miRNAs appears valid.

Given this result, we return to the original question of miRNAs’ diffuse action, i.e., weak and broad repression. Several contending hypotheses attempt to explain this phenomenon, as recently summarized by Liufu et al. (2017) and Zhao et al. (2017) [16, 72]. Some of these hypotheses suggest that weak repression is generally devoid of function [11]. In this view, only a few of the targets are truly functional where small expression differences can lead to significant phenotypic consequences [12, 15, 76]. While this view can explain what is measurable, it leaves the bulk of miRNA repression unaccounted for. Against this backdrop, the “canalization” view which posits miRNAs’ functions in stabilizing the transcriptome and, hence, canalizing the phenotypes [44, 45, 70, 77] deserves more attention. The merit of the canalization view is that it takes into account the broad actions of weak repression, since all weak repression events cumulatively contribute to the stability of the transcriptome [72]. Recently, this canalization view has been expanded from the motif structure consisting of a few nodes [71, 78, 79] to the entire RNA network [80]. This current study, by corroborating the extent of weak repression, should help to re-energize the debate.

Conclusion

Gene regulation is expected to be direct, specific and sufficiently powerful in order to exert non-trivial phenotypic effects. Under this expectation, microRNAs (miRNAs) in metazoans are an enigmatic class of regulatory molecules. A central conundrum about miRNA function is the weak target repression. Since the new mRNA synthesis is usually not accurately accounted for, this repression effect may have been under-estimated in previous studies. We measured repression effects directly by turning off new target transcript synthesis and found that the repression effect is indeed as weak as the conventional assays suggest. Our data therefore rekindle the debate on the diffusive actions of miRNAs.

Methods

Fly culture and the miR310s knock-out stock

Fly larvae were reared at 25 °C and were fed a normal diet of corn and soybean meal, agar, and molasses supplemented with yeast. miR310s were knocked out by P-element transposition excision and maintained as a stock in the lab [81].

ActD treatment

We applied ActD (Sigma) that was soluble in the Phosphate Buffer Saline (PBS) at a raw concentration of 0.5 mg/ml. To get a suitable ActD concentration, different doses of ActD were tried and 80 μg/ml ActD was applied to sufficiently inhibit synthesis. L3 larvae were soaked in 800 μl of 80 μg/ml ActD in a 35 mm cell culture plate, with the solution just covering larval bodies. The larvae could breath under this volume of water but could not climb out of the liquid. RNA-Seq samples from WT and KO lines were collected at 0, 2, 4, and 8 h with two biological replicates. High correlation of all mRNAs and targets were observed between the two biological replicates (Additional file 1: Figure S1).

Reverse transcription – Quantitative real-time polymerase chain reaction (qRT-PCR)

The isolated RNA was reverse-transcribed into cDNA using the PrimeScript II 1st Strand cDNA Synthesis Kit (TaKaRa). cDNA was amplified using the primer sets listed in the Additional file 8: Table S3. RP49 was used as an internal control. The SYBR Premix Ex Taq (TaKaRa) was used in accordance with the manufacturer’s instructions. Quantitative real-time reverse transcription PCR analysis was performed using an applied biosystems 7900HT Real Time System (ThermoFisher Scientific). Decay rates estimated by qPCR were highly correlated with those derived from RNA-seq (Additional file 8: Figure S5). Primers are listed in Additional file 9: Table S4.

Library preparation and RNA-Seq

Total RNAs were extracted using the TRIzol Reagent (Ambion). RNA quality was assessed using 1% agarose gel electrophoresis. Five microgram of total RNA was used and polyA positive RNAs were isolated using the Dynabeads mRNA DIRECT Kit (Invitrogen). RNA-seq libraries were prepared according to the standard Illumina RNA-seq Library preparation protocol. Libraries were barcoded, pooled, and sequenced using Illumina HiSeq V4 with 125-bp paired reads. The raw sequence data were deposited in the Genome Sequence Archive [82] of the BIG Data Center [83], which is maintained by the Beijing Institute of Genomics (BIG) of the Chinese Academy of Sciences. The website is http://bigd.big.ac.cn/gsa and the accession number is PRJCA000381.

RNA-seq data analysis

The gene annotation file was downloaded from the Ensembl database. RNA-seq reads were mapped to D. melanogaster Reference Genome (BDGP6.83) using tophat [84] and transcripts’ expression value was quantified as “Fragments Per Kilobase of exon per Million fragments mapped” (FPKM) using Cufflinks [85].

Identification of targets

It has been illustrated that the conservation of miRNA’s target site sequences at 3′ untranslated regions (3’UTR) indicates functionality [86] and miRNA’s target sites at coding sequences (CDS) are also conserved and functional between species [61]. Since miR310s are relatively young (predating the D. melanogaster and D. pseudoobscura), we identified targets based on the conservation of seed sequence between D. melanogaster and two relatively closely related species: D. simulans and D. yakuba. miR310s CDS targets are conserved between D.melanogaster and D.simulans. Transcripts targeted by miR310s were identified using the TargetScan algorithm [7, 47, 87, 88]. 3’UTR annotation was extracted from BDGP6.83. The whole genome alignments were dm6, droSim1, and droYak3 available at the UCSC Genome Browser. miRNA expression during the L3 larval stage was ranked according to previous studies [47, 48].

Decay rate calculation

Although the lifetime of rRNA is measured in days [89] and the expression samples are collected at relatively early time points, normalization is indispensable between different time points. After we used FPKM to normalize for gene length and sequencing depth and chose 10 most stable ribosomal protein mRNAs that were relatively and increasingly expressed at later time points to perform normalization across samples. Due to experimental variance, long half-lives computed for stable transcripts were often inaccurate. Therefore, mRNAs with half-lives > 15 h were filtered out before further analysis. Moreover, after normalization, a subset of transcripts’ expression became very low and could lead to misleading inference of D. Therefore, transcripts with FPKM lower than 5 in any one sample were also filtered out. In eq. (5) and (6), normalized FPKM were used for linear fitting. D was estimated as the regression slope. Since the range of D is small (0–0.4, Fig. 1d), the mode (the most frequent value, 0.12 in WT and 0.14 in KO), in addition to the median, can be used to estimate the background decay rate pattern. To enable direct comparisons, we divided each individual value by the mode (0.12 in WT, 0.14 in KO) in each genetic background separately. All calculations were conducted in R [90].

Gene ontology analysis

The top 200 most stable and unstable transcripts were extracted to examine their Gene Ontologies. GO analysis enrichment was performed by GO:TermFinder with the 0.05 false discovery rate cutoff [91, 92].

mRNAs with ARE, GRE and URE motifs analysis

ARE, GRE and URE motifs refer to the sequence, ATTTA, AWTAAA, GTTTG, TTTGTTT, WTTTW, WWTTTWW, WWWTTTWWW, WWWWTTTWWWW, and WWWWWTTTWWWWW (W is A/U). Genes were extracted from the ARE2 website (http://nibiru.tbi.univie.ac.at/AREsite2/genes) [93]. Genes with the above listed motifs in 5’UTR, Exon, Intron, CDS and 3’UTR simultaneously were considered as potentially unstable mRNAs. We found about 490 transcripts fitting the criteria.

Abbreviations

*:

The state for KO

3’UTR:

3′ untranslated regions

ActD:

Actinomycin D

a ij :

Regulation strength of gene j on gene i

ARE:

AU-rich elements

b:

Constant basal transcription rate

B:

Synthesis rate

d:

Degradation rate for KO

D:

Degradation rate for WT

D. melanogaster :

Drosophila melanogaster

D. pseudoobscura :

Drosophila pseudoobscura

D. simulans :

Drosophila simulans

D. yakuba :

Drosophila yakuba

Eqs:

Equation

FPKM:

Fragments Per Kilobase of exon per Million fragments mapped

GRE:

GU-rich elements

GRN:

Gene regulatory network

KO:

Knockout

L3:

The third instar larvae

miR310s:

MiR310 cluster

miRNAs:

MicroRNAs

mRNA:

Message RNA

PBS:

Phosphate Buffer Saline

RNA-seq:

mRNA high-throughput sequencing

rRNA:

Ribosomal RNA

S:

Other genes’ effect on gene i

TFs:

Transcription factors

UCSC:

University of California Santa Cruz

URE:

U-rich elements

WT:

Wild-type

X:

Equilibrium state

xt :

mRNA concentration of a transcript at time t

References

  1. 1.

    Kilpinen H, Waszak SM, Gschwind AR, Raghav SK, Witwicki RM, Orioli A, Migliavacca E, Wiederkehr M, Gutierrez-Arcelus M, Panousis NI, et al. Coordinated effects of sequence variation on DNA binding, chromatin structure, and transcription. Sci. 2013;342:744–7.

    Article  CAS  Google Scholar 

  2. 2.

    Spitz F, Furlong EEM. Transcription factors: from enhancer binding to developmental control. Nat Rev Genet. 2012;13:613–26.

    Article  PubMed  CAS  Google Scholar 

  3. 3.

    Barrera LA, Vedenko A, Kurland JV, Rogers JM, Gisselbrecht SS, Rossin EJ, Woodard J, Mariani L, Kock KH, Inukai S, et al. Survey of variation in human transcription factors reveals prevalent DNA binding changes. Science. 2016;351:1450–4.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  4. 4.

    Kasowski M, Grubert F, Heffelfinger C, Hariharan M, Asabere A, Waszak SM, Habegger L, Rozowsky J, Shi MY, Urban AE, et al. Variation in transcription factor binding among humans. Sci. 2010;328:232–5.

    Article  CAS  Google Scholar 

  5. 5.

    Lim LP, Lau NC, Garrett-Engele P, Grimson A, Schelter JM, Castle J, Bartel DP, Linsley PS, Johnson JM. Microarray analysis shows that some microRNAs downregulate large numbers of target mRNAs. Nature. 2005;433:769–73.

    Article  PubMed  CAS  Google Scholar 

  6. 6.

    Friedman RC, Farh KKH, Burge CB, Bartel DP. Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 2009;19:92–105.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. 7.

    Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005;120:15–20.

    Article  PubMed  CAS  Google Scholar 

  8. 8.

    He L, Hannon GJ. Micrornas: small RNAs with a big role in gene regulation. Nat Rev Genet. 2004;5:522–31.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  9. 9.

    Bushati N, Cohen SM. MicroRNA functions. Annu Rev Cell Dev Biol. 2007;23:175–205.

    Article  PubMed  CAS  Google Scholar 

  10. 10.

    Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136:215–33.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  11. 11.

    Pinzon N, Li B, Martinez L, Sergeeva A, Presumey J, Apparailly F, Seitz H. microRNA target prediction programs predict many false positives. Genome Res. 2017;27:234–45.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. 12.

    Ecsedi M, Rausch M, Grosshans H. The let-7 microRNA directs Vulval development through a single target. Dev Cell. 2015;32:335–44.

    Article  PubMed  CAS  Google Scholar 

  13. 13.

    Flynt AS, Lai EC. Biological principles of microRNA-mediated regulation: shared themes amid diversity. Nat Rev Genet. 2008;9:831–42.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  14. 14.

    Seitz H. Redefining MicroRNA targets. Curr Biol. 2009;19:870–3.

    Article  PubMed  CAS  Google Scholar 

  15. 15.

    Hunter SE, Finnegan EF, Zisoulis DG, Lovci MT, Melnik-Martinez KV, Yeo GW, Pasquinelli AE. Functional genomic analysis of the let-7 regulatory network in Caenorhabditis elegans. PLoS Genet. 2013;9:e1003353.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. 16.

    Liufu ZQ, Zhao YX, Guo L, Miao GX, Xiao J, Lyu Y, Chen YX, Shi SH, Tang T, Wu CI. Redundant and incoherent regulations of multiple phenotypes suggest microRNAs' role in stability control. Genome Res. 2017;27:1665–73.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. 17.

    Chen YX, Shen Y, Allesina S, Wu CI. From foodwebs to gene regulatory networks (GRNs) – weak repressions by microRNAs confer system stability. bioRxiv. 2017; https://doi.org/10.1101/176701

  18. 18.

    Zhang P, Huang BH, Xu XR, Sessa WC. Ten-eleven translocation (Tet) and thymine DNA glycosylase (TDG), components of the demethylation pathway, are direct targets of miRNA-29a. Biochem Biophys Res Commun. 2013;437:368–73.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  19. 19.

    Ruckerl D, Jenkins SJ, Laqtom NN, Gallagher IJ, Sutherland TE, Duncan S, Buck AH, Allen JE. Induction of IL-4R alpha-dependent microRNAs identifies PI3K/Akt signaling as essential for IL-4-driven murine macrophage proliferation in vivo. Blood. 2012;120:2307–16.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. 20.

    Hafner M, Landthaler M, Burger L, Khorshid M, Hausser J, Berninger P, Rothballer A, Ascano M, Jungkamp AC, Munschauer M, et al. Transcriptome-wide identification of RNA-binding protein and MicroRNA target sites by PAR-CLIP. Cell. 2010;141:129–41.

  21. 21.

    McCann KS, Hastings A, Huxel G. R. Weak trophic interactions and the balance of nature. Nature. 1998;395:794–8.

    Article  CAS  Google Scholar 

  22. 22.

    Emmerson M, Yearsley JM. Weak interactions, omnivory and emergent food-web properties. Proc R Soc B-Biol Sci. 2004;271:397–405.

    Article  Google Scholar 

  23. 23.

    Gross T, Rudolf L, Levin SA, Dieckmann U. Generalized models reveal stabilizing factors in food webs. Science. 2009;325:747–50.

    Article  PubMed  CAS  Google Scholar 

  24. 24.

    Allesina S, Tang S. Stability criteria for complex ecosystems. Nature. 2012;483:205–8.

    Article  PubMed  CAS  Google Scholar 

  25. 25.

    Karres JS, Hilgers V, Carrera I, Treisman J, Cohen SM. The conserved microRNA MiR-8 tunes atrophin levels to prevent neurodegeneration in drosophila. Cell. 2007;131:136–45.

    Article  PubMed  CAS  Google Scholar 

  26. 26.

    Miska EA, Alvarez-Saavedra E, Abbott AL, Lau NC, Hellman AB, McGonagle SM, Bartel DP, Ambros VR, Horvitz HR. Most Caenorhabditis elegans microRNAs are individually not essential for development or viability. PLoS Genet. 2007;3:2395–403.

    Article  CAS  Google Scholar 

  27. 27.

    Chen YW, Song SL, Weng RF, Verma P, Kugler JM, Buescher M, Rouam S, Cohen SM. Systematic study of Drosophila MicroRNA functions using a collection of targeted knockout mutations. Dev Cell. 2014;31:784–800.

    Article  PubMed  CAS  Google Scholar 

  28. 28.

    Alvarez-Saavedra E, Horvitz HR. Many families of C. Elegans MicroRNAs are not essential for development or viability. Curr Biol. 2010;20:367–73.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. 29.

    Ge W, Chen YW, Weng R, Lim SF, Buescher M, Zhang R, Cohen SM. Overlapping functions of microRNAs in control of apoptosis during Drosophila embryogenesis. Cell Death Differ. 2012;19:839–46.

    Article  PubMed  CAS  Google Scholar 

  30. 30.

    Filipowicz W, Bhattacharyya SN, Sonenberg N. Mechanisms of post-transcriptional regulation by microRNAs: are the answers in sight? Nat Rev Genet. 2008;9:102–14.

    Article  PubMed  CAS  Google Scholar 

  31. 31.

    Strovas TJ, Rosenberg AB, Kuypers BE, Muscat RA, Seelig G. MicroRNA-based single-gene circuits buffer protein synthesis rates against perturbations. ACS Synth Biol. 2014;3:324–31.

    Article  PubMed  CAS  Google Scholar 

  32. 32.

    Schmiedel JM, Klemm SL, Zheng YN, Sahay A, Bluthgen N, Marks DS, van Oudenaarden A. MicroRNA control of protein expression noise. Science. 2015;348:128–32.

    Article  PubMed  CAS  Google Scholar 

  33. 33.

    Bartoszewska S, Kochan K, Madanecki P, Piotrowski A, Ochocka R, Collawn JF, Bartoszewski R. Regulation of the unfolded protein response by microRNAs. Cell Mol Biol Lett. 2013;18:555–78.

    Article  PubMed  CAS  Google Scholar 

  34. 34.

    Baek D, Villen J, Shin C, Camargo FD, Gygi SP, Bartel DP. The impact of microRNAs on protein output. Nature. 2008;455:64–U38.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  35. 35.

    Eichhorn SW, Guo HL, McGeary SE, Rodriguez-Mias RA, Shin C, Baek D, Hsu SH, Ghoshal K, Villen J, Bartel DP. mRNA destabilization is the dominant effect of mammalian MicroRNAs by the time substantial repression ensues. Mol Cell. 2014;56:104–15.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. 36.

    Selbach M, Schwanhausser B, Thierfelder N, Fang Z, Khanin R, Rajewsky N. Widespread changes in protein synthesis induced by microRNAs. Nature. 2008;455:58–63.

    Article  PubMed  CAS  Google Scholar 

  37. 37.

    Jonas S, Izaurralde E. NON-CODING RNA towards a molecular understanding of microRNA-mediated gene silencing. Nat Rev Genet. 2015;16:421–33.

    Article  PubMed  CAS  Google Scholar 

  38. 38.

    Dori-Bachash M, Shema E, Tirosh I. Coupled evolution of transcription and mRNA degradation. PLoS Biol. 2011;9:e1001106.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. 39.

    Goler-Baron V, Selitrennik M, Barkai O, Haimovich G, Lotan R, Choder M. Transcription in the nucleus and mRNA decay in the cytoplasm are coupled processes. Genes Dev. 2008;22:2022–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. 40.

    Shalem O, Dahan O, Levo M, Martinez MR, Furman I, Segal E, Pilpel Y. Transient transcriptional responses to stress are generated by opposing effects of mRNA production and degradation. Mol Syst Biol. 2008;4:223.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  41. 41.

    Scherrer K, Latham H, JE D. Demonstration of an unstable RNA and of a precursor to ribosomal RNA in HeLa cells. Proc Natl Acad Sci U S A. 1963;49:240–8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  42. 42.

    Aravin AA, Lagos-Quintana M, Yalcin A, Zavolan M, Marks D, Snyder B, Gaasterland T, Meyer J, Tuschl T. The small RNA profile during Drosophila melanogaster development. Dev Cell. 2003;5:337–50.

    Article  PubMed  CAS  Google Scholar 

  43. 43.

    Lu J, Fu YG, Kumar S, Shen Y, Zeng K, Xu AL, Carthew R, Wu CI. Adaptive evolution of newly emerged micro-RNA genes in Drosophila. Mol Biol Evol. 2008;25:929–38.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. 44.

    Tang TA, Kumar S, Shen Y, Lu JA, Wu ML, Shi SH, Li WH, Wu CI. Adverse interactions between micro-RNAs and target genes from different species. Proc Natl Acad Sci U S A. 2010;107:12935–40.

    Article  PubMed  PubMed Central  Google Scholar 

  45. 45.

    Wu CI, Shen Y, Tang T. Evolution under canalization and the dual roles of microRNAs-A hypothesis. Genome Res. 2009;19:734–43.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. 46.

    Ninova M, Ronshaugen M, Griffiths-Jones S. Fast-evolving microRNAs are highly expressed in the early embryo of Drosophila virilis. RNA. 2014;20:360–72.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  47. 47.

    Ruby JG, Stark A, Johnston WK, Kellis M, Bartel DP, Lai EC. Evolution, biogenesis, expression, and target predictions of a substantially expanded set of Drosophila microRNAs. Genome Res. 2007;17:1850–64.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. 48.

    Lyu Y, Shen Y, Li H, Chen YX, Guo L, Zhao YX, Hungate E, Shi SH, Wu CI, Tang T. New MicroRNAs in Drosophila-birth, death and cycles of adaptive evolution. PLoS Genet. 2014;10:e1004096.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. 49.

    Reich EFR, Shatkin AJ, Tatum EL. Effect of actinomycin D on cellular nucleic acid synthesis and virus production. Science. 1961;134:556–7.

    Article  PubMed  CAS  Google Scholar 

  50. 50.

    Perry R, Kelley D. Inhibition of RNA synthesis by actinomycin D: characteristic dose-response of different RNA species. J Cell Physiol. 1970;76:127–39.

    Article  PubMed  CAS  Google Scholar 

  51. 51.

    Gangishetti U, Breitenbach S, Zander M, Saheb SK, Muller U, Schwarz H, Moussian B. Effects of benzoylphenylurea on chitin synthesis and orientation in the cuticle of the Drosophila larva. Eur J Cell Biol. 2009;88:167–80.

    Article  PubMed  CAS  Google Scholar 

  52. 52.

    Thomsen S, Anders S, Janga SC, Huber W, Alonso CR. Genome-wide analysis of mRNA decay patterns during early Drosophila development. Genome Biol. 2010;11:R93.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. 53.

    Sharova LV, Sharov AA, Nedorezov T, Piao Y, Shaik N, Ko MSH. Database for mRNA half-life of 19 977 genes obtained by DNA microarray analysis of pluripotent and differentiating mouse embryonic stem cells. DNA Res. 2009;16:45–58.

    Article  PubMed  CAS  Google Scholar 

  54. 54.

    Yang E, van Nimwegen E, Zavolan M, Rajewsky N, Schroeder M, Magnasco M, Darnell JE. Decay rates of human mRNAs: correlation with functional characteristics and sequence attributes. Genome Res. 2003;13:1863–72.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  55. 55.

    Narsai R, Howell KA, Millar AH, O'Toole N, Small I, Whelan J. Genome-wide analysis of mRNA decay rates and their determinants in Arabidopsis thaliana. Plant Cell. 2007;19:3418–36.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. 56.

    Vlasova IA, Tahoe NM, Fan D, Larsson O, Rattenbacher B, John JRS, Vasdewani J, Karypis G, Reilly CS, Bitterman PB, et al. Conserved GU-rich elements mediate mRNA decay by binding to CUG-binding protein 1. Mol Cell. 2008;29:263–70.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  57. 57.

    Scheiba RM, de Opakua AI, Diaz-Quintana A, Cruz-Gallardo I, Martinez-Cruz LA, Martinez-Chantar ML, Blanco FJ, Diaz-Moreno I. The C-terminal RNA binding motif of HuR is a multi-functional domain leading to HuR oligomerization and binding to U-rich RNA targets. RNA Biol. 2014;11:1250–61.

    Article  PubMed  Google Scholar 

  58. 58.

    Lu YC, Chang SH, Hafner M, Li X, Tuschl T, Elemento O, Hla T. ELAVL1 modulates transcriptome-wide miRNA binding in murine macrophages. Cell Rep. 2014;9:2330–43.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. 59.

    Halees AS, Hitti E, Al-Saif M, Mahmoud L, Vlasova-St Louis IA, Beisang DJ, Bohjanen PR, Khabar KSA. Global assessment of GU-rich regulatory content and function in the human transcriptome. RNA Biol. 2011;8:681–91.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  60. 60.

    Lee JE, Lee JY, Wilusz J, Tian B, Wilusz CJ. Systematic analysis of cis-elements in unstable mRNAs demonstrates that CUGBP1 is a key regulator of mRNA decay in muscle cells. PLoS One. 2010;5:e11201.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  61. 61.

    Liu GJ, Zhang R, Xu J, Wu CI, Lu XM. Functional conservation of both CDS- and 3 ‘-UTR-located MicroRNA binding sites between species. Mol Biol Evol. 2015;32:3276.

    CAS  Google Scholar 

  62. 62.

    Androsavich JR, Chau BN. Non-inhibited miRNAs shape the cellular response to anti-miR. Nucleic Acids Res. 2014;42:6945–55.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  63. 63.

    Boutet SC, Cheung TH, Quach NL, Liu L, Prescott SL, Edalati A, Iori K, Rando TA. Alternative polyadenylation mediates MicroRNA regulation of muscle stem cell function. Cell Stem Cell. 2012;10:327–36.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  64. 64.

    Grimson A, Farh KKH, Johnston WK, Garrett-Engele P, Lim LP, Bartel DP. MicroRNA targeting specificity in mammals: determinants beyond seed pairing. Mol Cell. 2007;27:91–105.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  65. 65.

    Boehm V, Haberman N, Ottens F, Ule J, Gehring NH. 3 ‘ UTR length and messenger ribonucleoprotein composition determine Endocleavage efficiencies at termination codons. Cell Rep. 2014;9:555–68.

    Article  PubMed  CAS  Google Scholar 

  66. 66.

    Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. elife. 2015;4:e05005.

    Article  PubMed Central  Google Scholar 

  67. 67.

    Rodriguez A, Vigorito E, Clare S, Warren MV, Couttet P, Soond DR, van Dongen S, Grocock RJ, Das PP, Miska EA, et al. Requirement of bic/microRNA-155 for normal immune function. Science. 2007;316:608–11.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  68. 68.

    Ebert MS, Neilson JR, Sharp PA. MicroRNA sponges: competitive inhibitors of small RNAs in mammalian cells. Nat Methods. 2007;4:721–6.

    Article  PubMed  CAS  Google Scholar 

  69. 69.

    Thomson DW, Dinger ME. Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet. 2016;17:272–83.

    Article  PubMed  CAS  Google Scholar 

  70. 70.

    Hornstein E, Shomron N. Canalization of development by microRNAs. Nat Genet. 2006;38:S20–4.

    Article  PubMed  CAS  Google Scholar 

  71. 71.

    Cui QH, Yu ZB, Purisima EO, Wang E. Principles of microRNA regulation of a human cellular signaling network. Mol Syst Biol. 2006;2:46.

    Article  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Zhao YX, Shen X, Tang T, Wu CI. Weak regulation of many targets is cumulatively powerful-an evolutionary perspective on microRNA functionality. Mol Biol Evol. 2017;34:3041–6.

    Article  PubMed  CAS  Google Scholar 

  73. 73.

    Shalem O, Groisman B, Choder M, Dahan O, Pilpel Y. Transcriptome kinetics is governed by a genome-wide coupling of mRNA production and degradation: a role for RNA pol II. PLoS Genet. 2011;7:e1002273.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  74. 74.

    Harel-Sharvit L, Eldad N, Haimovich G, Barkai O, Duek L, Choder M. RNA polymerase II subunits link transcription and mRNA decay to translation. Cell. 2010;143:552–63.

    Article  PubMed  CAS  Google Scholar 

  75. 75.

    Xu S, He Z, Zhang Z, Guo Z, Guo W, Lyu H, Li J, Yang M, Du Z, Huang Y, et al. The origin, diversification and adaptation of a major mangrove clade (Rhizophoreae) revealed by whole-genome sequencing. Natl Sci Rev. 2017; https://doi.org/10.1093/nsr/nwx065

  76. 76.

    Wen H, Wang HY, He X, Wu CI. On the low reproducibility of cancer studies. Natl Sci Rev. 2018; https://doi.org/10.1093/nsr/nwy021

  77. 77.

    Waddington C. Canalization of development and the inheritance of acquired characters. Nature. 1942;150:563–5.

    Article  Google Scholar 

  78. 78.

    Levy SF, Siegal ML. Network hubs buffer environmental variation in Saccharomyces cerevisiae. PLoS Biol. 2008;6:2588–604.

    Article  CAS  Google Scholar 

  79. 79.

    Rutherford SLLS. Hsp90 as a capacitor for morphological evolution. Nat. 1998;396:336–42.

    Article  CAS  Google Scholar 

  80. 80.

    Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell. 2011;146:353–8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  81. 81.

    Zhao T, Gu TT, Rice HC, KL MA, Roark KM, Lawson K, Gauthier SA, Reagan KL, Hewes RS. A Drosophila gain-of-function screen for candidate genes involved in steroid-dependent neuroendocrine cell remodeling. Genetics. 2008;178:883–901.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  82. 82.

    Wang Y, Song F, Zhu J, Zhang S, Yang Y, Fang X, Lei H, Zhang Z, Zhao W. GSA: genome sequence archive. Genomics, Proteomics & Bioinformatics. 2017;15:14–8.

    Article  Google Scholar 

  83. 83.

    Zhang Z, Zhao WM, Xiao JF, Song SH, Hao LL, Li RJ, Ma LN, Sheng X, Sang J, Wang YQ, et al. The BIG data center: from deposition to integration to translation. Nucleic Acids Res. 2017;45:D18–24.

    Article  CAS  Google Scholar 

  84. 84.

    Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  85. 85.

    Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2014;9:2513.

    Article  CAS  Google Scholar 

  86. 86.

    Xu J, Zhang R, Shen Y, Liu GJ, Lu XM, Wu CI. The evolution of evolvability in microRNA target sites in vertebrates. Genome Res. 2013;23:1810–6.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  87. 87.

    Ruby JG, Jan CH, Bartel DP. Intronic microRNA precursors that bypass Drosha processing. Nat. 2007;448:83–6.

    Article  CAS  Google Scholar 

  88. 88.

    Kheradpour P, Stark A, Roy S, Kellis M. Reliable prediction of regulator targets using 12 Drosophila genomes. Genome Res. 2007;17:1919–31.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  89. 89.

    Defoiche J, Zhang Y, Lagneaux L, Pettengell R, Hegedus A, Willems L, Macallan DC. Measurement of ribosomal RNA turnover in vivo by use of deuterium-labeled glucose. Clin Chem. 2009;55:1824–33.

    Article  PubMed  CAS  Google Scholar 

  90. 90.

    R [https://www.r-project.org/]. Accessed 02 Feb 2017.

  91. 91.

    Boyle EI, Weng SA, Gollub J, Jin H, Botstein D, Cherry JM, Sherlock G. GO::TermFinder - open source software for accessing gene ontology information and finding significantly enriched gene ontology terms associated with a list of genes. Bioinformatics. 2004;20:3710–5.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  92. 92.

    Blake JA, Christie KR, Dolan ME, Drabkin HJ, Hill DP, Ni L, Sitnikov D, Burgess S, Buza T, Gresham C, et al. Gene ontology consortium: going forward. Nucleic Acids Res. 2015;43:D1049–56.

    Article  CAS  Google Scholar 

  93. 93.

    Fallmann J, Sedlyarov V, Tanzer A, Kovarik P, Hofacker IL. AREsite2: an enhanced database for the comprehensive investigation of AU/GU/U-rich elements. Nucleic Acids Res. 2016;44:D90–5.

    Article  PubMed  CAS  Google Scholar 

Download references

Acknowledgments

We would like to thank Suhua Shi, Tian Tang, Jin Xu, Zhongqi Liufu and Guangan Lu for their insightful comments.

Funding

This study was supported by grants from the National Key Basic Research Program of China to C.I.W. (2014CB542006), the Strategic Priority Research Program of the Chinese Academy of Sciences to X.M.L and C.I.W. (XDB13040300) and the Natural Science Foundation of China to X.M.L. (91531305).

Availability of data and materials

The raw sequence data has been deposited in the Genome Sequence Archive of the BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences. The website is http://bigd.big.ac.cn/gsa and the accession number PRJCA000381.The open software used in this study is cited and materials can be requested from the authors.

Author information

Affiliations

Authors

Contributions

C.I.W and X.M.L designed the study. F.Q.M performed the experiment. C.I.W, Y.Z, F.Q.M, P.L and Q.J.C analyzed the data. C.I.W, Y.Z and F.Q.M wrote the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Yong E. Zhang or Chung-I Wu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S1. Correlation of all transcript and target expression (log2 transformed) between two biological replicates at 0 h, 2 h, 4 h, and 8 h. (a-d) results from the WT line, (e-h) results from the KO line. (TIF 368 kb)

Additional file 2:

Figure S2. Half-lives of whole transcripts in D. melanogaster third instar larvae. (TIF 157 kb)

Additional file 3:

Table S1. Relating mRNA stability to gene function using a Gene Ontology (GO) analysis in the WT line (the 200 least stable genes). (PDF 133 kb)

Additional file 4:

Table S2. Relating mRNA stability to gene function using a Gene Ontology (GO) analysis in the WT line (the 200 most stable genes). (PDF 163 kb)

Additional file 5:

Figure S3. Cumulative distribution of mRNA decay rates of with or without AU-rich, GU-rich and U-rich elements in 3’UTR, 5’UTR, Intron, Exon, and CDS simultaneously in (a) the WT and (b) KO line. AU-rich, GU-rich and U-rich element: ATTTA, AWTAAA, GTTTG, TTTGTTT, WTTTW, WWTTTWW, WWWTTTWWW, WWWWTTTWWWW, and WWWWWTTTWWWWW. W: A/U. P values are from Wilcoxon rank sum tests. (TIF 98 kb)

Additional file 6:

Table S3. The decay rates of three targets that belong to the RNA binding protein (RBP) family. (PDF 196 kb)

Additional file 7:

Figure S4. (a-c) Distribution of the equilibrium expression and degradation change of 3’UTR targets conserved between D. melanogaster and D. yakuba. Targets’ degradation rates were normalized by the mode of the background decay rate. The median of equilibrium expression value is 15.3 for WT and 16.4 for KO. The median of degradation rate is 1.19 for WT and 1.14 for KO. (a) The cumulative distribution of equilibrium expression change. (b) The cumulative distribution of degradation change. (c) The contour distribution of miR310s targets’ degradation between WT and KO. The red dot on the contour plot marks the densest point (1.04, 0.92). (TIF 212 kb)

Additional file 8:

Figure S5. The relationship between the decay rate calculated by qRT-PCR and mRNA-seq. RP49 is used as the reference. (TIF 72 kb)

Additional file 9:

Table S4. Real-time polymerase chain reaction primers used for mRNA quantification. (PDF 198 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ma, F., Lin, P., Chen, Q. et al. Direct measurement of pervasive weak repression by microRNAs and their role at the network level. BMC Genomics 19, 362 (2018). https://doi.org/10.1186/s12864-018-4757-z

Download citation

Keywords

  • mRNA degradation
  • MiR310s
  • Targets
  • Weak repression
  • Network