Dosage and Temporal Thresholds in microRNA Proteomics*

MicroRNAs (miRNAs) modulate protein and mRNA expression through translational repression and/or mRNA decay. In this study, we combined SILAC-based proteomics and RNAseq to identify primary targets based on measurements of protein and mRNA repression and analysis of transcript 3′UTR sequences. The primary target set was used to compare different prediction algorithms, revealing higher stringency of selection by Targetscan and PITA compared with miRanda, at the expense of higher false negatives. A key finding was that significant and unexpected variations occurred in the kinetics of repression as well as the sensitivity to exogeneous miRNA concentration. Bimodal thresholds were observed, which distinguished responses to low (10 nm) versus high (50–100 nm) miRNA, as well as the onset of repression at early (12–18 h) versus late (36–48 h) times. Similar behavior was seen at the transcript level with respect to kinetics of repression. The differential thresholds were most strongly correlated with ΔΔG, the net free energy of miRNA-target interactions, which mainly reflected inverse correlations with ΔGopen, the free energy of forming 3′UTR secondary structures, at or nearby the miRNA seed matching sites. Thus, our working model is that protein binding or other competitive mechanisms variably interfere with the accessibility of miRISC to the transcript binding site. In addition, biphasic responses were observed in a subset of proteins that were partially down-regulated at early times, and further down-regulated at later times. Taken together, our findings provide evidence for varying modes of miRNA target repression, which lead to different thresholds of target responses with respect to kinetics and concentration, and predict that certain transcripts will show graded responses in sensitivity and fold-change under cellular conditions that lead to varying steady state miRNA levels.

