The regulatory impact of RNA-binding proteins on microRNA targeting

Argonaute is the primary mediator of metazoan miRNA targeting (MT). Among the currently identified >1,500 human RNA-binding proteins (RBPs), there are only a handful of RBPs known to enhance MT and several others reported to suppress MT, leaving the global impact of RBPs on MT elusive. In this study, we have systematically analyzed transcriptome-wide binding sites for 150 human RBPs and evaluated the quantitative effect of individual RBPs on MT efficacy. In contrast to previous studies, we show that most RBPs significantly affect MT and that all of those MT-regulating RBPs function as MT enhancers rather than suppressors, by making the local secondary structure of the target site accessible to Argonaute. Our findings illuminate the unappreciated regulatory impact of human RBPs on MT, and as these RBPs may play key roles in the gene regulatory network governed by metazoan miRNAs, MT should be understood in the context of co-regulating RBPs.

I t is broadly accepted that metazoan miRNA targeting (MT) is governed by a ternary interaction of Argonaute (AGO), miRNA, and its mRNA target [1][2][3][4] . While AGO is the main mediator of MT, there is an increasing amount of evidence that other RNA-binding proteins (RBPs) also play important regulatory roles in MT 5 . For instance, phosphorylated Pumilio protein opens up the local hairpin structure of the CDKN1B 3′ UTR, where the miR-221/222 target site is otherwise inaccessible, and lead to productive MT by miR-221/222 6 . There are a few other examples of such an MT enhancer as PCBP2 and FUS 7,8 . In contrast, Dnd1, RBM38, and IGF2BP1 have been reported to function as MT suppressors where these RBPs suppress MT mostly by making the miRNA target site more inaccessible to AGO [9][10][11] . HuR and PTBP have been reported to be enhancers of MT by either recruiting AGO or opening the secondary structure to increase accessibility to AGO as well as suppressors of MT by competing against AGO [12][13][14] .
On the other hand, >1,500 human RBPs have been identified 15 and each RBP is estimated to have on average 22,000 3′UTRbinding sites 16 , leading to >33 million interactions that can occur between human RBPs and 3′UTRs (see below). Despite the enormous number of possible interactions, only a handful of aforementioned interactions have been examined so far, leaving almost all other interactions unexamined. Accordingly, our current understanding towards the global regulatory impact of RBPs on MT remains severely limited.
In this study, by analyzing the transcriptome-wide binding sites for 150 human RBPs and large-scale datasets that monitored the whole-transcriptome response to ectopically introduced or deleted miRNAs, we attempted to systematically evaluate the quantitative effect of RBPs on MT and thus to help gain a comprehensive insight into the gene regulatory network of metazoan miRNAs and their co-regulating RBPs.

