Stability of miRNA 5′terminal and seed regions is correlated with experimentally observed miRNA-mediated silencing efficacy

MicroRNAs (miRNAs) are key regulators of sequence-specific gene silencing. However, crucial factors that determine the efficacy of miRNA-mediated target gene silencing are poorly understood. Here we mathematized base-pairing stability and showed that miRNAs with an unstable 5′ terminal duplex and stable seed-target duplex exhibit strong silencing activity. The results are consistent with the previous findings that an RNA strand with unstable 5′ terminal in miRNA duplex easily loads onto the RNA-induced silencing complex (RISC), and miRNA recognizes target mRNAs with seed-complementary sequences to direct posttranscriptional repression. Our results suggested that both the unwinding and target recognition processes of miRNAs could be proficiently controlled by the thermodynamics of base-pairing in protein-free condition. Interestingly, such thermodynamic parameters might be evolutionarily well adapted to the body temperatures of various species.

experiments 2,3,27-33 suggested that siRNA could generate off-target effects through a mechanism similar to that of target silencing by miRNA 2,3,32 . Thus, miRNA-mediated silencing and siRNA-based offtarget effects may use similar machinery as downstream target recognition processes.
Previously, we found that the efficacies of seed-dependent offtarget effects of siRNAs are strongly correlated with the calculated thermodynamic stabilities of seed-target duplexes 33 . However, unlike siRNA off-target effects, the efficacy of miRNA-mediated gene silencing was not simply correlated with seed-target duplex stability. Here, we demonstrated using mathematization of base-pairing stability that the efficacy of miRNA-mediated gene silencing was determined principally by the combinatorial thermodynamic parameters, which might reflect the easiness in unwinding in addition to the base-pairing stability in the seed-target duplex. Furthermore, because temperature is a key regulator of base-pairing stability, the thermodynamic properties of miRNAs of various species with different body/rearing temperatures were evaluated. Interestingly, we found that the thermodynamic stability between the miRNA seed region and target mRNA is well correlated to the body temperatures of various species.