microRNAs (miRNAs) 1 are small noncoding RNAs that modulate cellular protein levels, impacting a wide range of processes ranging from cellular development and differentiation to cancer, metabolic disorder, and other human diseases (1)(2)(3)(4)(5). Mature miRNAs interact with Argonaute (AGO) family proteins to form miRNA-induced silencing complexes (miRISCs) (6,7), which in turn attenuate protein expression through translational repression and/or mRNA decay mechanisms. In animal cells, miRISC interacts with mRNA-bound poly(A) binding protein C (PABPC) through a complex between AGO and the GW182 trinucleotide-repeat-containing protein. This interaction recruits deadenylation, 5Ј-decapping, and mRNA degrading enzymes, and also blocks eukaryotic translation initiation factor G (eIFG) binding to confer translational repression (8 -11). Although the degree of protein or mRNA repression by miRNAs is often moderate (1.5-to 4-fold), fine-tuning at this level can be enough to shift cell and organismal phenotypes, or to set thresholds for gene and protein expression changes which control the robustness of cellular states (12)(13)(14)(15).
The specificity by which proteins are targeted by miRISC is determined by base-pairing between a "seed-sequence" at the 5Јend of a miRNA and a "seed-matching sequence" in a targeted mRNA, typically located at the 3ЈUTR. Seed sequences are only 6 -8 nt, therefore, a single miRNA can potentially target hundreds of different mRNAs. But many gene products containing seed-matching sequences are not true targets, leading to ambiguities in distinguishing primary targets from off-target (false positive) responses. Several algorithms have been developed to improve the accuracy of identifying miRNA primary targets (16 -21). These quantify scores based on additional features beyond seed sequence complementarity, determined from and validated by microarray evidence for mRNA repression. Such features include: 1) thermodynamic stability of miRNA-mRNA base pairing, 2) position and accessibility of seed matching sequences within the 3ЈUTR (21), 3) miRNA-mRNA base pairing in regions flanking the seed sequence (19), and 4) evolutionary conservation of seed-matching sequences, the latter is predictive for miRNAs that are conserved across mammalian and non-mammalian species (17,18,20). Although several studies have used SILAC-based LC-MS/MS methods to examine miRNA responses, it is not known how different prediction algorithms compare when evaluated by largescale proteomics.
Data sets used to train and validate scores used by prediction algorithms prioritize responses measured within 24 h after transfection (19). A few studies have examined mRNA responses with respect to time and miRNA concentration. For example, responses to let-7 miRNA showed repression of cell cycle genes within 16 -24 h, and postulated that early responding genes favor functional targets, whereas later responses are more likely to be indirect (22). In addition, mRNA repression was maximal at 20 nM siRNA added exogenously to HeLa cells (23,24), suggesting that repression of primary targets saturates at low miRNA concentrations. However, in these studies, only one or a few target proteins were typically analyzed, and on a global scale, miRNA responses might in fact vary widely with respect to time and concentration. This has been rarely addressed at the mRNA level and has been unexplored at the protein level.
Here, we present a systematic, proteome-wide analysis of the kinetics and dose responses of protein changes in response to a single miRNA. Our findings show unanticipated variations among primary miRNA targets with respect to their sensitivity to exogenous miRNA concentration and the time at which they reach the threshold for significant change. Targets respond to miRNA concentration in a bimodal manner, where the more sensitive targets show higher free energies of miRNA-mRNA hybridization and reduced 3ЈUTR structural stability. Likewise, targets can be grouped into early versus late responding groups, which are also correlated with 3ЈUTR stability. Similar behavior is seen with mRNA transcripts, where mRNA destabilization and decay appears mostly concomitant with protein repression. The results show that among primary targets, wide variations in kinetics and dose response occur that reflect differences between transcripts with respect to miRNA-mRNA and 3ЈUTR stability. This expands the understanding of how the regulation of primary miRNA targets may vary under different cellular states.
Cell Culture-WM239A metastatic melanoma cells were cultured in RPMI supplemented with dialyzed 10% FBS. Static SILAC experiments were carried out as described (25,26), by growing cells in RPMI/FBS containing heavy-or light-labeled ArgϩLys for at least five passages, before transfection with miRNA mimics or NTs. Pulsed SILAC (pSILAC) experiments were carried out as described (27), by maintaining cells in RPMI/FBS containing light ArgϩLys, which upon transfection was changed to either RPMI/FBS containing heavy or medium ArgϩLys.
Cell Transfection-miRNA mimics or nontargeting controls were mixed with Dharmafect 1 in serum-free RPMI and incubated at room temperature for 30 min. The miRNA-lipid mixture was added to 1 ϫ 10 6 cells on a 6 cm plate in RPMI/FBS. In static SILAC experiments, the RPMI/FBS was changed 24 h after transfection and cells were harvested 24 h later. In pSILAC experiments, RPMI/FBS was not changed during time courses to avoid inconsistent media exposure. Cells transfected with nontargeting miRNAs had growth rates and adherence comparable to untransfected cells, thus showing no evidence for cytotoxicity.
Cells were harvested by adding lysis buffer (0.15 ml) containing 4% SDS, 20 mM DTT, and 0.1 M Tris, pH 8.6, to each plate. Lysates were incubated at 95°C for 5 min, sonicated, and then clarified by centrifugation (14,000 ϫ g, 15 min). Proteins were quantified using the microplate BCA protein assay, and 125 g each of miRNA-treated cell lysates and NT-treated cell lysates were mixed and processed using the filter-aided sample preparation (FASP) protocol (28). Briefly, cell lysates were loaded onto Vivacon 500 spin columns and washed with 8 M urea, 0.1 M Tris, pH 8.6, (0.2 ml). Proteins were alkylated with iodoacetamide (50 mM), washed with urea/tris solution, and equilibrated with 0.1 M ammonium bicarbonate (0.2 ml) followed by overnight proteolysis with 1% (w/w) trypsin at 37°C. Tryptic peptides were recovered through the filter by centrifugation, dried by speedvac centrifugation, and stored at Ϫ80°C. Peptides were dissolved in 36 L 0.1% (v/v) formic acid and 10 L (28%) was analyzed by 2D-LC-MS/MS.
Data Analysis-Raw data were processed using MaxQuant (version 1.0.13.13) (29) and Mascot (version 2.2.0, Matrix Science) used to search the peak list against IPI human database version 3.52 (73,928 entries) (30) with 262 common contaminant entries (31). The search allowed trypsin specificity with two missed cleavages, and included Cys carbamidomethylation as fixed and Met oxidation as variable modifications. MaxQuant used 7 ppm maximum initial mass deviation for the precursor ion, and 0.5 Da MS/MS tolerance, searching six top MS/MS peaks per 100 Da. False discovery rates were 0.01 for both peptide and protein identifications, with 0.5 maximum posterior envelope probability (PEP) for peptides, six amino acid minimum peptide length, and two minimum total peptides, of which at least one peptide must be unique (i.e. non-razor) for each protein identification. Protein quantifications used both razor and unique peptides, allowed unmodified peptides and peptides with oxidized Met, and required at least two ratio measurements with ratio variability Ͻ 75% for each protein. Protein ratios were reanalyzed and accepted when more than two-thirds of peptide ratios were within 0.15 from the median peptide ratio. The representative ID, which was supported by the greatest number of peptides in a protein group was used for further analysis. Protein abundance ratios between experimental and control conditions were accepted only when ratios of at least two different peptides per protein were measured. Normalized protein ratios were used in all experiments, except for experiments that measured protein half-life. The identified proteins and protein ratios from each analysis are summarized in supplemental Table S1.
RNASeq Experiments-WM239A cells were transfected with 50 nM miR-22 or NT, and cells were harvested at 12, 18, 24, 36, and 48 h post-transfection. Total RNA was purified using RNAeasy Mini kit (Qiagen), of which 1 g RNA was processed using Illumina TruSeq RNAseq Sample Preparation kit v2 (Cat# RS-122-2001) to generate libraries. Library quality was verified using the Bioanalyzer DNA 1000 assay (Agilent Technologies), and sequencing was performed using the Illumina HiSeq 2000 system with single read (1 ϫ 50 bp) option. Reads were processed using TopHat2 (version 2.0.16) (35,36) and Cufflinks (version 2.2.0) (37) programs with the Tuxedo protocol (38) for alignment of short reads, transcript assembly, and analysis of differential expression. On average, there were ϳ13 million reads per sample, with Ͼ11,600 genes quantified at each time point with an acceptable quality by Cufflinks. Genes were quantified only when fragment per kilobase of transcript per million mapped fragments (FPKM) were greater than or equal to one under two conditions (e.g. miR-22 and NT).

Protein Abundance Changes in Response to miR-22-
WM239A cells derived from metastatic melanoma were labeled by SILAC with heavy (H) or light (L) isotopically-labeled ArgϩLys, and transfected with 50 nM of nontargeting (NT) miRNA or miR-22, previously identified as an miRNA whose expression is regulated by oncogenic mutant B-Raf signaling in melanoma (39). Proteins from whole cell lysate were trypsinized and the resulting peptides were analyzed by LC-MS/MS (Fig. 1A). Using a 2D-UPLC RP-RP LTQ-Orbitrap configuration with step gradient elutions in the first dimension followed by linear gradient elution in the second dimension, more than 4000 proteins were typically quantified after 56 h of data collection (supplemental Table S1).
Log 2 (protein ratio) normalized to nontargeting (NT) miRNA versus no miRNA (mock) showed good correlation (Fig. 1B, Pearson's r ϭ 0.73), indicating similarity between different control conditions. Therefore, only NT was used as the control in further experiments. Protein ratios were in good agreement when miR-22-and NT-transfected cells were reversely labeled ( Fig. 1C, r ϭ 0.86), or when protein ratios were measured in biological duplicate experiments (Fig. 1D, r ϭ 0.83). This supported reproducibility with no bias introduced by labeling order. Based on comparisons between duplicate experiments, protein ratios corresponding to log 2 (miR-22/NT) Ͼ 0.3 were accepted as significant with 1% false discovery rate (Fig. 1E). Of 3444 proteins quantified reproducibly, 350 were significantly up-regulated and 400 down-regulated in response to miR-22. Ninety percent of the significant protein changes showed ratios less than 2-fold, consistent with previous proteomics analyses of miRNA responses (27, 40 -44). Thus, protein ratio measurements were accurate and reproducible, and enabled quantification of 1.23-fold changes in proteins responsive to miR-22.
Classifications of mRNA seed matching sequences included: seed-matches to miRNA positions 2-8 with adenosine at target position 1 (8mer); seed-matches to miRNA positions 2-8 (7mer m8); seed-matches to miRNA positions 2-7 plus adenosine at target position 1 (7mer A1); and seed-matches to miRNA positions 2-7 (6mer) (19,45,46). Genomics and proteomics studies have shown that proteins corresponding to mRNAs with at least one 7mer or 8mer matching sequence within the 3ЈUTR are most strongly repressed by miRNAs (40,46). We observed the same behavior in response to miR-22 (Fig. 1F). In total, 156 proteins were both responsive to miR-22 and corresponded to mRNAs containing at least one 7mer or 8mer seed-matching sequence, and can be considered primary targets. They are referred to as "seed-matches" and are the focus of this study.
miRNA Target Prediction Algorithms-We examined the degree to which different methods predicted the targets of miR-22 that we identified by proteomics. Of the proteins observed in biological duplicates, 534 corresponded to mRNAs containing seed-matching sequences (8mer or any 7mer), of which 156 were significantly repressed by miR-22 ( Fig. 2A). The remaining 378 proteins were unaffected, providing an estimate of 71% false positives predicted based on seed-matching sequence alone. Of the 400 proteins that were significantly down-regulated, 244 (61%) lacked corresponding mRNA seed-matching sequences, suggesting that they are regulated indirectly. Cumulative distribution plots showed greater fold-changes in proteins corresponding to mRNAs with 8mer or 7mer seed-matches, compared with proteins with no seed-match (p Ͻ 10 Ϫ6 , one-sided Mann-Whitney U test, "M-W") ( Fig. 2B).
The primary targets were then used to analyze predictions by different algorithms, each applied using their default cutoffs and/or stringency thresholds. miRanda (http://www. microrna.org) combines evolutionary conservation of seedmatching sequences and calculated binding energies for miRNA-mRNA interactions with a flexible algorithm for sequence complementary and position within mRNA 3ЈUTR or coding sequences (16,17). The predicted targets include those with imperfect base pairing between seed:seed-matching sequences. Of the proteins that were quantified by proteomics, miRanda predicted 259 as miR-22 targets. Of these, 96 were down-regulated at the protein level, yielding 163 (63%) estimated false positives (Fig. 2C). Of the 156 primary targets, 69 were excluded by miRanda, yielding 44% false negatives. We then filtered predicted targets according to the presence of 8mer or 7mer seed matching sequences within corresponding transcripts. A subset of 194 miRanda-predicted targets had perfect 3Ј UTR seed matches, of which 87 were down-regulated, indicating 107 (55%) false positives. In contrast, 65 miRanda-predicted targets had 6mer sequences, imperfect seed base-pairing, and/or were located outside of the 3ЈUTR. Of these, only nine were down-regulated, indicating higher (86%) false-positives in this category. Cumulative distributions were comparable between all miRanda-predicted targets, the subset with perfect base-pairing, and targets predicted solely based on seed-matching sequences 99% of measurements were within log 2 (difference from mean) Ͻ 0.3 (red lines), from which protein changes of log 2 (miR-22/NT) Ͼ 0.3 were accepted as significant with 1% false discovery rate. F, Cumulative distribution showing proteins with transcripts with 8mer seed matching sequence (red), 7mer m8 (blue), 7mer A1 (cyan), 6mer only (green), or no site (black).
( Fig. 2D), suggesting that miRanda performs comparably with predictions based on seed-matching sequence alone.
Targetscan (http://targetscan.org) uses a model for site efficacy based on correlations between target repression and various sequence "contexts" (e.g. additional base pairing in the miRNA 3Ј region, local AU content, and position within the 3Ј UTR), as well as evolutionary conservation of 3ЈUTR seed matching sequences (19,20). Targetscan predicted 69 proteins as miR-22 targets (Fig. 2E), of which 36 were downregulated at the protein level, yielding 33 (48%) false positives. Of 156 primary targets, 120 were excluded by Targetscan, yielding 77% false negatives. Targetscan-predicted proteins showed stronger fold-repression compared with those predicted solely by seed-matching sequence (Fig.  2F), showing that the algorithm filters targets with low magnitude fold-change while accepting higher numbers of false negatives.
The Probability of Interaction by Target Accessibility (PITA) algorithm evaluates the net free energy of miRNA-target interactions (⌬⌬G), calculated as the difference between the free energy gain from miRNA-mRNA hybridization (⌬G duplex ) minus the free energy loss from 3ЈUTR secondary structure melting (⌬G open ) (21). PITA predicted 52 proteins as miR-22 targets (Fig. 2G), of which 28 were down-regulated at the protein level, yielding 24 (46%) false positives. Of the 156 primary targets, 128 were excluded by PITA, yielding 82% false negatives. Fold-changes in protein repression were higher among PITA-predicted proteins compared with those predicted solely from seed-matching sequence (Fig. 2H). Thus, PITA and Targetscan both filter targets with low magnitude fold-changes at the expense of high false negatives. In total, miRanda, PITA, and Targetscan identified 110 primary targets, of which only 13 were predicted by all algorithms. Approximately 50% of primary targets identified by Targetscan overlapped with targets identified by PITA and vice versa, whereas ϳ25% of targets identified by miRanda overlapped with either PITA or Targetscan (supplemental Fig. S1). Estimates of false positives and false negatives identified by each program are summarized in Table I. miRNA responses were then correlated with scores evaluated by prediction algorithms. Fig. 2I shows small but significant correlations (Pearson's r ϭ Ϫ0.32, p Ͻ 10 Ϫ6 ) between log 2 (protein ratio) and the Targetscan "probability of preferentially conserved targeting score (P CT )," which calculates a Bayesian estimate of the probability that a seed-match is conserved because of selective pressure during evolution (20). Correlations were also seen (r ϭ 0.34, p Ͻ 10 Ϫ6 ) with the Targetscan "context score," which estimates seed-match efficacy based on features of the 3ЈUTR (seed-match type contribution, base-pairing with miRNA 3Ј sequences, local AU content, and position within the 3ЈUTR), regardless of sequence conservation (19). In contrast, Fig. 2J shows low correlation between protein ratio and the PITA-calculated net free energy of miRNA-target interaction (⌬⌬G; r ϭ Ϫ0.11, p ϭ 0.04), which reflects the free energy gain from miRNA-mRNA hybridization (⌬G duplex ; r ϭ Ϫ0.01, p ϭ 0.88) minus the free energy of transcript 3ЈUTR folding (⌬G open ; r ϭ Ϫ0.13, p ϭ 0.02). This reflects a bias of Targetscan but not PITA for predicting targets with higher fold-changes.
Variations in Sensitivity to miRNA Concentration-We asked whether primary targets showed biases in their sensitivity to exogenous miRNA, by monitoring the dose-response of protein ratios with varying miRNA concentrations. Cells were transfected with 10, 25, 50, or 100 nM miR-22 or NT and harvested 48 h post-transfection. Levels of intracellular miR-22 increased linearly with the concentration of transfected miRNA, reaching up to 42,000 copies/pg at 50 nM miR-22 (supplemental Fig. S2), an amount comparable to physiological levels of other miRNAs (47,48). Plots comparing log 2 (protein ratios) showed strong correlations between experiments at different miR-22 concentrations (Pearson's r ϭ 0.77-0.88) (Fig. 3A). However, an unexpected discontinuity was observed, where the slope of the correlation plots equaled Ϸ1 when comparing 10 nM versus 25 nM or 50 nM versus 100 nM (Fig. 3A, left and right panels, respectively), but increased to Ϸ1.8 when comparing 25 nM versus 50 nM (Fig.  3A, middle panel). This indicates that a transition occurs above 25 nM and below 50 nM, which increases the magnitude of protein ratios. The transition between 25 and 50 nM was also apparent in plots of cumulative distribution versus log 2 (protein ratio), for the subset of proteins corresponding to mRNAs with 8mer or 7mer seed-matches (Fig. 2B). Thus, primary targets show discontinuous responses to miRNA concentration, in parallel with those seen proteome-wide.
We focused on 137 of the 156 primary seed-matched targets, excluding those with larger errors. Proteins were classified by the miRNA concentration at which repression reached the significance threshold of log 2 (protein ratio) Ͻ Ϫ0.3 (Fig.  3C), yielding 49 (36%), 10 (7%), 50 (37%), and 28 (20%) proteins repressed at 10, 25, 50, and 100 nM, respectively. Similar numbers of targets were responsive within the more sensitive (10 nM) and less sensitive (50 and 100 nM) groups, whereas fewer targets were responsive with intermediate sen- a miRanda predicted an additional 56 proteins as potential targets, which had imperfect base-pairing within seed-matching sequence. Of these, nine were down-regulated by miR-22. These miRanda-predicted targets were excluded from the number of true positives. sitivity (25 nM). Different dose response patterns were observed in the more sensitive group, where half of the proteins showed maximal repression at 10 nM, whereas the other half reached a plateau between 10 -25 nM followed by further repression above 25 nM (Fig. 3D). Thus, the analysis revealed an unexpected bimodality with respect to the dose response behavior of primary miRNA targets.
Primary targets classified with greater or less sensitivity to miRNA concentration were then examined for correlations with prediction scores. Targets with greater sensitivity (i.e. repressed at 10 nM miR-22) showed significantly higher P CT (p Ͻ 0.002, M-W) and context scores (p Ͻ 0.003) than targets with less sensitivity (repressed at 50 -100 nM) (Fig. 3E). Sensitive targets were also correlated with the net free energy for miRNA-target binding (⌬⌬G; p Ͻ 0.02), which was inversely correlated with the free energy of 3ЈUTR folding (⌬G open ; p Ͻ 0.003), and uncorrelated with the free energy of miRNA:mRNA binding (⌬G duplex ; p ϭ 0.30) (Fig. 3F). This reveals an energetic basis underlying bimodality, and indicates that Targetscan and PITA scores favor targets responding with higher sensitivity to exogeneous miRNA.
Variation in Kinetics of Target Repression-Pulsed SILAC (pSILAC) was used to examine the kinetics of primary target repression on a proteome-wide scale (27). Here, cell labeling began at the onset of miR-22 treatment so that only newly synthesized proteins were labeled (Fig. 4A) in contrast to SILAC, where proteins were labeled prior to treatment (Fig.  1A). Measurements between pSILAC versus parallel SILAC experiments were strongly correlated (r ϭ 0.96; Fig. 4B), indicating comparable abilities of the two methods to report protein changes. This argues that previously unexplained discrepancies in the literature that showed poor correlations (r ϭ 0.3) between seed-matched proteins responsive to miR-1 in HeLa cells (27,40), were not caused by the different labeling methods used in the two studies. However, the slope of the plot Ϸ1.2 (R 2 ϭ 0.89) indicating that pSILAC reports higher protein ratios compared with SILAC, which can be attributed to the lower background of newly synthesized proteins quantified by pSILAC.
We examined how rapidly proteins responded to miR-22 on a proteome-wide scale. Cells were transfected with 50 nM miR-22 or NT, then harvested at 12, 18, 24, 36, and 48 h. Among 4733 quantifiable proteins in total, 1842 (39%) were quantifiable at all five time points. Log 2 (protein ratio) measured at 12 versus 48 h were correlated with slope ϭ 0.22 (r ϭ 0.52), reflecting smaller responses at 12 h compared with 48 h (Fig. 4C). In contrast, log 2 (protein ratio) at 36 versus 48 h were correlated with slope ϭ 0.94 (r ϭ 0.95), indicating that the changes in newly synthesized proteins reached steady-state by 36 h. Thus, fewer than 5% of proteins (85 of 1842) responded at 12 h, whereas 25% (449 of 1842) responded significantly by 48 h. Of the 156 primary targets, 95 proteins were quantifiable across the time course, of which 21,20,9,24, and 21 reached the significance threshold at 12, 18, 24, 36, and 48 h, respectively (Fig. 4D). Thus, primary targets were distributed across the time course, with similar numbers responding at early (Ͻ18 h) versus late (Ͼ36 h) times, with half the number responding at intermediate times (24 h).
We asked whether differential protein turnover might explain the differences in kinetics, where targets undergoing faster turnover might be down-regulated more rapidly than those undergoing slower turnover. To address this, protein half-lives were measured proteome-wide by pulse-labeling cells for 24 h with heavy ArgϩLys. Assuming first-order decay, half-lives were quantified for 1804 proteins, with Յ 20% error in at least two of three biological replicate experiments (supplemental Table S2). All half-lives measured were long (median t 1/2 ϭ 112 h), and primary targets repressed early (12-18 h) indeed showed significantly shorter half-lives (median t 1/2 ϭ 107 h) than those repressed at later times (36 -48 h; median t 1/2 ϭ 140 h) ( Fig. 4E; p Ͻ 0.02, M-W). However, the median t 1/2 values for the primary targets in the two categories (107 h and 140 h) were too long to explain the difference in kinetics of protein repression on the timescale of 12-48 h, indicating that variations in protein turnover are not likely to be the main cause of the variable kinetics of target repression.
Targets repressed at early versus late times were examined for correlations with P CT , context score, and free energy. Interestingly, targets responding early showed significantly higher P CT (p Ͻ 0.0003, M-W) and context score (p Ͻ 0.0004, M-W) than targets responding later (Fig. 4F). Likewise, targets responding early showed significantly higher ⌬⌬G (p Ͻ 0.0002) and ⌬G duplex (p Ͻ 0.04) as well as lower ⌬G open (p Ͻ 0.0007) (Fig. 4G). Thus, parameters related to seed sequence matching and surrounding features, seed sequence conservation, and stability of miRNA-mRNA and 3ЈUTR secondary structure stability, all contribute toward distinguishing the varying kinetics of miRNA responses.
Next, we examined relationships between target sensitivity to miRNA, the kinetics of repression, and variations in log 2 (miR-22/NT) (Fig. 4H). Targets repressed at 10 nM showed higher protein fold changes as well as a preferential association with early response times. In contrast, targets repressed at 50 -100 nM responded at both early and late times, with comparable protein ratios across all time points. The results indicate that primary targets with highest sensitivity to miRNA concentration are down-regulated more rapidly, whereas the remaining targets show little correlation between sensitivity, kinetics, and fold-change.
Primary Targeted Proteins versus mRNA-Finally, we examined variations in repression at the mRNA level. Cells were transfected with 50 nM miR-22 or NT and harvested at 12, 18, 24, 36, and 48 h, followed by transcript sequencing and quantification by RNAseq. Among 12,023 genes quantifiable with FPKM greater than 1, as many as 11,088 were quantifiable across all five time points (Table S3). Weak correlations were observed between log 2 (mRNA ratio) measured at 12 versus 48 h (slope ϭ 0.39, r ϭ 0.25), in contrast to strong correlations between 36 versus 48 h (slope ϭ 0.85, r ϭ 0.75) (Fig. 5A). Thus, mRNA responses were small at 12 h and reached steady state by 36 h, as observed at the protein level. mRNAs with 3ЈUTR seed-matches were then selected and sorted based on the time at which they reached a significance threshold of log 2 (miR-22/NT) Ͻ Ϫ0.3 (Fig. 5B). In total, 617 transcripts containing 8mer or 7mer seed-matches were down-regulated, of which 283, 99, 47, 131, and 57 genes were repressed at 12, 18, 24, 36, and 48 h, respectively. Thus, like proteins, significant numbers of primary targets responded at early (12-18 h) or late (Ͼ36 h) times, but fewer responded at intermediate times (24 h).
Plots of mRNA versus protein ratios were correlated at all time points (r ϭ 0.43-0.76), with slope greater than 1 between 12-24 h and approximately equal to 1 at 36 -48 h (Fig. 5C). Thus, at earlier time points, mRNA ratios exceeded their corresponding protein ratios, whereas at steady state, the ratios closely matched. This suggests that genome-wide, mRNA repression generally precedes protein repression. Time courses for mRNAs and corresponding proteins were compared for primary targets, including 93 of the 95 protein targets shown in Fig. 4D that were quantifiable at the transcript level (Fig. 5D, supplemental Fig. S3). Of these, mRNA downregulation preceded or occurred concurrently with protein down-regulation in 87 (94%) cases (Fig. 5D, left and middle panels), whereas protein changes preceded mRNA changes in six (6%) cases (Fig. 5D, right panel). Overall, the majority of protein responses to miR-22 occurred concomitantly with responses of mRNA destabilization and decay, while few but finite numbers reflected predominant protein responses.
Next, we examined early versus late transcript responses for correlations with P CT , context score, and free energy. As observed with proteins, primary target mRNAs that responded early showed significantly higher P CT (p Ͻ 10 Ϫ6 ) and context score (p Ͻ 10 Ϫ6 ) (Fig. 5E). Like early proteins, early mRNA targets showed significantly higher ⌬⌬G (p Ͻ 3 ϫ 10 Ϫ6 ) and lower ⌬G open (p Ͻ 10 Ϫ6 ) compared with later targets (Fig. 5F), whereas ⌬G duplex showed little correlation (p ϭ 0.32). Taken together, global effects on mRNA measured by RNAseq recapitulated the kinetic behavior of proteins measured by the proteome-wide screen, with similar physical features distinguishing early versus late targets.
Analyses of miRNA reporters have suggested that target repression depends on a balance between cellular concentrations of miRNA and mRNA, where stronger repression occurs when mRNA targets are expressed at low abundance (12). We tested the hypothesis that targets responding with higher sensitivity and/or at earlier times might reflect lower copy number transcripts where higher repression efficacy might be expected. Transcript abundances were quantified by FKPM and correlated against primary targets that varied with respect to kinetics and dose responses. For targets responding at early times, the median transcript abundance was higher than for targets responding at later times (Fig. 5G, p ϭ 0.07). Likewise, for targets responding with higher sensitivity to miRNA, the median transcript abundance was greater than targets responding with less sensitivity (Fig. 5G, p ϭ 0.07). Therefore, the trends in each case reflected higher transcript abundance among the targets with greater sensitivity to miRNA and earlier responses, contradicting our initial prediction. DISCUSSION Our study represents the first proteome-wide analysis of primary miRNA responses as a function of treatment time and concentration. Primary targets were identified based on the criteria that 1) they are down-regulated in response to miR-22, and 2) their mRNA 3ЈUTRs contain 8mer or 7mer seedmatches. We observe significant and unexpectedly large dif-ferences in the response thresholds of primary targets with respect to time and concentration. A key finding is that variations in threshold appear bimodal, with greater numbers of proteins repressed at low (10 nM) and high (50 -100 nM) exogenous miRNA, compared with an intermediate (25 nM) concentration. Primary targets also show variable thresholds for the onset of repression, showing greater numbers with rapid (12-18 h) or slow (36 -48 h) kinetics, compared with intermediate (24 h) times. Similar behavior is seen at the transcript level. Proteins responding to low miRNA concentration also showed two distinct patterns, one where repression was maximal at 10 nM, and another where repression plateaued at 10 -25 nM and then increased further above 50 nM. Our results expand the understanding of miRNA target repression, by revealing heterogeneity in mechanisms that reflect distinct features of transcript sequences, and modulate thresholds in sensitivity and timing of protein responses.
Although the explicit mechanisms are unknown, certain possibilities can be eliminated. Our initial hypothesis was that rapid versus delayed repression may reflect proteins with shorter versus longer half-lives, respectively. This was clearly not supported by pulse-chase proteomics experiments, because nearly all proteins showed half-lives that were longer than the early versus late timescale of repression. We conclude that variable kinetics are better explained by differences in rates of translational repression and/or mRNA decay. Likewise, transcript abundance measurements lent no support to the hypothesis that transcripts found in low copy number are more sensitive to miRNAs. Furthermore, no bias was observed either in the number or the type of seed sequence sites, when comparing high versus low sensitivity, or fast versus slow responses (supplemental Fig. S4). This argues against models where repression mechanisms differ by the number of miRISC complexes that can bind to a single target transcript.
Both the kinetics and dose responses of protein and mRNA changes turned out to be most strongly correlated with ⌬⌬G, the net free energy of miRNA-target interactions, which mainly reflected inverse correlations with ⌬G open , the free energy of forming secondary structure at or nearby the miR-22 seed sequence sites. Thus, our working model is that variability somehow reflects stabilization of secondary structures within the 3ЈUTR, which in turn interferes with the accessibility of miRISC to seed matching sequences, blocking miRNA-mRNA hybridization. The heterogeneity in concentration dependence between different transcripts means that relevant binding sites are not all saturated at low concentrations of miR-22. More sensitive targets -i.e. those responding to 10 nM miR-22responded with faster kinetics, suggesting that one mode of repression reflects facile interactions with target sites and rapid responses. In contrast, the less sensitive targets are equally distributed between fast and slow kinetics, suggesting that energetic barriers to unfolding shift the threshold for saturation to higher concentrations, but are not rate-limiting.
Varying thresholds for saturation imply mechanisms involving competition for miRNA binding sites. For example, proteins binding to 3ЈUTR regulatory elements at or near the seed sequence site might compete for miRISC binding. This has been observed in several systems where RNA binding proteins (RNPs) compete for miRNA binding sites (49 -52). Alternatively, proteins may indirectly perturb transcript secondary structure, as in the case of HuR, an RNP which, in response to cell stress, reverses miRNA-mediated repression by oligomerizing in complex with 3ЈUTR sequences at a distance from the seed sequence site (53,54). Thus, the observed gap in concentration threshold might be explained if RNP interference is weak toward targets with high miRNA-mRNA affinity, leading to faster and more sensitive responses to miRNAs, but efficient toward targets with lower affinity, requiring higher miRNA concentrations or longer times, enough to compete for binding. Other mechanisms involving competing miRNA sites within the cell might also create gaps by increasing the threshold for miRISC binding to lower affinity transcripts (55,56). Our findings provide evidence for the potential involvement of cellular factors that separate miRNA targets into distinct classes with respect to timing and sensitivity.
Intriguingly, we found a set of transcripts that were downregulated at 10 nM, reached a plateau between 10 -25 nM, and were further down-regulated above 50 nM miRNA. This might be explained by the presence of multiple sites within 3ЈUTRs with varying affinities, which lead to different concentration thresholds for the entire population of each transcript. However, inspection of these transcripts showed no differences from others in numbers or types of seed sequence sites (not shown), arguing against this model. A more likely explanation is that certain transcripts are heterogeneous within each population, with respect to their ability to function as miRNA targets. This could reflect transcript isoforms such as alternative splice variants, or else variations in structure that lead to heterogeneity in RNP interactions and interference with binding at low miRNA concentrations. The finding predicts that certain transcripts will show graded responses in sensitivity and fold-change, under cellular conditions that lead to varying steady state miRNA levels.
The mechanism for miRNA repression involves miRISC-GW182-PABC recruitment of enzymes catalyzing deadenylation, 5Јdecapping, and exonuclease-catalyzed degradation of mRNA, as well as sequestration of eukaryotic initiation factors leading to translational repression (7,57). Current models cannot readily distinguish between mechanisms involving protein loss primarily caused by mRNA decay, versus situations where mRNA decay and translational repression occur independently. In principle, both processes could be strongly coupled, whereby triggering the same repression/decay mechanism induces both processes, which are biased toward loss of protein versus mRNA depending on cellular conditions. Thus, it was interesting that we observed three types of patterns when comparing the kinetics of mRNA versus protein repression. The overwhelmingly largest numbers of cases (77, 84%) were consistent with strong coupling between mRNA and protein regulation, showing mRNA decay occurring concomitantly with protein down-regulation. Only 10 cases showed mRNA preceding protein down-regulation (BTF3, CLCN7, EPB41L2, GBA, GNPDA1, MAP1B, PTGFRN, PURB, STX4, and TACC1, supplemental Fig. S3), while two cases showed mRNA following protein down-regulation (MYO5A and RBM17, supplemental Fig. S3). In four cases, proteins were down-regulated with no evidence for loss of mRNA (GLG1, GNAS, SPARC, and SYNE2, Fig. 5D, supplemental Fig. S3). Our findings reveal heterogeneity in protein versus mRNA repression mechanisms, which may be useful for identifying features of transcripts that influence the repression mechanism. Our findings differ from those in recent studies in Drosophila S2 cells and zebrafish embryos, where translational repression appears to be the dominant mechanism for miRNA-mediated protein repression (58,59). This could be caused by species differences and/or differences in timing (e.g. 6 -12 h in zebrafish). In any case, our observations suggest that once mRNA degradation starts, it becomes the dominant mechanism in mammalian cells. The primary target set revealed substantive differences between prediction algorithms (Table I). On one hand, miRanda predicted 259 targets, of which 63% were false positives and 44% were false negatives. On the other hand, TargetScan predicted only 69 targets, which decreased the false positives to 48% but increased false negatives to 77%. Likewise, PITA predicted 52 targets, yielding 46% false positives and 82% false negatives. The proteins identified by Targetscan overlapped those identified by PITA by ϳ50%, indicating greater similarity in performance between these algorithms compared with miRanda, where the overlap was only ϳ25% (supplemental Fig. S1). This reveals partial interdependence between sequence context and evolutionary conservation features evaluated by Targetscan and the free energy of miRNA-mRNA hybridization evaluated by PITA.
As in previous reports (27, 40 -44), the responses to miRNA reflected mild repression at the protein level, ranging from 1.25-to 4-fold. The magnitude of protein fold-changes, which is often used to evaluate efficacy, correlated well with seed matching sequence classification, P CT , and context scores. To some extent, such correlations reflect biases in the methods used to develop the scores. For example, context scores were developed from transcript microarray measurements assayed 12-24 h post-transfection (19), and therefore might be expected to favor targets that respond earlier, whereas excluding primary targets that respond more slowly. Likewise, PITA is optimized based on luciferase expression responses to endogenous miRNAs, measured 20 h after transfection of reporter plasmid. The bias in prediction programs toward early responses could be observed within the set of 156 primary targets, which showed that although the majority of Targetscan-or PITA-predicted targets mostly responded at early times (12-24 h), targets missed by these programs (false negatives) were more equally represented at early and late times (supplement Fig. S5). Proteins may show larger foldchange because they reflect targets that respond earlier or with higher sensitivity, whereas additional repression might be observable if experimental conditions were extended to longer times and higher concentrations. Viewed in this way, variations in fold-changes may partially reflect the heterogeneity in mechanisms affecting kinetics and dose response among primary targets. An implication of our study is that important targets that respond late or at high concentrations may be underestimated when efficacy is evaluated based on fold-change or scores that are correlated with magnitude of responses.
In summary, our study documents heterogeneity in primary targets of a single miRNA, predicting that in a normal biological context, waves of target responses may occur with respect to kinetics and concentration that vary depending on cellular state. Estimates of miRNA half-life range from 1 to Ͼ200 h in animal cells (60 -62), arguing that targets with extended kinetics are likely to be biologically relevant. In addition, measurements of intracellular miRNA concentrations show 60-to 250-fold increases in response to transfection of 15-50 nM exogenous miRNA (39,63,64), matching the more than 1000-fold changes in endogenous miRNAs quantified under varying cellular (e.g. disease) conditions (65,66). This argues that miRNA responses cover a dynamic range that includes a wide range of exogenous miRNA concentrations. We propose that primary targets that respond late and with lower sensitivity, although often excluded from analysis, are important to consider as physiological regulators of miRNA-controlled processes. It will be important to compare our findings with miR-22 against other miRNAs, to see whether miRNAs behave differently with respect to dose-response and kinetics of target repression.