Introduction MicroRNAs (miRNAs) are small noncoding RNAs whose complementary pairing to target mRNAs potentially regulates expression of more than 60% of genes in many and perhaps all metazoans [1]–[6]. Destabilization of mRNA and translational repression have been suggested as the mechanisms of action for miRNAs [1],[3],[7]–[15], and recent work directly measuring endogenous protein levels in response to altered miRNA expression levels found that specific miRNAs modestly inhibit the production of hundreds of proteins [16],[17]. The importance and functional range of miRNAs are evidenced by the diverse and often dramatic phenotypic consequences when miRNAs are mutated or misexpressed, leading to aberrant development or disease [7],[18]–[24]. Although regulation by miRNAs is an integral component of the global gene expression program, there is currently no consensus on either the mechanism by which they decrease the levels of the targeted proteins or even the steps in gene expression regulated by miRNAs [3],[25]–[29]. The proposal that miRNAs decrease protein levels without affecting mRNA stability arose from the observation that the miRNA lin-4 down-regulates lin-14 expression in the absence of noticeable changes in lin-14 mRNA abundance in Caenorhabditis elegans [7],[30]–[33]. Subsequent studies in mammalian cell culture provided further support for this model [34]–[37]. Several studies have found that repressed mRNAs as well as protein components of the miRNA regulatory system accumulate in P-bodies, suggesting that repressed mRNAs may be sequestered away from the translation pool [38]–[44]. Other evidence points to deadenylation of miRNA-targeted mRNAs, an effect that can inhibit translation [14],[45]–[53]. Some studies have argued that initiation of translation is blocked at either an early, cap-dependent stage or later during AUG recognition or 60S joining [10],[44],[52],[54]–[58]. Others have argued that a postinitiation step is targeted, resulting in either slowed elongation, ribosome drop-off, or nascent polypeptide degradation [7],[59]–[62]. One factor contributing to the lack of a consensus model for miRNA function is the evidence that miRNA targeting of an mRNA significantly reduces message levels (despite previous reports to the contrary) [9],[11],[12],[14],[52],[63],[64]. Indeed, very recent studies from Baek et al. and Selbach et al. found that the changes in mRNA abundance are not only correlated with the repression of many targets, but also can account for most of the observed reduction in protein expression [16],[17]. mRNA targets of the same miRNA can either be translationally repressed with little change in mRNA abundance, translationally repressed and have concordant changes in mRNA abundance, or have little translation repression with large changes in mRNA abundance [52],[65],[66]. That miRNAs can affect both protein production and abundance of their mRNA targets raises the question of to what extent these outcomes of miRNA regulation are mediated by a common mechanism or by competing or complementary processes. The regulatory consequence of a particular miRNA–mRNA interaction might be influenced by miRNA-independent factors such as cellular context or by additional information encoded by the target mRNA, e.g., presence of binding sites for other RNA-binding proteins and miRNAs, secondary structure around miRNA binding sites, or the intrinsic decay rate of the mRNA [25],[51],[67],[68]. Experiment-specific effects of in vitro translation assays, reporter constructs, or procedural differences that alter properties of gene expression could account for some of the wide variation in the apparent mechanisms by which miRNAs alter expression [25],[27]. To date, most studies on translational regulation by miRNAs have used reporter assays. Although assays that rely on engineered reporter transcripts are powerful, assay-specific anomalies are a concern; artificial mRNAs may lack key pieces of regulatory information, overexpression of reporter mRNAs could mask subtle regulatory functions, and DNA transfection can lead to indirect effects on cell physiology [26]. Indeed, recent reports have found that differences in experimental setup, such as the method of transfection, type of 5′-cap, or the promoter sequence of the DNA reporter construct can drastically alter the degree or even the apparent mode of regulation by miRNAs [59],[69]. In addition, some models have been based on studies in which only one or a few targets were studied, which introduces the possibility of generalizing the behavior of a single miRNA–mRNA interaction that may not represent the dominant biological mechanism. Two recent studies avoided many of these caveats by overexpressing, inhibiting, or deleting specific miRNAs and systematically measuring changes in endogenous mRNA and protein levels using DNA microarrays and stable isotope labeling with amino acids in cell culture (SILAC), respectively [16],[17]. Both studies found mostly concordant changes in mRNA levels and protein levels, with changes in mRNA levels accounting for much, but not all, of the changes in protein abundance. With data for hundreds of endogenous targets, these studies were the first to provide genome-wide evidence that mRNA degradation accounts for much of the reduction in protein levels. And whereas these results suggest that translation inhibition accounts for some of the observed changes in protein abundance of miRNA targets, they do not provide direct evidence of this, nor do they provide insight into which steps in translation are regulated, the extent this regulation contributes to reduced gene expression of specific mRNAs, or its possible links to mRNA decay. To investigate how miRNAs regulate gene expression, we systematically identified direct targets of the miRNA miR-124 by measuring the recruitment of target mRNAs to Argonaute (Ago) proteins, the core components of the miRNA effector complex, as previously described [70]–[72]. We then measured, in parallel, mRNA abundance and two indicators of translation rate, ribosome occupancy and ribosome density, for more than 8,000 genes, using DNA microarrays and a novel polysome encoding scheme. This strategy allowed us to directly investigate the behavior of miRNA–mRNA target pairs with respect to both mRNA fate and translation, on a genomic scale. Results Systematic Identification of mRNAs Recruited to Argonautes by miR-124 To study the effects of miR-124 on expression of mRNA targets, we first had to identify those targets. Recruitment to Ago complexes in response to the expression of a particular miRNA appears to be the most reliable criterion for target identification [70]. To this end, we lysed human embryonic kidney (HEK) 293T cells transfected with miR-124 and isolated Ago-associated RNA by immunopurification (IP) using a monoclonal antibody that recognizes all four human Ago paralogs [73]. We measured mRNA enrichment in Ago IPs by comparative DNA microarray hybridization of samples prepared from immunupurified RNA and total RNA from cell extracts. Three replicates of Ago and control IPs were performed from both miR-124 and mock-transfected cells (Datasets S1 and S5). To examine the enrichment profiles of the IPs, we first clustered the microarray results by their similarity and visualized the results as a heatmap, with the degree of enrichment of each RNA shown on a green (least enriched) to red (most enriched) scale (Figure S1). The Ago IP enrichment profiles were reproducible as evidenced by an average Pearson correlation coefficient between mRNA enrichment profiles of Ago IPs in mock-transfected cells and miR-124–transfected cells of 0.90 and 0.94, respectively. Thousands of mRNAs were reproducibly enriched in the Ago IPs from mock-transfected cells (Figures S1 and S2, and Text S1). We found that the presence of sequence matches to two highly expressed microRNA families, miR-17-5p/20/92/106/591.d and miR-19a/b, in the 3′-untranslated regions (UTRs) of mRNAs significantly correlated with Ago IP enrichment (Text S2), suggesting that association with Ago is in large part a reflection of the relative occupancy of each mRNA with the suite of miRNAs endogenously expressed in HEK293T cells. High-confidence Ago-associated mRNAs (at least 4-fold enriched over the mean, 1,363 mRNAs) disproportionately encode regulatory proteins (409, p = 0.001), with roles including “transcription factor activity” (95, p = 0.01), “signal transduction” (230, p = 0.02), and “gene silencing by RNA” (7, p = 0.02). To identify RNAs specifically recruited to Agos by miR-124, we compared the mRNA enrichment profiles of Ago IPs from miR-124–transfected cells to Ago IPs from mock-transfected cells using the significance analysis of microarrays (SAM) modified two-sample unpaired t-test (Datasets S1 and S5). At a stringent 1% local false-discovery rate (FDR) threshold, we identified 623 distinct mRNAs significantly enriched in Ago IPs from lysates of miR-124–transfected cells compared to Ago IPs from mock-transfected cells (Figure 1A). 10.1371/journal.pbio.1000238.g001 Figure 1 miR-124 recruits hundreds of specific mRNAs to Argonautes. (A) Supervised hierarchical clustering of the enrichment profiles of putative miR-124 Ago IP targets (1% local FDR) in Ago IPs from miR-124–transfected cells (blue) and mock-transfected cells (black). Rows correspond to 789 sequences (representing 623 genomic loci with a RefSeq sequence), and columns represent individual experiments. The color scale encompasses a range from 0.06- to 16-fold relative to global mean IP enrichment for each mRNA (−4 to +4 logs on log base 2 scale ). (B) Enrichment of seed matches to miR-124 in the 3′-UTRs and coding sequences of miR-124 Ago IP targets (1% local FDR). The significance of enrichment of seed matches in Ago IP targets was measured using the hypergeometric distribution function. Previous work established that the 5′-end of the miRNA, the “seed region,” is particularly important for interactions with mRNA targets [4],[11],[37],[74]–[77]. In most cases, there is a 6–8 bp stretch of perfect complementarity between the seed region of the miRNA and a “seed match” sequence in the 3′-UTR of the mRNA [4],[11],[37],[74]–[77]. We reasoned that if the mRNAs specifically recruited to Agos by miR-124 transfection were physically associated with miR-124, seed match sequences would be significantly enriched in miR-124–specific IP targets compared to nontargets. Indeed, we found strong enrichment of 6–8 base seed matches to miR-124 in the 3′-UTRs of miR-124 Ago IP targets (Figure 1B). We also found enrichment within the coding sequences of miR-124 Ago IP targets, as previously reported (Figure 1B) [11],[16],[17],[70],[71],[78],[79]. For instance, 60% of miR-124 Ago IP targets contain a perfect match to positions 2–8 of miR-124 (called 7mer-m8) in their 3′-UTRs, compared to 10% of nontargets (p 0.001) to miR-124 IP targets that change less than 40% in mRNA abundance. The red curve shows a normal distribution with mean (0.14) and standard deviation (0.04) from the 10,000 permuted sets. The red arrow shows the Pearson correlation of miR-124 IP targets that change less than 40% in mRNA abundance (r = 0.30, p<10−5). (0.39 MB PDF) Click here for additional data file. Figure S8 Concordant changes in mRNA abundance and translation of miR-124 Ago IP targets with 7mer 3′-UTR seed matches and miR-124 Ago IP targets that lack a 7mer 3′-UTR seed match. (A) Scatterplot between changes in mRNA abundance (x-axis) and the estimated translation rate (y-axis) for miR-124 Ago IP targets with 7mer 3′-UTR seed matches following transfection with miR-124 compared to mock. The gray line is a least-squares linear regression fit of the data, and the black line is a moving average plot (window of 10). The slope of the least-squares fit of the data = 0.23 (in linear space = 0.37) and the Pearson correlation = 0.59. (B) Scatterplot between changes in mRNA abundance (x-axis) and the estimated translation rate (y-axis) for miR-124 Ago IP targets that lack 7mer 3′-UTR seed matches following transfection with miR-124 compared to mock. The gray line is a least-squares linear regression fit of the data, and the red line is a moving average plot (window of 10). The slope of the least-squares fit of the data = 0.21 (in linear space = 0.24) and the Pearson correlation = 0.42. (0.35 MB PDF) Click here for additional data file. Figure S9 Changes in abundance and translation of miR-124 Ago IP targets with seed matches in 3′-UTRs, coding sequences and 5′-UTRs. (A) Cumulative distribution of the change in mRNA levels following transfection with miR-124 compared to mock. This analysis compares miR-124 Ago IP targets (1% local FDR) with at least one 3′-UTR 7mer seed match, but no coding sequence or 5′-UTR 7mer seed matches (red, 244), IP targets with at least one 3′-UTR 6mer seed match (green, 47), but no 3′-UTR, coding sequence, or 5′-UTR 7mer seed matches, IP targets with at least one coding sequence 7mer seed match, but no 3′-UTR or 5′-UTR 7mer seed matches (blue, 70), IP targets that lacked a 6mer seed match in the 3′-UTR, coding sequence, or 5′-UTR (orange,23), and nontargets (7385, black). This analysis compares Ago IP targets (red) versus nontargets (black). (B) Cumulative distribution of the change in translation following transfection with miR-124 compared to mock. This analysis compares miR-124 Ago IP targets (1% local FDR) with at least one 3′-UTR 7mer seed match, but no coding sequence or 5′-UTR 7mer seed matches (red), IP targets with at least one 3′-UTR 6mer seed match (green), but no 7mer seed matches in the 3′-UTR, coding sequence, or 5′-UTR, IP targets with at least one coding sequence 7mer seed match, but no 7mer seed match in the 3′-UTR or 5′-UTR (blue), IP targets that lacked a 6mer seed match in the 3′-UTR, coding sequence, or 5′-UTR (orange), and nontargets (black). This analysis compares Ago IP targets (red) versus nontargets (black). (C) Bar plot of the average change in mRNA abundance (blue) and translation rate (red) of miR-124 Ago IP targets following transfection with miR-124. The average change in mRNA abundance and translation of targets was calculated by subtracting the average change of nontargets for the mRNA abundance and translation rate measurements following transfection with miR-124. This analysis compares miR-124 Ago IP targets (1% local FDR) with at least one 3′-UTR 7mer seed match, but no coding sequence or 5′-UTR 7mer seed matches, IP targets with at least one 3′-UTR 6mer seed match, but no 7mer seed matches in the 3′-UTR, coding sequence, or 5′-UTR, IP targets with at least one coding sequence 7mer seed match, but no 7mer seed match in the 3′-UTR or 5′-UTR, IP targets with at least one 7mer seed match in the 5′-UTR, but no 7mer seed match in the 3′-UTR or coding sequence, and IP targets that lacked a 6mer seed match in the 3′-UTR, coding sequence, or 5′-UTR. This analysis compares Ago IP targets (red) versus nontargets (black). The error bars represent 95% confidence intervals in the mean difference estimated by bootstrap analysis. (0.65 MB PDF) Click here for additional data file. Figure S10 Efficiency of recruitment to Argonautes by miR-124 seed matches correlates with effects on both mRNA abundance and translation. (A) Scatterplot between changes in Ago IP enrichment (x-axis) following transfection with miR-124 compared to mock and estimated changes in protein production (Equation 2) (y-axis) for mRNAs with either 8mer seed matches (red dots) or 7mer seed matches (blue dots) to miR-124 in their 3′-UTRs. For 8mer seed matches, the slope of the least-squares fit of the log2 data is −0.46 (in linear space, −0.03), and the log2 Pearson correlation is −0.72. For 7mer seed matches, the slope of the least-squares fit of the log2 data is −0.39 (in linear space, −0.05), and the log2 Pearson correlation is −0.72. (B) Same as (A) except for mRNAs with seed matches to miR-124 in their coding sequences and no 7mer seed matches in their 3′-UTRs. For 8mer seed matches, the slope of the least-squares fit of the log2 data is −0.13 (in linear space, −0.008), and the log2 Pearson correlation is −0.39. For 7mer seed matches, the slope of the least-squares fit of the log2 data = −0.14 (in linear space, −0.02), and the log2 Pearson correlation is 0.38. (C) Same as in (A) except for mRNAs with 6mer seed matches to miR-124 in their 3′-UTR, but no 7mer seed match in their 3′-UTR or coding sequence (red), and mRNAs that lack 6mer seed matches in their 3′-UTR or coding sequence (blue). For 6mer seed matches, the slope of the least-squares fit of the log2 data is −0.21 (in linear space, −0.09), and the log2 Pearson correlation is −0.40. For mRNAs without 6mer seed matches, the slope of the least-squares fit of the log2 data is −0.13 (in linear space, −0.07), and the log2 Pearson correlation is −0.22. (1.71 MB PDF) Click here for additional data file. Table S1 Summary of miR-124 targets for Western blot analysis. (0.02 MB PDF) Click here for additional data file. Table S2 Exogenous doping control information. (0.05 MB XLS) Click here for additional data file. Text S1 miRNA-effector complexes appear to nonspecifically bind streptavidin-coated Dynal beads. (0.03 MB DOC) Click here for additional data file. Text S2 Enrichment of seed matches to highly expressed miRNAs in Ago IPs from mock-transfected cells. (0.03 MB DOC) Click here for additional data file. Text S3 Relationship between ribosome occupancy in mock-transfected cells and changes in ribosome occupancy following transfection of miR-124. (0.03 MB DOC) Click here for additional data file. Text S4 Evaluation of the significance of the correlation between changes in mRNA abundance and translation of miR-124 Ago IP targets following transfection with miR-124. (0.05 MB DOC) Click here for additional data file.
