RNAs competing for microRNAs mutually influence their fluctuations in a highly non-linear microRNA-dependent manner in single cells

Abstract

Background: Distinct RNA species may compete for binding to microRNAs (miRNAs). This competition creates an indirect interaction between miRNA targets, which behave as miRNA sponges and eventually influence each other's expression levels. Theoretical predictions suggest that not only the mean expression levels of targets but also the fluctuations around the means are coupled through miRNAs. This may result in striking effects on a broad range of cellular processes, such as cell differentiation and proliferation. Although several studies have reported the functional relevance of this mechanism of interaction, detailed experiments are lacking that study this phenomenon in controlled conditions by mimicking a physiological range.

Results: We used an experimental design based on two bidirectional plasmids and flow cytometry measurements of cotransfected mammalian cells. We validated a stochastic gene interaction model that describes how mRNAs can influence each other's fluctuations in a miRNA-dependent manner in single cells. We show that miRNA-target correlations eventually lead to either bimodal cell population distributions with high and low target expression states, or correlated fluctuations across targets when the pool of unbound targets and miRNAs are in near-equimolar concentration. We found that there is an optimal range of conditions for the onset of cross-regulation, which is compatible with 10-1000 copies of targets per cell.

Conclusions: Our results are summarized in a phase diagram for miRNA-mediated cross-regulation that links experimentally measured quantities and effective model parameters. This phase diagram can be applied to in vivo studies of RNAs that are in competition for miRNA binding.

Keywords: Bimodality; MicroRNA target synchronization; Post-transcriptional cross-regulation; Single cell; Stochastic modelling.

Fig. 1

Model and predictions. a Sketch of the minimal model of miRNA–target interactions. One miRNA s and two targets r ₁ and r ₂ are independently transcribed with rates k _s, $k_{r_1}$, and $k_{r_2}$, respectively. Each transcript can then degrade with rate g _s, $g_{r_1}$, or $g_{r_2}$, respectively. Each miRNA s can interact with targets r ₁ or r ₂ with effective binding rates g ₁ or g ₂. α measures the probability of miRNA recycling. If not bound to a miRNA, targets r ₁ and r ₂ can be translated into proteins p ₁ and p ₂, respectively, which could then degrade with rates $g_{p_1}$ and $g_{p_2}$. b–d Predictions from the stochastic model of interactions sketched in (a) as a function of p ₀ (which is the constitutive value of p ₁ when g ₁ tends to 0) in terms of b the mean amount of p ₁ free molecules, c the p ₁ coefficient of variation ${\mathit{\text{CV}}}_{p_1}$, and d the Pearson correlation coefficient between p ₁ and p ₂. In (b–d), the red curve is the reference curve for a given set of parameters while the red line identifies the threshold. Blue and green curves show how the red curve would move when increasing the interaction strength with the second target g ₂ or the pool of miRNA via the miRNA transcription rate k _s, respectively. e Schematic representation of the two bidirectional plasmids coding for the four fluorophores. miRNA microRNA, UTR untranslated region

Fig. 2

Titration-induced threshold determines the optimal crosstalk. a, b mCherry mean fluorescence (a proxy for p ₁ in the model, Fig. 1 b) is plotted against eYFP (a proxy for the constitutive expression p ₀ in the model). Error bars are evaluated on the biological replicates. Continuous lines are model fits. The gray curves in (a) and (b) are the model prediction with the parameters fitted from the data and miRNA/target effective interaction strength g ₁→∞. The black arrow points to the model-predicted threshold. A threshold (or non-linear behavior) emerges when increasing mCherry MRE (a) while it disappears when increasing mCerulean MRE (b). The onset of the threshold is very close to the origin of the plot, indicating a relatively small amount of free miRNA. The intensity of crosstalk (measured in terms of fold-repression F with respect to the unregulated fluorophores) depends on the particular combination of MRE on both exogenous targets (c–e). F is the ratio between the value of mCherry in the absence of miR-20a MREs and its value in the presence of MREs for each eYFP bin and for each N on mCerulean. Purple and cyan circles in legends represent the plasmids coding for the mCherry and mCerulean fluorophores. a.u. arbitrary units, eYFP enhanced yellow fluorescent protein, MRE miRNA regulatory element

Fig. 3

