Composition of seed sequence is a major determinant of microRNA targeting patterns

2.1 Data retrieval

The CLASH dataset was available as Supplementary Data in (Helwak et al., 2013). This dataset was downloaded from the journal’s website and analysed in our study. In this published study, HEK293 cells overexpressing human AGO1 were crosslinked by ultraviolet irradiation. The crosslinked RNA-induced silencing complex (RISC), containing both miRNAs and cognate mRNA target transcripts, were immunoprecipitated using an antibody that specifically binds to the tagged hAGO1 in the RISC. RNA–RNA ligation was then performed for the isolated RISC complexes to covalently ligate the miRNAs and cognate mRNA target transcripts. Ligated RNA products were finally used as a template for RNA-seq library construction, and the identity of both miRNAs and their paired mRNA target transcripts were determined by sequencing chimeric miRNA-target RNA. Altogether, 18 514 miRNA-target pairs were identified in this way.

Microarray data from the photoactivatable-ribonucleoside (PAR)-CLIP study (Hafner et al., 2010) were used to evaluate the impact of miRNAs on target RNA expression regulation. In this microarray study, 25 miRNAs were concurrently suppressed in HEK293 cells by a pool of antisense oligonucleotide inhibitors, and the impact on target expression regulation at the transcriptional level was determined by Affymetrix Human U133Plus2 chips. Raw microarray data were downloaded from the NCBI GEO database (Barrett et al., 2013), and array signals were normalized using the robust multiarray averaging method from the Bioconductor package (http://www.bioconductor.org). Signals from transcripts of the same gene were combined and averaged using Gene mapping index files downloaded from the NCBI ftp site (Maglott et al., 2011). The log2-transformed fold change of each gene after miRNA depletion was calculated by comparing with the negative control arrays (mock transfections).

Proteomic data from (Selbach et al., 2008) were retrieved from the online Supplementary File at the journal’s website, and used to evaluate the impact of miRNA on target protein expression. In the Selbach study, five miRNAs (let-7b, miR-1, miR-16, miR-30a and miR-155) were individually overexpressed in HeLa cells and changes in protein synthesis were determined with mass spectrometry-based proteomic analysis. Protein expression changes were represented by log2-transformed fold change as compared with the negative control (mock transfections).

2.2 Data analysis

Thermodynamic stability of miRNA-target RNA duplex, as represented by the ΔG value, was calculated with RNAfold (Hofacker, 2003). The stability of perfect base pairing between miRNA seed/non-seed sequence and the target binding site was calculated with the nearest neighbor method and thermodynamic parameters for RNA–RNA base pairing (Xia et al., 1998). Canonical seed is defined as any 6-mer sequence within positions 2–8 or 7-mer sequence within positions 1–8 of the miRNA; extended seed is defined as any 6-mer sequence within miRNA positions 4–10; and non-seed sequence is any sequence downstream of the extended seed region. Sequences of the miRNAs and their cognate target sites from the CLASH dataset were compared to identify base pairing patterns, including perfect pairing, imperfect pairing with a GU mismatch and imperfect pairing with a non-GU mismatch. As negative control, the target transcript from each miRNA-target pair was replaced by a randomly selected transcript from the CLASH dataset (shuffled control), and stochastic pairing pattern between the miRNA and random transcript was determined. The enrichment ratio of each base-pairing pattern between the miRNAs and paired target sites was calculated by comparing with the negative control. Statistical difference between various targeting patterns and the negative control for thermodynamic stability or gene expression change was calculated with Student’s t-test.

3.1 Canonical seed pairing was absent in most miRNA target sites

Using the CLASH method, 18 514 pairs of miRNAs and cognate target binding sites were identified (Helwak et al., 2013). Here, we reanalysed this dataset to identify sequence features that are significantly associated with miRNA targets. One well-characterized feature is the perfect pairing between canonical miRNA seed (positions 2–8) and the target binding site. The 3′ portion of the miRNA was generally considered to play an auxiliary role, with perfect seed pairing as the dominant determinant (Bartel, 2009). To reassess the relative importance of both seed and non-seed sequences in target recognition, all possible 6-mers from the miRNA sequences were compared with the corresponding paired target sites. Specifically, for each miRNA-target pair, the entire miRNA sequence was scanned to identify all position-specific 6-mers that were matched to the target sequence. As negative control, each miRNA sequence was compared with a randomly assigned transcript from the CLASH dataset (shuffled control). Overall, 6-mers from positions 3–8 were most enriched, which fall within the canonical seed region (Fig. 1). The other 6-mer position from the canonical seed region (nucleotides 2–7) was ranked as the third most enriched along the miRNA sequences. Thus, our analysis was consistent with previous findings that perfect canonical seed pairing is an important parameter for miRNA target recognition.

Interestingly, 6-mers from all positions of the miRNA sequence were paired to the target sites at significantly higher rates as compared with the negative control (Fig. 1). In particular, 6-mers from two positions, nucleotides 4–9 and 5–10 were the second and fourth most enriched, respectively, although they fall outside the canonical seed region. Thus, nucleotides 4–10 were considered as the “extended” seed region in our analysis as 6-mers from this region partly overlap with canonical seeds, matching to the target sites at similarly high rates to canonical seeds. Further, 6-mers from all other positions in the miRNA sequence were also significantly enriched, although to a lesser extent as compared with canonical or extended 6-mer seeds (Fig. 1). As shown in Table 1, perfect matches to canonical 6-mer seeds were identified in 22.4% of all target sites, which was 6.8-fold enriched as compared with the shuffled control. Similarly, 6-mers from the extended seed region were matched to 20.9% of all target sites, which was 6.3-fold enriched. Seed pairing was even more significant when 7-mer seeds were considered, with 11.0- and 11.5-fold in enrichment for canonical and extended seed regions, respectively. When 6-mers from both canonical and extended seed regions were considered collectively, seed matches were found in 32.9% of all target sites. Thus, the majority of the target sites did not have perfect pairing to any miRNA seed even when both canonical and extended seeds were considered.

Similar to perfect base pairing, G-U wobble pairs can also provide thermostability in the RNA duplex. However, a G-U mismatch in the seed region may potentially interfere with target binding based on structural analysis (Bartel, 2009). The impact of G-U mismatch on seed-target binding is still a debatable topic, as some studies suggested that the presence of a G-U mismatch would mostly abolish the miRNA targeting ability, while some other studies identified functional target sites with G-U mismatch to the miRNA seed (Brodersen and Voinnet, 2009; Vella et al., 2004). This issue was revisited using the CLASH data to analyse any 6-mer seed with a single G-U mismatch to the target site. As shown in Table 1, there were significantly more 6-mers from both canonical and extended seed regions, with 2.3- and 2.4-fold in enrichment, respectively, for target pairing with a G-U mismatch as compared with the shuffled control. Thus, our analysis indicated that imperfect seed pairing with a G-U mismatch to the target site may still be biologically relevant. However, the enrichment ratios were much lower than those of perfect seed match, suggesting significantly less important role of G-U mismatched seed pairing in miRNA target recognition. Furthermore, the impact of non-GU mismatches was also examined. A 7-mer with an internal non-GU mismatch was only slightly enriched in the target sites (1.3- and 1.4-fold in enrichment, Table 1), indicating that a non-GU internal mismatch in seed pairing would mostly abolish miRNA target recognition. Combined together, 51.1% of all target sites were paired perfectly to canonical/extended seeds, or imperfectly with a G-U mismatch.

We further analysed the relevance of non-seed sequence of miRNA in target recognition. Among all seed-pairing target sites (9465 in total), 43.5% of them were also paired to non-seed region of miRNA, which was a significant enrichment over shuffled control (Table 2). Furthermore, as to non-canonical target sites without any seed match, 40.9% were paired to the non-seed region of the miRNA. Thus, pairing to non-seed region of the miRNA was a significant feature for both canonical and non-canonical targets.

3.2 GC content of miRNA seed was strongly correlated to target recognition patterns

All 18 514 target sites from the CLASH dataset were analysed to identify sequence features that characterized perfect seed pairing. First, canonical 6-mer seeds, involving either positions 2–7 or 3–8, were characterized. Interestingly, there was a very strong correlation between GC content and perfect pairing of 6-mer seeds within the canonical region, with correlation coefficients (r) of 0.91 and 0.96 for nucleotides 2–7 and 3–8, respectively (Fig. 2A and B). Few targets were paired to any canonical seed when the GC content of the seed was low. On the other hand, a much higher percentage of seed-pairing targets were found for miRNAs with high GC content of the seed. GC content is commonly used to represent thermodynamic stability of RNA duplex. Thus, the thermostability of seed-target binding was also evaluated by calculating the free energy (ΔG) of seed-target RNA duplex with RNAfold (Hofacker, 2003). As expected, the ΔG of seed-target duplex was strongly correlated with the percentage of targets with seed pairing (r = −0.90 and −0.95 for nucleotides 2–7 and 3–8, respectively; Fig. 2C and D). Thus, seed-based target recognition was dependent on the thermostability of seed pairing to the target site. For miRNAs with weak canonical seed-pairing stability, non-canonical targeting was the dominant mechanism for target recognition.

Besides targets involving canonical seed pairing, we also analysed non-canonical targets including pairing to extended seeds, imperfect pairing with a G-U mismatch or pairing to non-seed sequences. Similar to canonical seed pairing, 6-mer matches to the extended seed region (positions 4–10) were strongly correlated to sequence GC content (r = 0.99, Fig. 3A), which reflected the stability of seed-target duplex as represented by ΔG (r = −0.98, Fig. 3B). Combined together, the GC content of both canonical and extended seeds within the full seed region was a strong predictor of perfect seed pairing to the target site. In addition to perfect seed pairing, imperfect seed pairing with a G-U mismatch was also analysed. As shown in Figure 3C, low GC content (≤50%) of the seed was also a significant predictor of poor seed pairing (r > 0.99). Interestingly, no such correlation was observed when the GC content was higher (>50%). Pairing between non-seed 6-mers in the miRNA sequence and the target sites was also strongly dependent on the GC content of the 6-mers, in a similar way to seed pairing of the target site (r = 0.98, Fig. 3D). Thus, binding stability of the paired nucleotides was a major determinant of miRNA targeting patterns, whether involving seed or non-seed sequences.

3.3 Canonical and non-canonical targeting had similar thermostability but distinct impacts on target expression

Overall thermostability of the miRNA-target RNA duplex was assessed for both canonical and non-canonical targeting patterns. The following non-canonical target types were included in the analysis: targets pairing to extended seed region, pairing to non-seed region and with no pairing anywhere in the miRNA sequence. These non-canonical targets were compared with canonical targets as well as randomly assigned non-target transcripts from the CLASH dataset (shuffled control). Specifically, overall distribution of thermostability for each type of target binding, as represented by ΔG, was determined with RNAfold (Hofacker, 2003). As shown in Figure 4, binding of miRNAs to their cognate targets was significantly more stable than binding to randomly matched transcripts (P < 10⁻³⁰⁰ with Student’s t-test for all target types). The mean binding free energies were −24.6, −24.8, −24.4 and −20.9 kcal/mol for targets involving canonical seed pairing, extended seed pairing, non-seed 6-mer pairing and no pairing anywhere in the miRNA sequence, respectively. There was no difference among the first three target types, all of which had significantly lower free energy than the last type. In contrast, binding between miRNAs and randomly matched transcripts were much more unstable, with mean ΔG of −17.1 kcal/mol.

It was shown previously that binding of miRNA to the target generally resulted in suppressed gene expression at both the transcriptional and translational levels (Baek et al., 2008; Selbach et al., 2008). However, most previous studies were focused on canonical miRNA targeting involving perfect seed pairing, and the impact of miRNA binding to non-canonical targets has not been well characterized. In our analysis, the CLASH dataset was combined with a public microarray dataset, which studied the impact of concurrent depletion of 25 miRNAs on mRNA levels in HEK 293 cells (Hafner et al., 2010). As negative control of the analysis, each target transcript in the miRNA-target pairs from the CLASH dataset was replaced by a randomly assigned non-target transcript with similar expression level to the target as determined by microarrays. Specifically, the expression difference between a target transcript and its random control was <0.5 (log2-transformed expression value). Consistent with the Hafner study, canonical targets with perfect seed pairing sites were significantly upregulated by miRNA depletion as compared with the negative control (mean fold change of 0.044 and P < 10⁻¹⁷ by Student’s t-test, Fig. 5A). In contrast, all three types of non-canonical targets were upregulated only moderately as compared with the negative control. The mean fold changes for individual non-canonical target types were in the range of 0.014 to 0.021 and P-values in the range of 0.003 to 0.0002. Thus, miRNA binding to target transcripts led to distinctively different outcomes for gene expression regulation at the mRNA level, which were dependent on the type of miRNA targets.

Furthermore, the impact on target protein expression by different miRNAs was evaluated with proteomic data. Selbach et al. performed a proteomic study to identify the global impact of miRNA overexpression on protein synthesis (Selbach et al., 2008). In the Selbach study, five miRNAs were individually transfected into HeLa cells and changes in protein synthesis due to miRNA overexpression were determined by proteomic studies. In our analysis, different miRNA targeting types were evaluated with this proteomic dataset. Specifically, all miRNA-target pairs from the CLASH dataset were cross-checked with the proteomic dataset, and 330 common miRNA-target pairs presented in both datasets were identified. From these identified pairs, different miRNA targeting patterns were correlated to changes in target protein expression due to overexpression of individual miRNAs. As shown in Figure 5B, canonical targets with perfect 6-mer seed match were significantly suppressed by miRNA overexpression (mean fold change of −0.114 and P = 0.0004 by comparing with the negative control by Student’s t-test). In contrast, none of the other miRNA targeting patterns were significantly associated with suppression of target protein expression (mean fold change ranging from −0.027 to 0.001 and P-values ranging from 0.36 to 0.87 by Student’s t-test).