Challenges in microRNAs' targetome prediction and validation.

MicroRNAs (miRNAs) are small RNA molecules with important roles in post-transcriptional regulation of gene expression. In recent years, the predicted number of miRNAs has skyrocketed, largely as a consequence of high-throughput sequencing technologies becoming ubiquitous. This dramatic increase in miRNA candidates poses multiple challenges in terms of data deposition, curation, and validation. Although multiple databases containing miRNA annotations and targets have been developed, ensuring data quality by validating miRNA-target interactions requires the efforts of the research community. In order to generate databases containing biologically active miRNAs, it is imperative to overcome a multitude of hurdles, including restricted miRNA expression patterns, distinct miRNA biogenesis machineries, and divergent miRNA-mRNA interaction dynamics. In the present review, we discuss recent advances and limitations in miRNA prediction, identification, and validation. Lastly, we focus on the most enriched neuronal miRNA, miR-124, and its gene regulatory network in human neurons, which has been revealed using a combined computational and experimental approach.

The Canonical MicroRNA Processing Machinery

In recent years, the central role of microRNAs (miRNAs) as post-transcriptional repressors has been extensively studied in development and disease. In animals, miRNAs are non-coding 22 nucleotide (nt) sequences that regulate gene expression by binding either partially- or fully-complementary RNA sequences - predominantly messenger RNA (mRNA), but also but also ribosomal RNAs, transfer RNAs, long intergenic non-coding RNAs, and other miRNAs (Helwak et al., 2013). The main determinants of target recognition reside within the miRNA 5′ end in positions 2 to 7, which together are commonly known as the “seed” sequence (Ha and Kim, 2014). Computational predictions suggest that individual miRNAs regulate large genetic networks by targeting hundreds of mRNA molecules (Agarwal et al., 2015), although experimental data on entire networks within specific cell types are rare. Most experimental studies only validate one or a few miRNA targets at once. In turn, an individual mRNA molecule can possess binding sites for several miRNA species (Friedman et al., 2009; Chou et al., 2018), though the latter might compete against each other for specific binding sites (Saetrom et al., 2007) and thereby lead to non-synergistic outcomes (Xu et al., 2011). Additionally, the existence of competing endogenous RNAs (Tay et al., 2014; Denzler et al., 2016), i.e. natural miRNA sponges such as long non-coding and circular RNAs, is further evidence of the intricacy of miRNA regulatory networks.

The biogenesis and maturation of miRNAs is highly conserved across species. miRNAs are generated from partially-complementary regions of primary RNA transcripts (pri-miRNA) produced mainly by RNA polymerase II, but also by RNA polymerase III. These transcripts anneal and give rise to short RNA hairpins or stem-loops (Ha and Kim, 2014). Initially contained within non-coding RNAs or within the introns of coding sequences (host genes), these hairpins are recognized and cleaved in the mammalian cell nucleus by the miRNA-processing complex, also known as microprocessor, formed by the Drosha ribonuclease and the DiGeorge critical region 8 (DGCR8) protein. By cutting 11 base pairs (bp) away from the base of the hairpin stem, the miRNA-processing complex produces a 70-nt-long sequence known as a precursor miRNA (pre-miRNA) and characterized by a 5′ phosphate and a 2-nt 3′ overhang (Han et al., 2004). Pre-miRNA hairpins are subsequently transported by exportin 5 to the cytoplasm, where they are cleaved by the Dicer endoribonuclease 22 nt away from the Microprocessor cutting site. This second cut generates a new 5′ monophosphate and one more 2-nt 3′ overhang on the opposite side of the double-stranded RNA sequence (Chiang et al., 2010). The resulting ~22-nt RNA duplexes are incorporated into the Argonaute protein and unwound to form, together with other proteins, the RNA-induced silencing complex (RISC). RISC then mediates the miRNA-guided scanning of target transcripts and the subsequent miRNA-target pairing (Bartel, 2009; Ha and Kim, 2014). Ultimately, miRNA-targeted mRNA molecules are degraded or exhibit reduced translational efficiency (Krol et al., 2010).

Challenges in the Identification of Novel MicroRNAs