miRNA increase shifts the maximal crosstalk region. a mCherry mean fluorescence (a proxy for p ₁ in the model, Fig. 1 b) is plotted against eYFP (a proxy for the constitutive expression p ₀ in the model). Blue triangles and red circles are data from cotransfection with pre-miR20a and negative controls, respectively. Error bars are evaluated on the biological replicates. The gray curve is the model prediction with the parameters fitted from the data and miRNA/target effective interaction strength g ₁→∞. The black arrow points to the model-predicted threshold. According to the model, increasing the pool of available miRNAs (transfecting pre-miRNAs) shifts the threshold to higher constitutive expression values. b Different combinations of miR-20a MREs lead to different levels of fold-repression and crosstalk. Triangles and circles in the plot are data from transfections with pre-miR20a and negative controls, respectively. Purple and cyan circles in legends represent the plasmids coding for mCherry and mCerulean fluorophores, respectively. a.u. arbitrary units, eYFP enhanced yellow fluorescent protein, MRE miRNA regulatory element

Fig. 4

Retroactivity increases cell-to-cell variability. a, b mCherry total noise, quantified by its coefficient of variation (CV, a proxy for ${\mathit{\text{CV}}}_{p_1}$ in Fig. 1 c), is plotted against eYFP (a proxy for the constitutive expression p ₀ in the model). The black arrow identifies the model-predicted threshold shown in Fig. 2. Error bars are evaluated on the biological replicates. CV globally increases on increasing the number of mCherry MREs (a) while it decreases on increasing the number of mCerulean MREs (b). The competition between these two strengths results in lowering the noise even if the expected repression from the rough number of mCherry MREs is high. Histograms in the lower panels show mCherry data distributions for the shaded regions in (a) and (b). A strong miRNA target repression strength increases cell-to-cell variability with the eventual appearance of different phenotypes (bimodal distributions). Purple and cyan circles in legends represent the plasmids coding for mCherry and mCerulean fluorophores, respectively. a.u. arbitrary units, CV coefficient of variation, eYFP enhanced yellow fluorescent protein, MRE miRNA regulatory element

Fig. 5

Fold Pearson and p values. The Pearson ratio is measured for three different values of eYFP basal expression: below threshold (a), around threshold (b), and above threshold (c). p values are reported for each combination of miRNA MREs on the two plasmids. The regions inside the blue perimeters are statistically significant with p<0.01. As predicted by the model, the correlation is maximal around the threshold and could be even 12-fold higher than in the unregulated case. Blue-delimited areas are regions whose Pearson ratio (i.e., the ratio of the Pearson coefficients between mCherry and mCerulean possessing different MREs for the same measure in the absence of MREs) is statistically relevant with respect to the corresponding unregulated case. eYFP enhanced yellow fluorescent protein, miRNA microRNA, MRE miRNA regulatory element

Fig. 6

Interplay between transcriptional activity and miRNA–target interaction strength. The figure shows model predictions and experimental results obtained when investigating the effect on one target (say p ₁) of the interplay between the second target (say r ₂ and, thus, p ₂) and the miRNA. The interplay between r ₂ and miRNA is tuned both via the transcription rate $k_{r_2}$ of r ₂ and via the interaction strength g ₂ between r ₂ and the miRNA. p ₁ is plotted against p ₀ on a increasing the transcription rate $k_{r_2}$ of r ₂, b increasing the interaction strength g ₂ between miRNA and r ₂ when $k_{r_2}>k_s$ (excess of targets), and c increasing the interaction strength g ₂ between miRNA and r ₂ when $k_{r_2}<k_s$ (excess of miRNA). The model prediction for cases depicted in (a) and (b) are qualitatively very similar. d mCherry mean fluorescence (a proxy for p ₁ in the model) is plotted against eYFP (a proxy for the constitutive expression p ₀ in the model). The dashed black line corresponds to the unregulated case while the blue data points correspond to the reference case with four MREs on mCherry and one MRE on mCerulean. Either increasing the copy number of mCerulean (a proxy for $k_{r_2}$ in the model), black data points, or the number of MREs on its sequence (a proxy for g ₂ in the model), red data points, has the effect of decreasing the amount of miRNA available to target mCherry (which globally increases). e Fold-repression with respect to the unregulated case plotted against eYFP. Error bars are evaluated on the biological replicates. a.u. arbitrary units, eYFP enhanced yellow fluorescent protein, miRNA microRNA, MRE miRNA regulatory element

Fig. 7

Phase diagram for mCherry (the target product p ₁). The figure shows how the crosstalk between targets and bimodality on mCherry behave on varying the effective miRNA interaction strength and the mean numbers of target mRNA molecules. The effective miRNA interaction strength on target r ₁ (and, thus, p ₁) is measured theoretically through the ratio g ₂/g ₁ and experimentally with different combinations of miRNA binding sites on both synthetic constructs