Results
Variation in miRNA-mediated gene silencing activity. To determine the relationship between miRNA structures/sequences and their direct silencing efficacies, we performed reporter assays. Three tandem repeats of partially complementary sequences containing a seed-matched sequence of each of 20 arbitrarily chosen human miRNAs ( Figure S1) were introduced into the Renilla luciferase 39 UTR in the psiCHECK plasmid, hereafter called psiCHECK-SM (Figure 1a). The pGL3-Control, encoding firefly luciferase, and each psiCHECK-SM construct were transfected into human HeLa cells along with the corresponding miRNA. Twentyfour hours after transfection, the relative luciferase activity was measured as a function of miRNA concentration (Figure 1c). Little or no silencing effects were observed by transfection with any of the 20 miRNAs at 0.05 nM, and six miRNAs at 0.5 nM reduced the luciferase activity to below 50%. For six miRNAs (miR-373*, miR-548d-5p, miR-606, miR-335, miR-643, and miR-199b-3p), no appreciable silencing effects were seen even when the miRNA concentration was increased to 5 nM. We also performed the reporter assay mimicking the RNAi effect using psiCHECK-CM, which has a complete-matched sequence of each miRNA in the Renilla luciferase 39 UTR (Figure 1b). Previously, we reported that siRNAs can be divided into three groups; namely, class I: highly functional, class II: intermediate, and class III: ineffective siRNAs, or siRNAs with long GC stretches (.9) have little silencing activity 34 . Of the 20 miRNAs, 15 miRNAs were classified into class I and five into class II. Excluding miR-296-5p, all of the miRNAs strongly reduced the relative luciferase activities in the psiCHECK-CM reporter assay to less than 50% at 0.5 nM (Figure 1c). The same RISC may play a role in RNAi and miRNA-mediated gene silencing 35 . Thus, part of the miRNAs may have little or no gene silencing activity on the partially complementary target mRNAs, even when engaged by the RISC. In contrast, miR-296-5p reduced the luciferase activity of psiCHECK-SM, but not psiCHECK-CM. miR-296-5p is classified into class II and has a 10-nt GC stretch. Thus, miR-296-5p can be held on the RISC for miRNA-mediated gene silencing, but could not repress the complete-matched target, likely due to the long GC-stretch by which the cleaved target might not be easily released from the RISC. In a previous study, we reported that the content of luciferase mRNA produced within cells was about 300 copies/ng total RNA (one-hundredth that of b-actin mRNA) under our experimental conditions, and that the luciferase activities measured using different psiCHECK-SM constructs were almost proportional to the levels of mRNA 33 . Thus, under our conditions, the majority of the luciferase activity reduction was attributable to miRNAmediated luciferase mRNA degradation.
Correlation between miRNA silencing and combined stability. We previously demonstrated the seed-dependent off-target effects of siRNA measured using a luciferase reporter assay at a concentration of 50 nM, which was negatively correlated with thermodynamic stability in the duplex formed between the seed region of the siRNA guide strand and its target mRNA with a correlation coefficient (r) of 20.72 33 . Hence, because the siRNA seed region is a primary target-recognition region, it is possible that the highly stable seed-target duplex results in a strong siRNA seed-dependent off-target effect. This correlation was successfully calculated when melting temperature (Tm) was used as a measure of duplex thermodynamic stability. The off-target effect of siRNA is generated through a mechanism similar to that of target silencing by miRNA 2,3,32 . Unlike the siRNA, the efficacy of miRNA-mediated luciferase gene silencing using psiCHECK-SM was poorly correlated (r 5 20.50) with the calculated Tm values in the seed (positions 2-8)-target duplex (Tm 2-8 ) at 100 mM NaCl (Figure 2a).
We then considered that miRNA-specific features that are involved in the silencing process before target recognition may be responsible for the efficacy of miRNA-mediated silencing. miRNA has specific structural features such as an internal bulge or mismatch (see Figures 7 and S1), but siRNA is simply composed of perfectly complementary double-stranded RNAs. We examined the involvement of thermodynamic stability in the miRNA duplex (miTm) by calculating Tm values considering the internal bulge/mismatch. To determine the optimal region with a high correlation with silencing efficacy, miTm values in the most-to-least optimal regions in the miRNA duplex were calculated. Each of these values was subtracted from the Tm 2-8 value and their correlation coefficient was estimated with respect to silencing efficacy (luciferase activity) according to the formula Tm 2-8 2 k 3 miTm x-y , where x is the start nucleotide position (1 # x # 17), y is the end nucleotide position (2 # x # 18), and k is a multiplicative factor ( Figure 3). Surprisingly, the stability in the 59 terminal region had a significant effect on silencing efficacy. When the optimal k value was used for each region, the Tm 2-8 2 k 3 miTm x-y values in region A shown in Figure 3a starting from position x 5 1,5 and ending at position y 5 2,9 were closely correlated with silencing efficacy (r 5 20.51 to 20.77). The strongest negative correlation coefficient (20.77) was obtained when miTm 1-5 was used; the resultant formula was Tm 2-8 2 k 3 miTm [1][2][3][4][5] (Figures 2d and 3a). Weak but significant correlations were also observed in region B from nucleotides 13,16 to 16,17 (r 5 20.53 to 20.57). The regions from nucleotide 1,12 ending at 10,18 showed little or no correlation. The optimal k values were independently calculated for each region, and the results at positions x 5 1 to y 5 2,8 are shown in Figure 3b. The prominent peaks of correlation coefficients were obtained around 0.5. The closest relationship was found when k was 0.53. miTm 1-5 alone showed little correlation (r 5 0.50) with silencing efficacy (Figure 2b), as in the case with Tm 2-8 alone (r 5 20.50) (Figure 2a). There was no correlation between the values of miTm 1-5 and those of Tm 2-8 (r 5 0.14), although these regions partially overlapped, indicating that the 59 terminal structures of miRNAs are extraordinarily diversified independent of their nucleotide sequence ( Figure 2c). Thus, the silencing efficacy might be estimated based on the combinatorial parameters representing the stability of the miRNA 59 terminal duplex, miTm [1][2][3][4][5] , and base-pairing stability between the seed region and target mRNA, miTm [2][3][4][5][6][7][8] . The following formula should appropriately predict miRNA-mediated gene-silencing efficacy: Tm 2-8 2 0.53 3 miTm 1-5 .
A considerable deviation was also observed in luciferase activity measurements (Figure 2d). This may have been due in part to differences in the non-seed sequence and/or its counterpart in the  target mRNA (see Figure S1, right column) because the target sequences that correspond to the non-seed region make an appreciable contribution to target recognition by miRNAs and/or siRNAs in microarray profiling 2,3,28 .
Different silencing efficacies of miRNAs with common seed sequences or common nucleotide compositions. Although miRNAs recognize their target genes based on the complementarity to the seed region, our results suggested that miRNAs with the same seed sequences but different duplex structures may have different silencing efficacies according to the Tm 2-8 2 0.53 3 miTm 1-5 values. To evaluate this possibility, we investigated the silencing activities of miR-302a/372/373/520c-family miRNAs containing a common seed sequence (AAGUGCU), so their Tm 2-8 values are identical. Members of this miRNA family are known to induce miRNAinduced pluripotent stem (mirPS) cells 36 , and we have reported that the expression levels of many genes with seed complementary sequences in their 39UTRs are commonly regulated 20 . As shown in Figure 4a, miR-372 and miR-520c-3p have the same Tm 2-8 2 0.53 3miTm 1-5 values of 40.1uC, but the values of miR-373 and miR-302a are 35.8uC and 39.8uC, respectively. The reduction of luciferase activity of psiCHECK-SM was weak (38% of relative luc activity at 0.5 nM miRNA duplex) after treatment with miR-373, while rather strong activity (18%) was induced by the miR-302a duplex. Furthermore, miR-372 and miR-520c duplexes significantly reduced the luciferase activities of psiCHECK-SM to 8% and 6%, respectively. The results indicated that their silencing efficacies were correlated with their Tm 2-8 2 0.53 3 miTm 1-5 values. We also tried to evaluate the efficiencies of the other miRNAs including miR-1302-1/1302-2/ 1302-7/1302-8, miR-7-1/7-2/7-3, and artificially mutated miR-30c-1s and miR-643s, which have identical seed sequences but different duplex structures. However, we could not perform the experiments since part of these miRNAs were not successfully annealed, leading to failed miRNA duplex formation, under our annealing conditions (see Methods).
Thermodynamic properties of miRNAs in various species. Small RNA-mediated gene silencing is a conserved phenomenon in metazoans 37 . In this study, we found that the efficacy of miRNA-mediated silencing could be successfully determined based on the thermodynamic properties of protein-free RNA duplexes. Thermodynamic  propensity is naturally controlled by ambient temperature. Because the systemic or rearing temperatures of each species differs, the actual stability of the RNA duplex should differ by species. Thus, it is possible that functional miRNAs in each species are evolutionarily adapted according to temperature. We analyzed the thermodynamic parameters Tm 2-8 , miTm [1][2][3][4][5] , and Tm 2-8 2 0.53 3miTm 1-5 using miRNAs of 16 different species registered in miRBase (Table S1). The average Tm 2-8 values varied widely, from 29.9 to 37.6uC (Figure 5a). The planarian Schmidtea mditerranea and ascidian Ciona intestinalis are heterothermic animals that are usually maintained at approximately 10-15uC. The African frog Xenopus tropicalis, nematode Caenorhabditis elegans, silkworm Bombyx mori, Drosophila melanogaster, Drosophila pseudoobscure, lamprey Petromyzon marinus, and zebrafish Danio rerio are reared at 23-27uC. The body temperature of homothermic animals such as human Homo sapiens, mouse Mus musculus, dog Canis familiaris, horse Equus caballus, orangutan Pongo pygmaeus, and pig Sus scrofa is about 37uC. The body temperature of the chicken Gallus gallus is highest at about 42uC. In the present study, low-temperature animals preserved miRNAs with significantly low average Tm 2-8 values (Table S1 and Figure 5b). In contrast, the miRNAs of high-temperature animals had high average Tm 2-8 values. The cumulative fraction of the Tm 2-8 values of the miRNAs of S. mediterranea, C. intestinalis, D. melanogaster, H. sapiens, and G. gallus clearly varied according to their temperatures (Figure 5c). An apparent correlation (r 5 0.83) between Tm 2-8 values and temperature was observed (Figure 5b), suggesting that miRNA sequences in the seed region are evolutionarily adapted to temperature. 3 miTm x-y , which shows the highest correlation coefficient with the efficacy of miRNA-mediated gene silencing (relative luciferase activity), using the most-to-least optimal regions. The optimal multiplicative factor k, which was determined as shown in Figure 3b, was used to calculate Tm 2-8 2 k 3 miTm x-y for each region. Red lines indicate the region represented in Figure 3b  We further calculated the average miTm 1-5 value for the miRNAs of each species. The values varied significantly from 225.2 to 29.2uC (Table S1 and Figure 5d). The average miTm 1-5 values of B. mori (224.3uC) and C. intestinalis (225.2uC) were significantly lower than those of the others (Table S1 and Figure 5e). The cumulative fractions of miTm 1-5 also varied (Figure 5f), but were not correlated with temperature (r 5 0.33; Figure 5e), suggesting that evolutionary pressure on the miRNA 59 terminal regions is weak.
We calculated the predicted miRNA-mediated silencing efficacy using the formula Tm 2-8 2 0.53 3 miTm 1-5 for 1,902 human miRNAs registered in miRBase ( Figure 6 and Table S2). The values ranged from 4.6uC to 203.8uC, suggesting that human miRNAs have enormously divergent silencing efficacies.

Discussion
The efficacy of miRNA-mediated gene silencing was estimated based on the combinatorial thermodynamic parameters, Tm 2-8 2 0.53 3 miTm 1-5 , of protein-free RNA duplexes in the regions administrating unwinding efficacy (miTm 1-5 ) and base-pairing stability with target mRNA (Tm 2-8 ), as shown in Figure 7. Our results are in excellent agreement with the known silencing machineries. First, an RNA strand containing the thermodynamically less stable 59 end is preferentially entrapped on the RISC 33,38,39 . The internal bulge/mismatch is thought to form a less stable base-pairing, suggesting that miRNAs with such relaxed structures on the 59 end are easily retained by RISC. Second, miRNA recognizes target mRNA with seed-complementary sequences 2,3,32 . The stability between seed region and target mRNA is a determinant of the efficacy of siRNA off-target effects 33 . Thus, the high stability of the seed-target duplex might function as a positive factor, but that in the miRNA duplex the 59 end might be a negative regulator of target gene silencing. Although these regions overlap, the stability of the duplex formed between the miRNA seed and target mRNA is defined by the nucleotide sequence, while the stability within the miRNA duplex is largely attributable to structural features rather than the nucleotide sequence. Thus, our results suggest that these two regions may   coordinately but independently regulate the silencing machinery. Furthermore, the coefficient factor, 0.53, in the mathematical formula suggested that the contribution of the nucleotides 1-5 in miRNA for the silencing efficacy might be about half of that of nucleotides 2-8 in the seed region.
Seed pairing is known to be both necessary and sufficient for target regulation by microRNAs in some experimental contexts 40 . We performed experiments focused on the canonical base pairing in the seed region, and presented a model of predicting miRNA silencing efficacy. However, alternative modes of target recognition by miRNAs have been reported recently, including 39-compensatory sites 3 , centered sites 41 , or bulged sites 42,43 . In our experiments, it was also apparent that region B, corresponding to the 39-compensatory sites from nucleotides 13,16 to 16,17, contributed to the silencing efficacy to some extent (Figure 3a). The nucleotides from 13,16 are required for increased efficacy but are only slightly effective compared to those without the supplementary pairing 3 ; they play a modest role in target recognition. Furthermore, this site can compensate for a single-nucleotide bulge/mismatch in the seed region. We attempted to incorporate the miTm [13][14][15][16] values into the thermodynamic parameters, but no significant improvement was observed (data not shown), probably because only target mRNAs with seed-complementary sequences without internal bulges were used in our luciferase reporter assays. It was also reported that miRNA represses the expression of mRNAs with seed-complementary sequences with bulges in the mRNA side 43 . Although we did not examine the effects on the bulged targets in this study, the Tm 2-8 value in seed-target duplex is identical to that in the duplex formed between miRNA and a non-bulged target, suggesting that our model is applicable. However, Ha et al. 42 reported that miRNA also represses mRNAs when bulges are formed in the miRNA side of the seed-target duplex. In most of these cases, the Tm 2-8 values should change according to the bulged structures, possibly leading to different silencing efficacies. Furthermore, miRNAs that lack both perfect seed pairing and 39-compensatory pairing and instead have 11-12 contiguous Watson-Crick pairs in the center of the miRNA are also functional as miRNA cleavage substrates in vitro 41 . In addition, a given miRNA was shown to generate non-canonical functioning heteroduplexes with targets that do not contain the miRNA seed by molecular dynamics analyses, indicating that the spectrum of potential targets for a miRNA includes a wide-spectrum of seed-less targets and thus substantially differs from what is anticipated based on the canonical seed mode 44 . However, we did not consider centered sites or noncanonical seed-less targets in this study.
We previously reported that the efficiency of siRNA seed-dependent off-target silencing is strongly correlated with Tm 2-8 values 33 . Since siRNA forms perfectly complementary double-stranded RNA, the Tm values at positions 1 to 5 (Tm 1-5 ) showed a strong positive correlation with the reduction in luciferase activity (r 5 0.61; Figure S2b) similar to Tm 2-8 (r 5 0.80), as might be expected ( Figure  S2c). Thus, Tm 1-5 might not be a good parameter for estimating the strand separation (unwinding) efficiency of siRNA. The efficacy of the siRNA off-target effect was not correlated with Tm 2-8 2 0.53 3  target pathway and miRNA-mediated silencing pathway, respectively. The miRNA/siRNA duplex is unwound to a single-stranded RNA, which loads onto the RISC, and acts as a guide to recruit the silencing complex to the 39 UTR of the target mRNA to promote translational repression or cleavage. In the siRNA pathway, the efficiencies of the siRNA seed-dependent off-target effects were strongly correlated with the calculated thermodynamic stabilities in the seed-target duplexes, Tm 2-8 . In the miRNA pathway, the efficacies of miRNA-mediated silencing are determined by the combined thermodynamic parameters that might reflect their unwinding properties (miTm [1][2][3][4][5] ) in addition to their base-pairing stabilities in the seed-target duplex (Tm 2-8 ).
www.nature.com/scientificreports SCIENTIFIC REPORTS | 2 : 996 | DOI: 10.1038/srep00996 Tm 1-5 (r 5 0.47; Figure S2d). A possible explanation is that the region responsible for unwinding in the siRNA duplex differs from region 1-5 in miRNA. To determine the optimal region for the siRNA off-target effect, Tm values in the most-to-least optimal regions in siRNA duplexes were calculated in the same manner as for miRNA ( Figure S3a). However, non-significant improvement in the efficiency of the siRNA off-target effect was observed when the optimal values (x 5 5, y 5 14, and k 5 1.0) were used in the formula Tm 2-8 2 k 3 Tm x-y . Furthermore, no prominent peak in k value was detected for any position (Figure S3b), suggesting that the efficiency of the siRNA off-target effect is determined primarily by the stability of the seed-target duplex, as reported previously 33 . For siRNA-based off-target effects, the siRNA guide strand in the RISC cleaves the passenger strand with the Ago2 protein leading to its dissociation from RISC 11,12,45 . However, because most miRNA duplexes contain bulges to prevent cleavage, the miRNA* strand dissociates by unwinding 14,15,46 . Thus, instability in the 59 terminal region might be essential for unwinding of the miRNA duplex, but might not be necessary for cleavage of the siRNA passenger strand.
The thermodynamic properties of miRNAs vary significantly among organisms ( Figure 5). One causal factor of this may be natural selection in response to differences in temperature (body or rearing temperature). Considering the chemical characterization of nucleic acids, an RNA duplex with low thermodynamic stability would not be formed in the cells of an organism with a high temperature. The Tm 2-8 and Tm 2-8 2 0.53 3 miTm 1-5 values of the miRNAs of 16 different species were strongly correlated with the temperature of each species (Figure 5c and 5i), indicating that organisms with higher and lower temperatures possess miRNAs with higher and lower seed-target duplex stabilities, respectively. This suggests that the stabilities in the miRNA seed sequences are evolutionarily selected according to the adaptive temperature of each organism. In contrast, the diversity of miTm 1-5 values was marginally affected by temperature (Figure 5f). However, since the values varied significantly, positions 1-5 may produce variation of silencing efficiencies independent of temperature. Thus, the functions of miRNA 59 terminal and seed regions in miRNA-mediated gene silencing may differ. Some miRNAs, such as let-7 47,48 , miR-34 49 , miR-124 50 , and miR-125 49 , are known to be evolutionarily conserved and their target sites are conserved among various species 51 . Furthermore, the conserved miRNAs are ancient animal miRNAs whose localization in tissues are closely coupled in evolution 52 . The conserved miRNAs may be functional across species if their Tm 2-8 values are high. The Tm 2-8 values of let-7, miR-34, miR-124, and miR-125 were high at 37.6, 45.5, 40.8, and 45.2uC, respectively, suggesting that these conserved miRNAs play similar roles due to their high stability in seed-target duplexes.
Calculation of miRNA thermodynamic parameters. Melting temperature (Tm) of each miRNA duplex and seed-target duplex were predicted by means of nearestneighbor model 54 . The formula for the calculation is as follows. Statistical analysis. Student t-test was carried out for assessing the correlation between relative luciferase activity and Tm 2-8 , miTm 1-5 (Tm 1-5 ), or Tm 2-8 2 0.53 3 miTm 1-5 (Tm 1-5 ) value, and that between Tm 2-8 and miTm 1-5 (Tm 1-5 ). Kolmogorov-Smirnov (K-S) test was carried out to validate the difference of significance in Tm values of miRNAs of each species.