Results
RBPs have a large number of 3′UTR-binding sites. To accurately detect binding sites of RBPs, CLIP-seq (crosslinking and immunoprecipitation followed by sequencing) has been developed [17][18][19][20][21][22] . Taking advantage of this powerful technology, the ENCODE consortium has published a massive-scale dataset of an enhanced version, termed eCLIP-seq 16,23 . We obtained and analyzed the ENCODE eCLIP-seq dataset to identify the transcriptome-wide binding sites for 150 RBPs profiled in HepG2 and K562 cell lines.
Our analysis indicated that human RBPs not only bind to the 5′UTRs and coding regions, but also to the 3′UTRs substantially: the evaluated 150 RBPs have on average 22,000 3′UTRbinding sites ( Supplementary Fig. 1a) ranging from 1,000 to 73,000. Extrapolating the average number of 3′UTR-binding sites to 1,500 human RBPs, we estimate that >33 million interactions, several orders of magnitude larger than the previously evaluated interactions, can occur between human RBPs and 3′UTRs.
Because metazoan MT occurs primarily in the cytoplasm 24,25 , we quantitatively assessed the subcellular localization of RBPs by analyzing immuno-fluorescence images 26 . When measuring the cytoplasmic fraction compared to the nucleus fraction for each RBP, 91% of the evaluated RBPs exhibited >5% of cytoplasmic fraction ( Supplementary Fig. 1b), suggesting that almost all RBPs are localized in the cytoplasm to a detectable degree. These findings demonstrate that human RBPs have a large number of 3′ UTR-binding sites and also localized in the cytoplasm despite the remarkably varying fraction, justifying our hypothesis to evaluate whether RBPs may globally influence MT.
Strong association between RBP binding and enhanced miRNA targeting. To examine the association between RBP binding and MT, we have generated a large dataset of mRNA-seq that measured the whole-transcriptome response to overexpressed miRNAs in HepG2 cell line. By combining this dataset with the ENCODE RBP-binding site (RBS) dataset, we collected target mRNAs with a single 7, 8mer miRNA target site (MTS) and tested whether the distance between an MTS and the nearest RBS on the 3′UTR, denoted as d MTS-RBS , is correlated with MT efficacy (Fig. 1a, Supplementary Discussion), hypothesizing that the RBS located close to an MTS might influence MT. Although we examined whether d MTS-RBS is associated with MT efficacy either positively or negatively to potentially discover both MT enhancers and suppressors, the shorter d MTS-RBS was significantly correlated with the improved MT efficacy only (Fig. 1b, HepG2  panel). This association was significant even after correcting for known confounding factors (local AU content, 3′UTR size, target-site abundance, and seed-pairing stability) [27][28][29][30][31][32][33][34] by the multiple linear regression (Fig. 1c, HepG2 panel). To rule out these confounding effects more definitively, we selected a group of 3′UTRs that have different d MTS-RBS but have statistically indistinguishable confounding factors. Even after such a rigorous correction, our observation was still consistent and we were able to confirm the independent correlation between d MTS-RBS and MT efficacy on a global scale (Fig. 1d, HepG2 panel). When compared with known determinants of MT such as target-site abundance 29 , the overall impact of d MTS-RBS on MT was significantly stronger than those of previously reported determinants, emphasizing its potential role as a key determinant of MT ( Supplementary  Fig. 1c).
To confirm our observation is not limited to HepG2 cell line, we analyzed a large dataset that monitored the wholetranscriptome response to overexpressed miRNAs in various other cell lines 29 (Fig. 1b-d). For the RBS information of these cell lines, we used the ENCODE RBS information obtained from both HepG2 and K562 cell lines, since we confirmed that RBSs are robust enough to be preserved between these two cell lines and thus these RBSs can also be applied to other cell lines ( Supplementary Fig. 1d, e). Accordingly, we again observed a significant correlation between d MTS-RBS and MT efficacy with other independent datasets (Fig. 1b-d and Supplementary Fig. 1f) and with the preserved binding sites ( Supplementary Fig. 1g), supporting that the observed association is general enough to be extended in various biological contexts.
However, given that these results are based on the transcriptome response to ectopically introduced miRNAs, it was crucial to confirm the results in an endogenous condition as well. To do so, we used a dataset of DROSHA and DICER knockout cell lines where endogenous miRNAs are depleted 35 . In accord with our previous results, we observed a significant stronger derepression of mRNAs with shorter d MTS-RBS ( Fig. 1e and Supplementary  Fig. 1h), demonstrating that the observed effect of RBPs on MT is a phenomenon occurring in an endogenous environment.
Most RBPs are associated with enhanced miRNA targeting. To examine if this correlation between d MTS-RBS and MT efficacy is observed with individual RBPs, we performed a correlation analysis for each individual RBP after correcting for confounding factors in various cell lines (Fig. 2a-c). First, 86%, 93%, and 94% of the RBPs evaluated in HepG2, HeLa, and other human cancer cell lines, respectively, exhibited significant correlations between d MTS-RBS and MT efficacy. These results illustrate that the regulatory impact of RBPs on MT may be broad in contrast to the previously reported examples 6-8,10-12,36 . Second, similar to Fig. 1b, in 100% of these RBPs that showed significant correlation, the shorter d MTS-RBS was associated with stronger MT efficacy indicating that all of these RBPs function as MT enhancers rather than suppressors; this is another striking inconsistency with the previously reported instances of RBPs that function as MT suppressors [10][11][12]36 . Taken together, these results demonstrate that for most RBPs, if not all, their binding close to the MTS is associated with enhanced MT efficacy, while no RBP is detectably associated with suppression of MT on a global scale.
To evaluate the collective effect of RBPs on MT, we investigated the association between MT efficacy and the number of bound RBPs in close proximity to MTSs. As a result, the number of RBPs was positively correlated with improved MT efficacy ( Fig. 3a and Supplementary Fig. 2a, b), even when a small number of RBPs are bound near the MTS ( Supplementary  Fig. 2c). The correlation was also consistently observed for the MTSs with the preserved binding sites (Supplementary Fig. 2d) and regardless of whether the RBP has specific binding motifs 37 or not ( Supplementary Fig. 2e), suggesting that the observed impact of RBPs is quite general. To further examine whether the specific identities of RBPs instead of the number of bound RBPs determine MT efficacy, we have partitioned MTSs into two subgroups with similar numbers of the bound RBPs but with different identities of the RBPs (Fig. 3b, c and Supplementary  Figs. 3, 4). Consequently, two subgroups did not show a significant difference in MT efficacy ( Fig. 3c and Supplementary Figs. 3b, 4b, e). Therefore, instead of specific identities of RBPs, the overall number of bound RBPs to MTSs appears to be the primary factor that is associated with MT efficacy.
A potential mechanism by which RBP binding enhances miRNA targeting. One of the possible mechanisms that could explain the impact of RBPs on MT is protein-protein interaction between AGO and RBPs, as some RBPs have been previously Fig. 1 RBP binding close to a miRNA target site is associated with enhanced miRNA targeting. a Overview of the analysis to investigate the effect of RBP binding on the miRNA targeting (MT) efficacy. Genes were grouped together based on the distance between a miRNA target site (MTS) and the nearest RBP-binding site (RBS) on the 3′UTR, denoted as d MTS-RBS . The genes with an MTS that overlaps with an RBS were categorized as 'overlapped'. The association between MT efficacy and d MTS-RBS was analyzed by measuring the mRNA fold change. b Association analysis between d MTS-RBS and MT efficacy. The mean mRNA fold changes of 3′UTRs with a single 7, 8mer MTS obtained from mRNA-seq (HepG2) or microarray datasets (HeLa and other human cancer cell lines) that monitored the whole-transcriptome response to overexpressed miRNAs/siRNAs are shown. The 3′UTRs were grouped with respect to d MTS-RBS , as depicted in (a). The mRNA fold changes were compared between the 3′UTR groups (two-sided Wilcoxon's rank-sum test). The number of 3′ UTRs used for measuring d MTS-RBS is shown on top. mRNA fold changes are displayed in the log2 scale and the error bars represent 95% confidence intervals. c Residuals of the mRNA fold change for each group from (b) after correcting for four known confounding features of MT (local AU content, target abundance (TA), seed-pairing stability (SPS), and 3′UTR length). The regression residuals represent the remaining information after reducing the contribution from the confounding features. Otherwise as in (b). d Association analysis between d MTS-RBS and MT efficacy after a more rigorous correction for the confounding features. A subset of 3′UTRs analyzed in (b) were selected and split into four subgroups with respect to the d MTS-RBS . Each subgroup was carefully chosen to have statistically indistinguishable confounding features between each other (see Supplementary Fig. 1f for full versions). The mean values of confounding features, d MTS-RBS , and log 2 (mRNA fold change) are displayed ( *** P < 0.001). Otherwise as in (b). e Association analysis between d MTS-RBS and MT efficacy after deleting DROSHA or DICER. Five miRNAs whose targets show the strongest derepression in response to miRNA removal were chosen and the association between d MTS-RBS and MT efficacy was evaluated. See Supplementary Fig. 1h for full versions. Otherwise as in (d).
identified as direct interactors with miRISC 38,39 and have been revealed to regulate MT for several miRNA targets 9 . However, when partitioning RBPs into direct interactors with AGO and the others, both exhibited a significant correlation between d MTS-RBS and MT efficacy (Fig. 4a), indicating that both groups of RBPs may function as MT enhancers. Therefore, a more plausible mechanism to explain this general impact of RBPs on MT would be that RBP binding alters the local secondary structure of an MTS or its vicinity such that AGO can more readily access the MTS, in a similar manner that Pumilio and PCBP2 regulate the miRNA targets 6,7 .
To examine whether RBP binding leads to the opening of the local secondary structures, we obtained a dataset of DMS-seq that detects unpaired adenines and cytosines at a nucleotide resolution enabling us to accurately probe in vivo secondary structures of endogenous RNAs 40 . When selecting three groups of 3′UTR fragments that have none, lenient, and stringent RBS signals, the 3′UTR fragments with lenient and stringent RBS signals had substantially unpaired secondary structures than those with none and lenient RBS signals, respectively (Fig. 4b). It is noteworthy that when selecting these three groups we have carefully chosen 3′ UTR fragments that have statistically indistinguishable RNA secondary structure in vitro and mRNA expression levels among these groups. Thus, the elevated level of DMS score can be attributed to higher RBP-binding activities instead of RBPs already bound to the open structures.
This association between the RBP binding and elevated DMSseq was more pronounced when a larger number of RBPs bind to the 3′UTR (Fig. 4c), consistent with our previous observation that the overall number of bound RBPs determines MT efficacy ( Fig. 3a). When looking into individual RBPs, 99% of RBPs exhibited significant association, indicating that RBP binding generally results in an opening of the local secondary structure of the 3′UTR in vivo (Fig. 4d), and this could be a common mechanism by which RBPs enhance MT efficacy. Additionally supporting this hypothesis, analysis of AGO2 PAR-CLIP-seq in response to ectopically introduced miRNA 18 exhibited significant increased AGO2 occupancy for those MTSs that have enriched RBP binding (Fig. 5a).
Another possible hypothesis that could explain the correlation between RBP binding and MT efficacy is AGO-mediated recruitment of RBPs where the MTSs of the transfected miRNAs promote the recruitment of RBPs. However, we examined the MTSs of the ectopically introduced miRNAs, which were apparently absent in the WT cells where eCLIP-seq was performed and thus unable to affect eCLIP-seq results. Therefore, the enrichment of RBP binding near the effective MTSs for the transfected miRNAs cannot be explained by the AGO-mediated recruitment of RBPs.
Therefore, we propose a model that takes RBP binding into consideration when explaining MT (Fig. 5b). Compared to the conventional model of ternary interaction among AGO, miRNA, and target mRNA, our proposed model better explains MT efficacy for multiple RBPs (Supplementary Fig. 5a). These results demonstrate that one of the main mechanisms by which RBPs influence MT may be that RBPs open up local secondary structures close to the MTS such that AGO can more easily access the MTS thus improving MT. However, our proposed mechanism does not rule out the previously reported mechanism where protein-protein interaction between AGO and some RBPs  ZC3H11A  SUB1  TIA1  DDX55  DDX6  DGCR8  KHSRP  EIF3D  TIAL1  FASTKD2  PABPN1  AKAP1  CDC40  UCHL5  DROSHA  GRWD1  NKRF  DDX59  TROVE2  BCLAF1  WDR43  LIN28B  ZNF800  DDX52  DHX30  HNRNPUL1  GTF2F1  G3BP1  RBM5  SF3A3  AGGF1  SND1  U2AF1  HLTF  EIF3H  FXR2  DDX3X  RBM15  FAM120A  STAU2  IGF2BP3  PCBP2  TAF15  TRA2A  GRSF1  XPO5  HNRNPK  BUD13  RBFOX2  PRPF4  ILF3  XRCC6  SDAD1  CSTF2  FUS  IGF2BP1  SSB  PPIG  NCBP2  POLR2G  RPS3  SFPQ  U2AF2  UPF1  HNRNPC  PRPF8  SRSF9  SRSF7  HNRNPA1  XRN2  SUPV3L1  BCCIP  SRSF1  SLTM  FKBP4  TBRG4  SMNDC1  EFTUD2  DKC1  PCBP1  QKI  PTBP1  UTP18  NOLC1  HNRNPM  RBM22  FTO  SUGP2  LARP7  HNRNPU  YBX3  EXOSC5  AQR  SAFB  NIP7  SF3B4  MATR3  HNRNPL   0   20 40 q < 0.01 (86 / 100 RBPs, 86%)  For each of the 100 RBPs whose RBSs in HepG2 cell line were identified in the ENCODE eCLIP-seq dataset, the association between d MTS-RBS and MT efficacy was tested after correcting for the four known confounding features for MT by the multiple linear regression (two-sided t-test). The y-axis represents the statistical significance of the observed association: upward and downward directions indicate the MT enhancement (−log 10 (q value)) and the MT suppression (log 10 (q value)) by the RBP, respectively. The number of RBPs that exhibits a significant association between d MTS-RBS and MT efficacy after multiple test correction by the false discovery rate is shown on top. b, c Association analysis between d MTS-RBS and MT efficacy by using 59 microarrays that measured the wholetranscriptome response to ectopically introduced miRNAs or siRNAs into HeLa cell line (b), and for 75 microarrays that measured the whole-transcriptome response to ectopically introduced miRNAs into various human cancer cell lines (c). Otherwise as in (a).
mediates MT, and perhaps both mechanisms function together in a cooperative or independent manner depending on their cellular context.
RBP binding generally enhances miRNA targeting. To examine whether our proposed model is general enough to explain a wide spectrum of various RBPs, we partitioned RBPs into several subgroups with respect to mRNA stabilization function, RNA helicase activity, strand specificity, MTSs in 3′UTR or ORF, or by whether the RBP directly interacts with AGO (Figs. 4a, 6a-g, Supplementary Fig. 5b, c, and Supplementary Discussion). Accordingly, we observed a consistently significant association between d MTS-RBS and MT efficacy for all of subgroups. For instance, d MTS-RBS in ORF also exhibits a significant impact, albeit a modest degree ( Supplementary Fig. 5c). Double-stranded RBPs and nuclear RBPs also seem to have enhancing effect on MT efficacy (Fig. 6e, f). To eliminate the concern of our results confounded by binding sites that overlap with those of single-stranded RBPs or cytoplasmic RBPs, we only examined the RBPs with minimal overlapping binding sites based on the currently annotated data. Even when focusing on the subset of RBPs, the consistent results were observed (Supplementary Fig. 6 and Supplementary Discussion). Taken together, these results support that our findings are robustly general to be expanded to various RBPs and even to other regions such as ORF.
Next, we examined cases where an RBS overlaps with an MTS to inspect whether the competition between miRNA-loaded AGO and other RBPs hampers MT. Our previous results indicate that MT efficacy in such overlapping cases is greater than or equal to that of nonoverlapping cases (see white bar graphs in Fig. 1b-e), suggesting that miRNA-loaded AGO may easily outcompete RBPs. To more definitively investigate the potential competition, we have looked into various subset of MTSs separated by their site types, seed-pairing stabilities, and RBP functions including such cases where MT efficacy is expected to be very weak so that the potential competition between AGO and other RBPs gets more detectable. However, in all cases, MT efficacy was strongest for the MTS that overlaps with RBSs ( Fig. 6a-g and Supplementary Fig. 5b), illustrating that miRNA-loaded AGO can outcompete RBPs instead of RBPs competing against AGO (Fig. 6h). Although this observed lack of competition between RBPs and miRNA-loaded AGO may seem counterintuitive at first, our observation is consistent with previous biochemical studies where the binding affinity for a single RBP and its target RNA is on average >1,000-fold weaker (dissociation constant K D in nanomolar range 41  RBP binding enhances miRNA targeting by improving targetsite accessibility. To validate our proposed mechanism, we performed in vitro gel mobility-shift assays by using disrupted RBSs in three different 3′UTRs and recombinant RBPs (His-FUBP3 and His-PCBP2) (Fig. 7a, b and Supplementary Fig. 7a) and confirmed that the disruption of RBSs reduce RBP binding to mRNA targets (Fig. 7c). Successively, we performed gel mobility-shift assays with recombinant human AGO2 protein and RBPs to confirm whether the disrupted RBSs also reduce AGO binding to the MTS. As a result, miRNA-loaded AGO2 bound to its mRNA target to a much weaker extent when the nearby RBS is disrupted (Fig. 7a, d), suggesting that RBP binding can improve the accessibility of MTSs to AGO2 in vitro. Since FUBP3 has been reported to directly interact with AGO 38,44 while PCBP2 has not been 45,46 (Supplementary Fig. 7b), our results suggest that RBPs generally enhance AGO binding to its MTSs regardless of whether they directly interact with AGO. We also confirmed that the accessibility of AGO is enhanced only when it is loaded with a targeting miRNA ( Supplementary Fig. 7c) by the additional gel mobility-shift assay using AGO2 loaded with non-targeting miRNA (miR-1).
To monitor the transcriptome-wide structural change of RBSs upon RBP binding, we deleted IGF2BP1 in HEK293T cell line ( Supplementary Fig. 7d) and performed DMS-seq in parental and IGF2BP1 KO cell lines. In parental cells, the stringent RBSs of IGF2BP1 exhibited greater DMS reactivity than that of lenient RBSs, demonstrating that the structure of these RBSs opens up upon binding of IGF2BP1 (Fig. 7e). Conversely, this pattern was reversed in IGF2BP1-depleted cells and DMS reactivity decreased   TRA2A  PUM1  METAP2  NIPBL  HNRNPK  U2AF1  SF3B4  TBRG4  CSTF2T  UCHL5  AQR  TROVE2  TIA1  FUS  ZNF622  NSUN2  DROSHA  WDR3  GRWD1  AARS  DDX42  AKAP8L  EFTUD2  UPF1  KHDRBS1  U2AF2  XRCC6  FTO  RBM15  ILF3  FXR1  BUD13  MTPAP  QKI  SRSF7  HNRNPC  LIN28B  GPKOW  SF3B1  SND1  AGGF1  HNRNPUL1  IGF2BP1  FAM120A  PPIL4  IGF2BP2  TAF15  SRSF1  YBX3  SAFB2  MATR3  HLTF  HNRNPL  HNRNPU  as the RBS signal stringency increased (Fig. 7e). The analysis indicates that when IGF2BP1 is absent, the secondary structure of the RBSs is in more closed state and this structural change is highly specific to the IGF2BP1-binding sites. We also tested whether RBP binding improves AGO binding to an MTS in vivo by AGO2-IP followed by RT-qPCR (Fig. 7f). Accordingly, the mRNA level bound to AGO2 in IGF2BP1 knockout cells was significantly lower than that in the parental cells (Fig. 7g). When deleting another RBP, PCBP2 ( Supplementary Fig. 7d), a consistent result was observed (Fig. 7g). Collectively, these results demonstrate that RBP binding improves target-site accessibility of MTSs to AGO both in vitro and in vivo.
To more directly evaluate whether RBP binding leads to improved MT in vivo, we performed luciferase reporter assays for ten 3′UTRs, each of which contains an MTS and a nearby RBS (Fig. 8a). When disrupting an MTS, nine of the ten 3′UTRs displayed a significantly reduced MT efficacy: among these nine 3′UTRs, eight exhibited significant reduction of MT in response to the disrupted RBS, showing that RBP binding indeed improves MT efficacy in vivo.
We further assessed the regulatory impact of RBPs on MT by using RBP knockout cells. When PCBP2 is removed, all of six 3′ UTRs, each of which contains an MTS and a nearby PCBP2 RBS, exhibited reduced MT and then the reduced MT efficacy was restored after overexpressing PCBP2 protein (Fig. 8b). Similarly, in response to IGF2BP1 deletion, both of the examined 3′UTRs, each of which contains an MTS and an IGF2BP1 RBS, exhibited reduced MT efficacy and then rescued when overexpressing IGF2BP1 protein (Fig. 8b).
To examine the impact of the RBPs on MT on a transcriptomewide scale, we have performed mRNA-seq experiments after knocking out PCBP2 or IGF2BP1 and measured the wholetranscriptome response to overexpressed miRNAs. Accordingly, upon the deletion of IGF2BP1 or PCBP2, the miRNA targets were de-repressed for MTSs containing RBSs of the deleted RBP in the vicinity, while such derepression was not detected for non-targets or miRNA targets with MTSs located far from RBSs (Fig. 8c). KO of IGF2BP1, which had been previously reported as an MT suppressor, also led to the derepression of its target mRNAs ( Supplementary Fig. 8a, b), indicating that although some RBPs can suppress MT of specific target mRNAs in a certain cellular context, more generally they function as MT enhancers on a global transcriptome-wide scale. Taken together, our extensive validation experiments and analyses provide multiple solid lines of evidence that RBPs function as MT enhancers by improving the target-site accessibility to AGO.
Evolutionary perspective on the regulatory impact of RBPs on miRNA targeting. To gain a global insight into the impact of RBPs on MT, we investigated the locations of MTSs of 108 broadly conserved miRNA families 47 relative to RBSs on 3′UTRs (Fig. 9a). When considering RBPs profiled in HepG2 cell line, 90% of MTSs included ≥1 RBSs of these RBPs within 100 nt flanking regions. If focusing on conserved and highly conserved MTSs, this pattern was more prominent: 92% and 97% of conserved and highly conserved MTSs, respectively, had ≥1 RBSs within 100 nt flanking regions. Given that there are >1,500 human RBPs 15 , our estimation may well be an underestimation and therefore we extrapolated our analysis to 1,500 human RBPs. Accordingly, we estimate that 100% of conserved and nonconserved MTSs have RBSs in their close proximity. When iterating our analysis for 120 RBPs profiled in K562 cell line, similar results were observed, implying that almost all of human MTSs are likely to be influenced by nearby RBSs.
To gain an evolutionary insight into MT and RBPs, for each of 54 human tissues 48 , we examined whether RBSs of those RBPs expressed in a given tissue tend to co-occur or to mutually exclusively occur near MTSs of evolutionarily conserved miRNAs. As a result, significant enrichment of RBSs was observed near the MTSs in 33% of tested tissues, while significant depletion of RBSs was observed near the MTSs in none of tissues (Fig. 9b). This analysis indicates that locations of RBSs in 3′UTRs are evolutionarily selected to locate near MTSs perhaps to help enhance MT, providing an interesting perspective on MT and coregulating RBPs. Based on these results, we propose a new revised model of MT that takes >1,500 co-regulating RBPs into account (Fig. 9c).

Discussion
In our analysis for an RBS that overlaps with an MTS (Figs. 4a, and 6a-g), we suggest that binding affinity of miRNA-loaded AGO is much stronger (>1,000 fold) than that between an RBP and its mRNA target 41,42 and therefore miRNA-loaded AGO The mean values of log 2 (fold change) of AGO2 occupancy are displayed, and the error bars represent 95% confidence intervals. b Proposed model that takes RBP-binding information into account when analyzing miRNA targeting in comparison to ternary interaction model among AGO, miRNA, and target mRNA. In the ternary interaction model, the changes between miRNA-unbound (ΔG 0 ) and miRNA-bound (ΔG miR ) MFEs are compared assuming there is no RBP binding in the 3′UTR. In our proposed model, additional information of RBP binding is incorporated and therefore the changes between miRNA-unbound (ΔG RBP ) and miRNA-bound (ΔG RBP+miR ) MFEs are expected to better reflect in vivo interactions between the RBS and MTS.  50 and RBPs are considered to function together to accumulatively induce structural opening of RNAs 51 . For instance, an RBP with a helicase activity can initiate opening up a closed RNA structure and then other RBPs can be additionally recruited to induce further structural changes 52 (Supplementary Discussion). Consistent with this hypothesis, we observed that as the number of bound RBPs increases, the RNA structure tends to be more open (Fig. 4c). Once the RNA becomes less structured, miRNA-bound AGO readily accesses the MTS by replacing those RBPs already bound (Fig. 5a), eventually leading to productive MT (Fig. 3a).
Although the present study is based on rigorous bioinformatics analyses and experimental validations, it has the following limitations. First, the global analyses employed the wholetranscriptome datasets of microarray and mRNA-seq. This approach can be justified by the fact that mRNA destabilization rather than translational repression is the dominant mode of MT in steady state and therefore our approach is capturing most of relevant effects of MT 1,53,54 . However, to more conclusively address this issue, future efforts can be made to generate and analyze large-scale proteomics and/or ribosomal footprint datasets, attempting to dissect the contribution of RBPs to translational repression of MT compared to mRNA destabilization. Second, our study focuses on the role of RBPs in MT only in steady-state level and thus lacks an approach to investigating the dynamics of MT. As previously reported, the regulatory mode of MT in transient state can be different from that of steady state depending on the biological contexts 53,55 . Therefore, future studies that aim to revisit the role of RBPs in the context of dynamic regulation of MT such as the maternal to zygotic transition may  Fig. 1d were selected with respect to different site types (a), SPSs (b), function of RBPs on mRNA stability (c), helicase activity (d), strand specificity (e), subcellular localization (f), and previously reported function of RBPs on MT regulation (g). The selected 3′UTRs were split into four subgroups with respect to d MTS-RBS . Each subgroup was carefully chosen to have statistically indistinguishable confounding features among four subgroups, and mRNA fold changes were compared (see Supplementary Fig. 5b for full versions). Otherwise as in Fig. 1d. h Schematic illustration of the MT mechanism for the MTS that overlaps with an RBS. Before AGO-miRNA complex binds to the MTS, RBPs can bind to their RBS (top). One possible scenario afterwards is that the RBP binding competes against miRNA-bound AGO preventing it from productive MT (bottom left). However, our results support the other scenario where miRNA-bound AGO predominantly replaces mRNA-bound RBPs leading to productive MT (bottom right).
help draw a more complete picture of MT. Third, we employed a set of representative mRNA isoforms, which are not specific to a particular cell line throughout our analyses. Although use of mRNA annotation specific to the corresponding cell line can provide more accurate information on the existence of the target sites, the overall association observed by our analyses is not largely dependent on cell-line-specific conditions. Our consistent results using the representative isoforms in various biological contexts support our claim of the general regulatory impact of RBP binding on MT efficacy. Since 3′UTR isoforms vary between cells, these alternative isoforms can affect MT efficacy by inclusion or exclusion of MTSs between different cellular contexts 56 . Cell-type-specific MT regulatory mechanism by different 3′UTR isoforms is an interesting hypothesis worth investigating in the future. Fourth, although we reported that dsRBPs and nuclear RBPs also function as MT enhancers (Fig. 6e, f and Supplementary Fig. 6), the result should be still carefully interpreted because there are caveats of incomplete eCLIP-seq datasets and UV crosslinking bias to single-stranded regions 57 . These concerns may be alleviated when more comprehensive eCLIP-seq datasets and unbiased methods to capture RNA-RBP interactions become available, that would allow us to comprehensively study the effects of the RBPs on MT.
In this study, based on massive-scale analyses of binding sites for >100 human RBPs and on extensive validation experiments, we report that most RBPs enhance MT instead of suppressing it on a global scale, by making the local secondary structure of the MTS more readily accessible to AGO. Our study raises a challenging question about the broadly accepted model of MT that consists of a simple ternary interplay between AGO, miRNA, and mRNA target, proposing a largely revised model that takes a much broader context of >1,500 co-regulating RBPs into additional consideration (Fig. 9c). Our study illuminates the previously unappreciated regulatory impact of RBPs on MT, unraveling the complex nature of the gene regulatory network governed by metazoan miRNAs and their co-regulating RBPs. Undoubtedly, the RBP-binding information, if carefully combined with known determinants of MT, will help more accurately identify functional miRNA targets.
Processing of RBP-binding information of eCLIP-seq data. The ENCODE database provides both raw sequencing data and their processed data of the exact genomic location of RBP-binding sites (RBSs) detected from HepG2 and K562 cell lines 16 . Binding information of 103 RBPs in HepG2 cell line and 120 RBPs in K562 cell line, 150 RBPs in total, were analyzed. Since two replicates of binding information are provided for each RBP, the merged binding information of RBPs were used for analyses in corresponding cell line and thus the RBSs that have been detected only in one replicate was also included in our list of RBSs. When analyzing the data of other cell lines rather than HepG2 or K562, RBP-binding information of both cell lines was used. In case there exist data of RBPs profiled in both cell lines, RBSs were merged together.
To examine the impact of individual RBPs on miRNA targeting (MT) efficacy in Fig. 2a-c, following filtering steps were applied for selecting RBPs to be used in the analyses. Among the total profiled RBPs, SERBP1 was discarded since it has a low number of RBSs in the 3′UTRs (<300) while other RBPs have >1,000 RBSs in the 3′UTRs. When analyzing HepG2 and HeLa datasets, those RBPs whose expression levels belong to the lower 50% were additionally excluded. As a result, 100 out of 103 RBPs, 147 out of 150 RBPs, and 149 out of 150 RBPs were used for investigation on the datasets of HepG2, HeLa, and other human cancer cell lines, respectively.
When analyzing the association between the RBP binding and mRNA secondary structures (Fig. 4b-d), those RBSs that have strong statistical significance, which the ENCODE processed data provides, were defined as 'stringent RBS', while the others as 'lenient RBS'.
Multiple linear regression analysis. We used multiple linear regression (MLR) to correct for confounding factors known to influence MT (Figs. 1c, 2a-c, and Supplementary Fig. 5a). In order to correct for the effect of potentially confounding features contributing to MT (local AU content, target-site abundance, seed-pairing stability of corresponding miRNA, and 3′UTR length) [27][28][29][30][31][32][33][34] , MLR models with aforementioned features were fitted to the log 2 (fold change) of the target mRNAs for each dataset of HepG2, HeLa, and other human cancer cell lines. After fitting the MLR model, regression residuals were calculated by subtracting fitted values of the MLR model from the observed log 2 (fold change) for each mRNA. The regression residual is interpreted as a remaining information of MT efficacy after correcting for the effect by four known confounding factors and it was used to test whether the distance between a miRNA target site (MTS) and the nearest RBS on the 3′UTR, denoted as d MTS-RBS , is associated with MT efficacy (Fig. 1c). 'OLS.fit()' regression function in the 'statsmodels' package in Python 58 was utilized to fit the MLR models.
When investigating the association of RBP binding and MT efficacy for individual RBPs, MLR models were constructed for each RBP. From miRNA overexpression dataset of HepG2, HeLa and other human cancer cell lines, log 2 (fold changes) for the mRNAs which contain a single 7, 8mer MTS on their 3′ UTRs were collected, in order to clearly observe the impact of RBPs on a single MTS. For each RBP, the MLR model was fitted for the feature of RBP binding, either d MTS-RBS (Fig. 2a-c) or structural change of mRNA ( Supplementary Fig. 5a), and the previously reported confounding features of MT. The association between MT efficacy and d MTS-RBS was measured by the P value of t-test. To correct the multiple testing problem, the P values were corrected by the false discovery rate (Fig. 2a-c).
Correction for potentially confounding features of MT. The confounding features of 3′UTR length, local AU content, TA, and SPS might also be potentially confounded with the RBP binding so need to be carefully controlled. To rigorously correct these confounding effects, we have devised an in-house algorithm facilitating to sample 3′UTRs into four different groups only according to d MTS-RBS while keeping the confounding effect similar across the groups. 3′UTRs containing a single 7, 8mer MTS and one or more RBSs were defined as the input of sampling process. The sampling process comprises of multiple rounds of selection. For a round of selection, a set of four 3′UTRs was selected with the following criteria: (i) In a same set, the site type of MTSs should be identical to one of the 8mer, 7mer-m8, and 7mer-A1, (ii) The min-max range of confounding effect should fall within the specified cutoffs, that were determined by the Gaussian process implemented in the 'bayes_opt' Python package to maximize the number of selected 3′UTRs, (iii) One 3′UTR should have an overlapping MTS with the nearest RBS (d MTS-RBS = 0) and other three should have MTS separated from the nearest RBS (d MTS-RBS > 0). Fig. 7 Disrupted RBSs or RBP knockout reduce the target-site accessibility to AGO. a 3′UTR structures of RNA substrates (SSB, UBA1, and CRAT) used for gel mobility-shift assays shown in b-d. The 3′UTRs of SSB and UBA1 were also used as targets for luciferase assays in Fig. 8a, b. The miRNA target sites (MTSs) and RBP-binding sites (RBSs) are shown as red and blue boxes, respectively, with the mutated positions indicated by orange boxes. b Histidinetagged recombinant RBPs (His-FUBP3 or His-PCBP2) and GST proteins were used for gel mobility-shift assays shown in (c, d). c Gel mobility-shift assays for RBS WT and RBS MUT with His-FUBP3 or His-PCBP2 with GST protein used as a negative control (top). The free RNA and RNA:RBP complex bands are shown as black and blue rectangles, respectively. Mean fractions of the bound RNA:RBP complexes ±95% confidence intervals are displayed (bottom, n = 3). d Gel mobility-shift assays with (1) 3′UTR, (2) 3′UTR and rhAGO2, (3) 3′UTR and RBP, and (4) 3′UTR, RBP, and rhAGO2. Otherwise as in (c). e DMS reactivities on A and C nucleotides of 3′UTRs were measured by comparing DMS counts of DMS(−) and DMS(+) samples. Corresponding nucleotides were divided into three groups by the RBP-binding signal of IGF2BP1 eCLIP-seq dataset (No RBS, Lenient, and Stringent). DMS reactivities between WT and RBP KO were compared (two-sided Wilcoxon's rank-sum test). The mean DMS reactivity values ±95% confidence intervals are displayed. f Experimental procedure for AGO2-immunoprecipitation (IP) followed by western blot and RT-qPCR. g AGO2-IP followed by western blot and RT-qPCR in AGO2overexpressed HEK293T parental and RBP KO (IGF2BP1 or PCBP2) cells. Protein levels in the input and IPed samples were visualized by the western blotting (left). Relative RNA levels of each input and IPed samples were quantitated in parental and RBP KO cells and normalized to each input sample (right, twosided Wilcoxon's rank-sum test, *P < 0.05, **P < 0.01, ***P < 0.001, n = 9 for PCBP2 and n = 6 for IGF2BP1). The mean relative RNA level of the IPed samples ±95% confidence intervals are displayed. The RNA level of KATNA1 was used as a negative control and U6 snRNA-level served as a technical control of AGO2-IP. P values are provided in Source Data.
3′UTRs which met the criteria above were then assigned to the first to fourth groups in ascending order of the value of d MTS-RBS . Next rounds of selection of four 3′UTRs proceeded until the inputs were exhausted. If any confounding factor exhibited significant difference across the groups, the sampling process was retried with the narrower cutoffs mentioned at the second criterion. After the sampling process finished, comparison of the log 2 (fold change) or the confounding features across the groups was performed using Wilcoxon's rank-sum test of 'SciPy' package in Python. Accordingly, we can clearly divide 3′UTRs solely dependent on d MTS-RBS without any concern that the confounding effect involves across groups. By using the dataset of DICER or DROSHA knockout in HCT116 cell lines 35 , we performed the similar analysis to detect the derepression of target mRNAs. To select miRNAs whose targets show the strongest derepression in response to miRNA removal, a 2 × 2 contingency table was constructed for each miRNA by examining whether its MTS was included in the 3′UTR and whether the 3′UTR was highly de-repressed.  The most significant five miRNAs after χ 2 tests were chosen and used for the analysis.
DMS-seq analysis. To distinguish whether RBP binding leads to the 3′UTR structural change, we used a DMS-seq dataset generated in K562 cell line with treatment of DMS (Accession ID: GSM1297493) 40 and the binding information of 90 RBPs from the ENCODE eCLIP-seq in K562 cell line. We monitored mRNA expression profile in K562 cells 23 and used top 50% of the most highly expressed mRNAs for the analyses because mRNAs with low expression level tend to have depleted signal of DMS-seq, which can lead to an inaccurate detection of unpaired regions for the mRNAs. Across the highly expressed human mRNAs, we collected 3′UTR fragments with all possible 100 nt windows shifted by every 10 nts. The 3′ UTR fragments were separated into three groups: fragments that contain one or more RBSs with high confidence as explained above (stringent RBSs), fragments that contain RBSs with moderate signal only (lenient RBSs), and the other fragments without any RBS (No RBS). 3′UTR fragments were subsampled to have similar values of the predicted minimum free energy and the expression levels of mRNAs where the fragments are originated among the groups.
We normalized DMS read counts to RPM scale for each replicate of DMS-seq data. The DMS levels were then calculated by averaging the normalized counts of multiple replicates. We removed the background signal by comparing DMS levels between in vivo and denatured control. The resulting value, 'DMS score', is defined by Eq. (1). Also for in vitro sample, the DMS score was calculated by Eq. (2). For each selected fragment, ΔDMS score was measured by Eq. (3) where comparing DMS scores between in vivo and in vitro and the scores were compared among the groups (Fig. 4b). To assess the correlation between the number of bound RBPs and opening of the secondary structure of mRNAs, 3′UTR sequences with stringent RBSs were collected. These sequences were separated with respect to the number of bound RBPs, and their ΔDMS scores were normalized by the expression level of 3′ UTR (Fig. 4c). The enrichment of ΔDMS score was tested by comparing the corresponding score in the region without any RBP binding (Fig. 4d) The HEK293T parental and IGF2BP1 KO DMS-seq libraries were prepared using the protocol a previous study 40 with following modifications. The fragmented libraries were ligated to 3′ adapter from a recent study 59 with T4 RNA ligase2 truncated K227Q (NEB), and size-selected on a 10% UREA polyacrylamide gel. Next, the library was reverse transcribed and ligated to 5′ adapter 59 with RNA ligase 1 high concentration (NEB) at 22°C overnight, as described in another previous study 16 . Then, the library was size-selected on a 10% UREA polyacrylamide gel and amplified by PCR. The final size-selected DMS-seq libraries were sequenced using HiSeq 2500 (single-end, 51 bp and 101 bp, Illumina). All the information of the primers and adaptors used in this study are listed in Supplementary Table 4. DMS count was calculated by detecting 5′ end of mapped reads (RT-stops). DMS reactivity was measured on A and C nucleotides of 3′UTRs by comparing DMS counts of DMS-treated and DMS-untreated samples (Fig. 7e).
Gel mobility-shift assay. To assess if RBPs generally enhance AGO binding to MTSs regardless of whether they directly interact with AGO, FUBP3, and PCBP2 were chosen for our gel mobility-shift assay because FUBP3 is reported to interact with AGO2 38 while PCBP2 has not been reported to interact with AGO2 45 . As one of FUBP3 binding targets, SSB 3′UTR, whose RBP-binding signal is indicated in Supplementary Fig. 7a, was chosen for the gel mobility-shift assay (Fig. 7c, d). To confirm a non-direct interactor enhances AGO2 binding to the 3′UTR, we have performed gel mobility-shift assays with UBA1 and CRAT 3′UTRs that contain PCBP2-binding sites (Fig. 7c, d).
PCR fragments for 3′UTRs of SSB, UBA1, and CRAT containing T7 promoter were amplified from psiCHECK-2-SSB, -UBA1, and -CRAT plasmids, respectively. To prepare for RBS MUT sequences that we designed as described in Supplementary  Table 4, we performed site-directed mutagenesis from the plasmids containing RBS WT . The RNA substrates were synthesized with [ 32 P]-UTP by in vitro transcription and purified by resolving on a denature-PAGE. To prepare N-terminal Histidine-tagged protein, coding sequences of each FUBP3 and PCBP2 were amplified by specific primers (Supplementary Table 4) and inserted to pET-28b vector. Both recombinant histidine-tagged FUBP3 and PCBP2 proteins were purified by the HisTrap column (GE healthcare), according to manufacturer's instructions. The GST protein was purified by the GST affinity beads (Elpis Biotech) to be used as a control protein. The purified recombinant proteins were resolved on 12% SDS-PAGE and stained with Coomassie brilliant blue (BIORAD) (Fig. 7b).
RNA immunoprecipitation (IP) and qPCR. FLAG-tagged human AGO2 was cloned to pcDNA3.1 vector. HEK293T parental and RBP KO (PCBP2 and IGF2BP1) cells were seeded in a six-well plate at~60% confluency. After 24 h of incubation, the FLAG-hAGO2 plasmids were transfected into cells using Lipofectamine 2000 according to the manufacturer's protocol. After 24 h incubation, the cells were harvested, and then lysed by lysis buffer (20 mM HEPES-KOH at pH 7.4, 150 mM KOAc, 1.5 mM Mg(OAc) 2 , 0.1% Triton X-100) with 1× Protease inhibitor cocktail (Roche). The precleared lysate was incubated with anti-FLAG M2 affinity gel (Sigma) at 4°C overnight. Immunoprecipitated samples were washed four times with wash buffer containing 300 mM KOAc. Total RNAs were extracted by TRIzol (GeneAll) to be subjected to RT-qPCR. . The normalized fold repressions were compared between MTS WT RBS WT and MTS WT RBS MUT (two-sided Wilcoxon's rank-sum test, *P < 0.05, **P < 0.01, ***P < 0.001, n = 12, See Methods for normalization procedures). The changes in the minimum free energy of the 3′UTRs (ΔG miR − ΔG 0 and ΔG RBP+miR − ΔG RBP ) are listed at the bottom. KATNA1 was used as a technical control. The mean values ±95% confidence intervals are displayed. P values are provided in Source Data. b The design of constructs and the expected changes in the secondary structures of mRNAs are shown (top left): parental, RBP KO, and rescue conditions of the deleted RBP are depicted by blue, orange, and green colors, respectively, with striped colors representing the mutated MTS. The protein abundance was quantified by western blot (top right). For the rescue experiment of each RBP KO, FLAG-tagged RBP was quantified (top right). GAPDH levels serve as a loading control. The MTS WT values were normalized to MTS MUT and were compared between parental and KO and between KO and rescue (two-sided Wilcoxon's rank-sum test, *P < 0.05, **P < 0.01, ***P < 0.001, n = 12). Otherwise as in (a). c Transcriptome-wide response of miRNA targets after RBP removal. mRNA expression fold changes after miRNA overexpression were measured for IGF2BP1 or PCBP2 KO cells and compared to HEK293T parental cells, with the x-axis indicating the difference of the log 2 (mRNA fold change) values. Target 3′UTRs with a single 7, 8mer MTS were selected if the distance between an MTS and the nearest RBS, denoted as d MTS-RBS , is short (<100, orange) or long (≥100, blue). After controlling for confounding features of MT (right), distributions of log 2 (mRNA fold change) were compared between the subgroups (left, two-sided Kolmogorov-Smirnov test). Distribution of miRNA non-targets ('No-site', gray) was plotted for comparison. The mean values of confounding features are displayed and the error bars represent 95% confidence intervals.
RNAs were isolated, treated with DNase I (Takara), and then converted to cDNA by reverse transcription using oligo dT primers and PrimeScript reverse transcriptase (Takara). Quantitative PCR was performed using SYBR green I master mix (Roche) with the gene-specific primers (Supplementary Table 4). As a technical control of the IP experiment, U6 snRNA was detected using Taqman miRNA assay (Applied Biosystem). The qPCR was performed with LightCycler II 480 v.1.5.1 (Roche) following the manufacturer's protocol. The Ct values of IP samples were normalized to those of each input sample using the following equations.