Although miRNA identification was initially achieved by Sanger sequencing of complementary DNA clones produced from small-RNA libraries, the rapid advance in sequencing technologies towards higher-throughput methods has led to an exponential increase in the number of predicted miRNAs (Kozomara and Griffiths-Jones, 2011). RNA sequencing (RNA-Seq), in particular, has revolutionized the field of transcriptomics and spearheaded the discovery of all types of novel transcripts (Kukurba and Montgomery, 2015). Offering higher coverage and resolution than traditional Sanger- and microarray-based sequencing methods, RNA-Seq experiments have revealed many novel miRNAs (Kozomara and Griffiths-Jones, 2011): this has resulted in an enormous challenge in miRNA annotation and validation. The latter is further compounded by the usage of heterogeneous RNA-Seq methods and deep-sequencing parameters. Similarly, low RNA quality has a notoriously detrimental effect on sequencing results. Sequencing low-quality RNA has the potential to bias miRNA quantifications and might lead to assuming that RNA fragments resulting from degradation are real miRNAs (Kukurba and Montgomery, 2015; Ludwig et al., 2017). Indeed, hundreds of miRNAs can be erroneously annotated from a single, poorly-analyzed RNA-Seq dataset and decimate the fraction of legitimate miRNAs annotated in databases (Kozomara and Griffiths-Jones, 2014). In order to deal with these difficulties, guidelines for miRNA annotation have been devised and are continuously being refined.

The first set of guidelines for annotating miRNAs was published in 2003 (Ambros, 2003). These guidelines stated that cloning, sequencing, or northern blotting should be used to provide evidence for the existence of a 22-nt-long sequence, and that the predicted origin of such sequence should be a region with the potential to generate a stem-loop structure. These principles, however, are subject to certain limitations. For instance, millions of hairpin structures are predicted to exist in the genome, yet not all of them give rise to miRNAs (Kozomara and Griffiths-Jones, 2011). Hence, stem loops are necessary but insufficient to predict the existence of a “true” miRNA. Furthermore, even “true” miRNA candidates might escape detection in sequencing experiments if tissues or cells are examined at time-points in which their expression levels are low. For this reason, complementing criteria for the validation of animal miRNAs were established. These criteria are often incorporated in modern miRNA identification tools such as miRDeep (An et al., 2013) and include: a) the presence of matching 3′ overhangs up to 4 nt-long; b) 5′ ends compatible with the miRNA processing machinery; c) a maximum number of permitted mismatches between pre-miRNA strands; d) a minimum number of required mapped reads; and e) a specific pattern for such reads (Friedlander et al., 2008; Hendrix et al., 2010). With regard to criterion e), a high number of reads is expected to map to the mature miRNA, whereas much lower read numbers should map to its loop region and to the second strand of the miRNA duplex (miR*; also known as “passenger” strand). Additionally, reads mapping to both loop and miRNA (or to both loop and miR*) are expected to be almost undetectable as they contain the cutting sites for the miRNA processing machinery. Although these principles are certainly useful in aiming for high-confidence annotations, defining a minimum number of required mapped reads as a guideline could cause problems, particularly because some miRNA species are expressed only at low levels and in a temporally- and tissue-restricted manner (Karali et al., 2010). Numerous candidate miRNAs have been annotated in miRNA repositories even with a single mapped read from a deep-sequencing experiment as evidence (Berezikov et al., 2006). Since many of these miRNAs have not yet been validated, it is still possible that they are artifacts.

“True” MicroRNAs: Needles or Haystack?

Several of the hurdles inherent in the identification of novel miRNAs are evident in multiple miRNA databases. miRBase, for instance, represents the primary public repository and online resource for miRNA sequences and annotations, and describes 48,860 mature miRNAs from 271 organisms in its 22nd version (Kozomara et al., 2019). Nonetheless, many miRNA annotations in this and other databases correspond to predicted miRNAs and might thus be inaccurate. As an example, over 150 annotated miRNAs have failed to yield RNAs with miRNA features in response to the overexpression of their hairpin sequences flanked on both sides by the ~100-nt genomic sequence required for their correct processing (Chiang et al., 2010). Moreover, although the mapping patterns of over 5.5 billion reads gathered from almost 1500 small RNA-Seq datasets support the validity of the miRBase annotations to varying degrees (up to 65% depending on the species), over 200 miRNA sequences have been removed from this database in the past 5 years (Kozomara et al., 2019). In contrast, more than 10,000 miRNAs were added within the same timeframe.

For homo sapiens, in particular, the miRBase presently contains 2654 mature miRNA sequences predicted to originate from 1917 hairpin precursors. Yet, only approximately one quarter of these miRNAs are annotated with high confidence, and no more than 20 reads are associated with each strand of the RNA duplex for about half of them (Kozomara et al., 2019). However, validating all candidate miRNAs via overexpression assays (Chiang et al., 2010) is a titanic task likely to fail in identifying miRNAs produced via non-canonical pathways and hence not complying with the features that arise from conventional miRNA processing (Yang and Lai, 2011; Ha and Kim, 2014; Kozomara et al., 2019). Thus, improved annotation approaches, and the development and implementation of stringent validation pipelines are imperative for defining whether a miRNA is “true” or falsely-identified.

Experimental Approaches to Identify and Validate MicroRNAs

As most of the RNA species obtained from biological tissues or cells in culture correspond to ribosomal RNA, a variety of depletion and enrichment strategies have emerged over time to preferentially sequence RNAs of interest. A prominent strategy to enrich miRNAs from total RNA is to isolate small RNAs, i.e. species of 15–30 nt. This is commonly achieved either by excising gel fragments after electrophoresis in the region corresponding to the expected RNA fragment size, or by using silica columns that bind small RNAs (Kukurba and Montgomery, 2015). These approaches, however, are prone to miRNA misidentification in the case of low-quality or degraded RNA. Recently, a number of methods based on the crosslinking of miRNAs and their targets, together with immunoprecipitation of RISC proteins, have emerged for the validation of miRNAs and for the identification of miRNA-target interactions. High-throughput sequencing of RNAs isolated by crosslinking and immunoprecipitation of Argonaute, for instance, allows miRNAs contained within AGO and the binding sites of Argonaute-miRNA complexes within distinct mRNA molecules to be identified (Chi et al., 2009; Boudreau et al., 2014). Building on this method, photoactivatable ribonucleoside-enhanced crosslinking and immunoprecipitation uses randomly-incorporated photoactivatable nucleosides, commonly 4-thiouridine, prior to crosslinking RNA molecules with RNA binding proteins and ribonucleoprotein complexes (Hafner et al., 2010). This enables enhanced crosslinking and subsequent RNA recovery with reduced ultraviolet irradiation and, more importantly, leads to characteristic mutations, depending on the nucleoside utilized, which correlate with RNA-protein crosslinking sites. Thus, by coupling photoactivatable ribonucleoside-enhanced crosslinking and immunoprecipitation with AGO immunoprecipitation, miRNA binding sites can be identified by mutational analyses, and discriminated from contamination by other RNA species. Similarly, by crosslinking the RNA sequences bound to the AGO-miRNA complex after immunoprecipitation and sequencing the hybrids generated, CLASH (crosslinking, ligation, and sequencing of hybrids) allows miRNAs and their targets to be identified with high confidence (Helwak et al., 2013; Helwak and Tollervey, 2014). Importantly, CLIP-based sequencing approaches, including CLASH, are among those considered reliable for miRNA-target validation within the miRTarBase (version 7.0) (http://mirtarbase.mbc.nctu.edu.tw/php/index.php), the most comprehensively annotated database of experimentally-validated miRNA-target interactions (Chou et al., 2018). Yet, in spite of its benefits, crosslinking poses the risk of generating undesired mutations, introducing background, and failing to detect mRNA-miRNAs interactions occurring via non-canonical “seed” pairing. Further, CLIP-based methods often require large amounts of biological material and are biased towards the identification of miRNAs that are highly expressed and that contain specific sequences (Tan et al., 2014). To overcome these limitations, RNA immunoprecipitation (RIP) followed by next generation sequencing (RIP-Seq) enables specific RNA-binding proteins and their target-RNA complexes to be pulled down either in their native state or after cross-linking (Malmevik et al., 2015; Gagliardi and Matarazzo, 2016). In this context, immunoprecipitation of AGO-bound RNAs in the absence of crosslinking, followed by sequencing (AGO-RIPseq), allows miRNAs and bound mRNAs to be identified (Petri and Jakobsson, 2018). RIP-Seq datasets are considered as evidence of miRNA-mRNA binding in the DIANA-TarBase (version 8.0) (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=tarbase/index), a manually-curated database containing over 790,000 experimentally-validated miRNA-gene interactions (Karagkouni et al., 2018). This repository also incorporates the results of the IMPACT-Seq approach (identification of miRNA-responsive elements by pull-down and alignment of captive transcripts-sequencing), in which RNA targets are sequenced after co-purification with transfected miRNA mimics biotinylated at their 3’ end in the absence of crosslinking (Tan et al., 2014). Although these methods have enabled an enormous number of novel miRNA−mRNA interactions to be identified, such interactions remain limited by the expression levels of the miRNAs of interest within particular cells or tissues at the examined time-point.

Predicting the Targetome and Genetic Regulatory Network of a MicroRNA

The targetome of a miRNA refers to the set of RNA molecules that it targets. Several computational approaches have been developed to predict the targetome of specific miRNAs. To date, most of these prediction algorithms use the canonical rule of miRNAs binding to their target mRNAs mainly by base pairing between the seed region of the miRNA and fully-complementary sequences within the 3′ untranslated region of the mRNA (Tan et al., 2014). However, evidence indicates that bulges, G:U wobbles, “seedless” interactions, and binding regions within a gene’s coding sequence are admissible and not uncommon, with canonical miRNA-mRNA interactions accounting for only 25% (Helwak et al., 2013). Hence, computational prediction of binding regions within mRNA molecules, and thus of miRNA target genes, remains a major challenge. Although non-canonical miRNA-mRNA interactions seem to have no effect on mRNA or protein levels in mammals (Agarwal et al., 2015), incomplete miRNA-mRNA base-pairing is common and enhances Argonaute-mediated target slicing and the rate of product release (or target turnover) in other classes of living organisms (Wee et al., 2012; Helwak et al., 2013; Chen et al., 2017). Thus, although the biological relevance of accurately predicting partially-complementary miRNA-mRNA interactions remains controversial, there is no doubt such interactions are important to ensure undesired mRNA regulation is avoided upon administration of RNA interference-based therapeutics.

In general, the prediction accuracy for individual miRNA targetomes is relatively high if distinct parameters are considered. For example, to predict the probability of a specific miRNA-mRNA, TargetScan has developed the context++ model to go beyond simply searching for conserved multimers (more specifically, sequences of 6, 7, or 8 nt) that match the miRNA seed region. The TargetScan context++ model factors in site accessibility (Robins et al., 2005; Tafer et al., 2008; Hölscher, 2014), and allows mismatches in the miRNA seed region (Agarwal et al., 2015) by incorporating the concepts of conserved 3′ pairing (Friedman et al., 2009) and centered sites (Shin et al., 2010). Moreover, many online platforms now even incorporate machine-learning algorithms and advanced text-mining approaches using natural language processing methods to explore the literature and devise distinct miRNA-targeting mechanisms, the probability that they occur, and their effects in terms of translational repression and mRNA degradation. Offering high-quality miRNA targetome predictions for distinct vertebrate species and model organisms, such platforms include the user DNA Intelligent Analysis (DIANA)-microT (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=microT_CDS/index) (Paraskevopoulou et al., 2013), miRDB (mirdb.org), and MIRZA-G (Gumienny and Zavolan, 2015). Although targetome and large-scale genetic network predictions for individual miRNAs are increasingly reliable, accounting for the combined action of multiple miRNAs and predicting their joint effects on mRNA or protein levels remain daunting endeavors. The distance between miRNA binding sites in an mRNA molecule, for example, defines whether different miRNAs will act cooperatively or competitively, with distant sites favoring the former, but short distances promoting the latter (Grimson et al., 2007). In this sense, the DIANA-miRPath online platform enables predictions on the combined effect of modulating multiple miRNAs to be run (Vlachos et al., 2012, 2015). Altogether, only the combination of prediction algorithms and experimental approaches will ultimately lead to a detailed understanding of miRNA regulatory networks.

The Role of MicroRNAs in the Central Nervous System

miRNAs are implicated in virtually all aspects of cell biology, and play essential roles in several developmental and pathological processes. Dysfunctional miRNA regulatory networks correlate with the emergence and/or progression of diabetes, heart failure, cancer, and neurodevelopmental disorders (Allen, 1940; Blenkiron and Miska, 2007; Rajman and Schratt, 2017). In the brain, for instance, miRNAs act not only as “fine-tuners” but also as master regulators of neuronal circuit development, maturation, and function, and influence processes such as cell-fate determination, cell migration, neuronal polarization, cognition, and synapse formation and plasticity (Ye et al., 2016; Rajman and Schratt, 2017). Early studies aimed at assessing the role of miRNAs in the central nervous system revealed that deletions of key components of the miRNA-processing pathway, including Dicer, DGCR8, and Argonaute, lead to defects in neurogenesis, impaired brain development, and neuronal atrophy (Schaefer et al., 2007; Davis et al., 2008, 2015; Im and Kenny, 2012). Moreover, Cre-loxP-mediated conditional deletion of Dicer in specific neuronal cell types leads to their premature death, often as a consequence of a depleted neuronal progenitor pool or of enhanced apoptosis (Barca-Mayo and De Pietri Tonelli, 2014). Thus, miRNAs are likely to regulate cell-type specification in the central nervous system. Remarkably, transcriptomic profiling studies have identified miR-124, in particular, as the most abundant small RNA in the central nervous system (Lagos-Quintana et al., 2002; Landgraf et al., 2007).

While ectopic expression of miRNA-124 in mice depletes the pool of neuronal progenitors within the subventricular zone by forcing their differentiation, its partial inhibition blocks neurogenesis and promotes progenitor self-renewal (Cheng et al., 2009; Åkerblom et al., 2012). Nonetheless, as the cellular and molecular composition of the human brain is different to that of other mammals, including primates (Somel et al., 2011), investigating the role of miRNAs in human model systems is indispensable. In the HeLa human cell line, miR-124 overexpression has previously been reported to drive the acquisition of a neuronal-like transcriptional profile (Conaco et al., 2006), strongly hinting at the crucial role of this miRNA species in also defining neuronal identity in humans. In order to gain a deeper insight into the role of miR-124 in human neurogenesis, we used a highly homogeneous cellular model system known as iNGN cells. In this system, overexpression of the transcription factors neurogenin-1 and neurogenin-2 is driven by a doxycycline-inducible promoter in induced pluripotent stem cells. These cells acquire a homogeneous bipolar neuronal morphology and become electrically active 7 days after doxycycline administration (Busskamp et al., 2014). More importantly, as in the murine central nervous system, miR-124 is also the most abundant miRNA species in iNGN cells. In our study, CRISPR/Cas9-mediated disruption of the six miR-124 paralogs revealed this miRNA species to not be required for neurogenesis (Kutsche et al., 2018). However, miR-124-depleted iNGN cells acquired more complex morphologies, reduced their dendritic extension speed, and, unexpectedly, shifted from a glutamatergic identity towards the preferential use of acetylcholine as a neurotransmitter, which are findings consistent with previous reports (Yu et al., 2008; Franke et al., 2012; Volvert et al., 2014). Of note, miR-124-deficient cells also exhibited increased apoptosis in long-term cultures, as was previously shown in mouse neurons in vivo upon knock-out of a single miR-124 locus (Sanuki et al., 2011). Since the molecular machinery mediating these phenotypic alterations is so far largely unexplored, we decided to combine experimental and computational approaches to deepen our understanding of the miR-124 regulatory network.

Taking its limitations into consideration, we used AGO-RIP-Seq in miR-124-depleted iNGN cells at distinct time-points over the course of differentiation to identify miR-124 high-confidence targets (Kutsche et al., 2018). Although 4024 human transcripts have been described as possessing binding sites for miR-124, our approach allowed us to reliably identify 127 transcripts with reduced expression in miR-124-deficient cells compared to their wildtype counterparts. Remarkably, 98 of these 127 transcripts (77%) had been predicted to be miR-124 targets, and 38 of them had also been annotated as previously validated in the miRTarBase (Chou et al., 2018). Using luciferase reporter assays, which are considered a reliable validation method in the miRTarBase, we confirmed 43 novel bona fide miR-124 targets. Despite this approach signifying important progress in the validation of mRNAs targeted by miR-124, its interaction capacity with the remaining ~4000 annotated transcripts needs to be confirmed, potentially by also considering changes in protein levels. However, targeting 98 mRNAs, of which 81 have now been validated, demonstrated that miR-124 exerts its regulatory function at the network level. Another striking observation in our study was that 24 of these 98 predicted miR-124 targets were transcription factors (TFs) (Kutsche et al., 2018), hinting that these TFs are potential effectors of miR-124. As TFs modulate entire genetic programs and thereby have a strong impact in cell fate decisions (Hobert, 2008; Lambert et al., 2018), miRNA-mediated perturbations in the expression levels of TFs are particularly relevant.

Knowledge of the existence of regulatory feedback mechanisms between miRNAs and TFs is not new. The die-1 TF, required for controlling chemosensory laterality in Caenorhabditis elegans, was shown over a decade ago to perform its functions by activating the lsy-6 miRNA in left-sided chemosensory neurons but having its expression levels modulated in left-sided neurons by mir-273 (Chang et al., 2004). Similarly, multiple miRNA-TF feedback circuits have been described as essential for mammalian neurogenesis, including miR-133b/Pitx (Kim et al., 2007), miR-7/Pax6 (de Chevigny et al., 2012), miR-132/Nurr1 (Nr4a2) (Yang et al., 2012), and miR-17-3p/Olig2 (Chen et al., 2011). Moreover, specific miRNAs are known to regulate the expression of multiple TFs, as is the case for the neuron-enriched miR-9 in the murine brain (Shibata et al., 2011). As a consequence of their biological importance and widespread occurrence, multiple online platforms have been developed to predict miRNA-TF interactions, including TFmiR (Hamed et al., 2015) and DIANA-mirExtra (Vlachos et al., 2016). However, building an interaction network is not trivial. To generate gene regulatory networks, in particular, co-expression networks are widely used (Yang et al., 2014; van Dam et al., 2018). Such networks are characterized by connected pairs of nodes representing significantly correlated expression patterns which can also contain information regarding both the direction and the strength (or weight) of the interaction. Nonetheless, network analyses are often performed using defined subsets of entities, leading to biased results if relevant interactions via other intermediaries are not considered. To overcome this limitation, the weighted topological (wTO) network analysis method implicitly includes correlations among nodes which are subsequently omitted in the analysis (Nowick et al., 2009; Gysi et al., 2018). Thereby, the wTO approach enables a more reliable evaluation of the interactions among sets of nodes of interest relative to other network analysis pipelines. By implementing the wTO network analysis method on the differentially regulated TFs in miR-124-deficient iNGN cells, we identified a number of factors likely to mediate the biological functions of miR-124. To validate the results, the function of the zinc finger 787 (ZNF787) was additionally characterized in iNGN cells, confirming its role as a repressor of neuronal features and revealing that it modulates a big fraction of the targets regulated by miR-124 (Kutsche et al., 2018).

Lastly, the experimental deletion of a highly abundant miRNA species did not result in a vacuum of miRNA regulation, as we detected other miRNAs being upregulated (Kutsche et al., 2018). Hence, our data suggest that this de novo miRNA regulation upon manipulation affects the primary phenotype. This likely contributes to the difficulties in interpreting experimental miRNA knockout and knockdown phenotypes, and must therefore be considered.

Outlook and Perspective

With less than 2% of the human genome encoding protein-coding genes, non-coding RNAs represent most of the human transcriptome (Tay et al., 2014). The central roles of such RNAs, including miRNAs, in health and disease is now widely recognized. Consequently, steady progress has been made over recent years in the annotation and validation of novel miRNA. However, there are still a number of important challenges in collecting highly heterogeneous data from multiple experimental sources and in developing highly accurate prediction tools. The latter are essential for enhancing the medical potential of miRNAs (Jimenez-Mateos et al., 2012). We foresee that machine learning will make an important contribution to the refinement of miRNA repositories and prediction tools acting at the systems level.