Integrated Analysis of microRNA Targets Reveals New Insights into Transcriptional-Posttranscriptional Regulatory Cross-Talk

1. Introduction

Over the last 20 years, genome projects have shown that the complexity of the most evolved organisms is not explained when looking at protein-coding genes, as they are highly conserved even in distant species. Large genomic regions are occupied by genes for the so-called non-coding RNAs, which perform their function without being translated into amino acid chains, and that may contribute to modulate the variability among different organisms. Among them, microRNAs have been shown to be involved in regulatory functions since their first discovery in 1993, and to be an essential component of the intricate network that regulates gene expression in many cell processes [1,2]. Their mature form is a very short double strand RNA, composed of a passenger strand that is then misplaced, and a guide strand that identifies the target, commonly through base pairing between bases in positions 2-7 of the microRNA (the seed) and short regions in the mRNA, often located in 3’ UTR of the transcript [3,4]. The peculiar shape of microRNAs has promoted the development of many computational approaches, both to discover them in the genome, and to predict their targets. Currently, more than a hundred target prediction tools are available, based on sequence complementarity and other features, as described in several reviews [5,6]. According to these bioinformatic analyses, almost all mRNAs contain several sequences that are potentially recognised by microRNAs, however, only a small part of them have been proven to function for regulation [6]. The base pairing between each microRNA and the selected mRNA occurs in the miRNA-induced silencing complex (miRISC) and leads to decreased protein expression through translational repression and/or mRNA degradation. Briefly, when the microRNA is loaded onto an Argonaute protein, the passenger strand is released and the guide strand hybridises to the complementary region on the mRNA [3]. TNR6A (also known as GW182), binds to the Argonaute complex and helps to recruit the CCR4-NOT and PAN2-PAN3 deadenylase complexes as well as the decapping complexes [7,8]. Several studies report that microRNAs repress mRNA translation, however it has been suggested that the main mechanism of regulation is mRNA decay [8]. This is supported by the observation that the majority of the studies which demonstrate microRNA-mRNA interaction by luciferase assay [9], also observe an mRNA decreasing upon miRNA mimics transfection. Several RNA binding proteins also contribute to the outcome of the regulation process, by recognising specific sequences or secondary structures on the 3’ UTR of the mRNAs [10]. The binding may mask the seed sequence or may facilitate the microRNA binding. This partially explains the differences between the set of predicted versus real targets.

A growing body of evidence points to the key role of microRNAs in gene regulation, and some of them control multiple genes that are altered in pathological conditions, so that they are under investigation to develop new therapeutic strategies. Therefore, it is crucial to have the most reliable knowledge on the network of the targeted genes, to be able to foresee most of the effects. To this aim, several repositories collect microRNA interaction data, and considerable efforts have been made to standardise annotation vocabularies and RNA identificatives to support data integration [11,12,13]. Nonetheless there is still a significant discrepancy, among different databases, as to which experimental approaches demonstrate “true” interactions, especially for high-throughput assays. Immunoprecipitation of AGO2 in Clip-seq experiments [14] return lists of microRNAs and mRNAs co-immunoprecipitated, and some resources capture these results as interaction couples, while others recognise their value to improve the performance of prediction tools but not as proof of interaction [13]. An adaptation of this technique was proposed to specifically investigate microRNA targets by linking the ends of RNA duplexes before sequencing [15]. In the last years, hundreds of papers have been published, showing the physical binding of a microRNA to its partner by luciferase assay performed on the wild type and on the mutated mRNA. This approach is among the few that are generally accepted to demonstrate the binding between the two RNAs, together with the RNA-Ip which is more rarely utilized [13,16,17]. The UCL functional annotation group [16] and the IntAct database [17] have annotated a collection of such interactions, although the coverage is still limited. MirTarBase [18] and RNAInter [19] have a higher coverage, and it is possible to select the stronger evidence from the other annotations. To investigate the main features of the microRNA mediated regulation, we have built a network with the most reliable interactions and analysed the main characteristics of the post-transcriptional regulatory circuit.

2. Materials and Methods

2.1. Dataset Integration

Human microRNA-mRNA interactions were downloaded from four repositories in Jan 2024: miRTarBase [18] RNAinter [19] QuickGO [16,20], IntAct [17]. Interactions from QuickGO and IntAct are consistent to the International Molecular Exchange consortium common standards (https://www.imexconsortium.org/) [21], while interactions derived from miRTarbase and RNAinter were filtered for low-throughput, strong evidence. More precisely, to assure the physical binding of the microRNA to the mRNA, we selected “luciferase assays” or RNA-IP, and disregarded “western blot” and “qpcr” as these can be the result of a direct or indirect effect (Figure1,[13]) To integrate the data, the Uniprot ID mapping tool (https://www.uniprot.org/id-mapping) [22] was used when necessary. Each interaction was considered once and all duplications were deleted. Cytoscape tool [23] was then used to visualise and analyse the network. The list of DNA-binding transcription factors was retrieved from [24]. To obtain an estimation of the number of genes regulated by each transcription factor we downloaded gene regulation data (TF-TG) from TFlink “small scale dataset” filtering out inferred interactions [25]. This should contain only interactions verified in low throughput experiments. TF regulating zero genes according to this list were deleted.

2.2. Network Analysis

The microRNA-mRNA network was constructed and the number of edges for each node (defined as node degree) was calculated [23]. Degree classes (dc) were calculated based on how many times a certain degree appears (abundance). As an example, if there are two nodes with 94 interactions each one, then degree 94 belongs to the degree class 2. Three datasets were produced, listing the abundance of dc relative to mRNA, and miRNA and for the subclass of transcription factors (see below). To compare the abundance distribution of the datasets, the regression model best fitting each abundance distribution was computed on the basis of several models proposed in the R library ‘sads’ [26]. The best model was selected on AIC basis, and where the abundance distributions converged on the same model, then this was taken as an evidence of similarity. We compared the abundance distribution of transcription factors (TF) with that of all mRNAs. Transcription factors regulated by zero microRNAs were deleted and the others were subtracted from the mRNA dataset, and the two curves were compared by means of the AIC value.

Mutual Information (MI) analysis was applied for the effective estimation of the possible relationship between the number of microRNAs regulating each TF and the number of genes regulated by the same TF, by means of the R library “entropy”[27]. MI can be applied to assess the presence of an unspecified relationship among data sets, when a linear model is not applicable. It measures the amount of information obtained about one variable, when measuring the other [28,29]. MI returns an absolute value, which can be more easily interpreted if expressed as relative MI (rMI), i.e. the ratio between MI and the maximum MI that can be extracted for the dataset [30]. To evaluate confidence intervals (CI) and significance (p), MI and rMI were bootstrapped by means of the R libraries 'simpleboot', 'boot' and 'boot.pval'.

2.3. Gene Ontology and Pathway Enrichment Analysis

Genes interacting with more than 20 microRNAs were analysed for Gene Ontology terms enrichment against the whole genome. To this aim, Gprofiler [31], VLAD [32] (moved to https://www.informatics.jax.org/vlad/) and BinGO [33] were utilised, giving very similar results. The data presented in Figure 3 show the enrichment for molecular function and are obtained with BinGO, the significance was imposed below 0,0005 and the whole genome was considered as benchmarking. Supplementary figure 1 was drawn with VLAD. Pathway enrichment analyses were conducted with Reactome analysis tool [34].

2.4. Other Statistical Analysis

Intersections of common target genes (Figure 5) were calculated with https://bioinformatics.psb.ugent.be/webtools/Venn/.

To correlate 3’ UTR lengths with predicted and verified interactions, data tables were generated as follows: transcript lengths and 3’ UTR lengths were extracted from ensembl (https://www.ensembl.org/index.html) database and refers to the “main” transcript as defined by GIFT tool (https://www.ebi.ac.uk/gifts/). The number of interactions was extracted from the integrated dataset described above. Predictions were calculated with miRDB [35] and miRWalk [36]. miRwalk prediction was calculated on the same Refseq as for the UTR and length measures and filtered at 0.95. Pearson correlation coefficient was calculated using the ggpairs functions from GGally package [37] in R [38].

3. Results

3.1. General Characteristics of the microRNA-mRNA S

Although any gene is a potential target for microRNAs, as suggested from the numerous target prediction methods, low throughput experiments seem to reveal a significant discrepancy among transcripts. In a recent paper, the collection from the literature of microRNAs which regulate genes involved in rare diseases has shown that for the majority of mRNAs no microRNA binding has been demonstrated, while a small number of mRNAs is heavily regulated [13]. In order to gain more insights into the distribution of microRNA targets, we have integrated 4 manually curated datasets [13,16,18,19] and filtered them to only provide strong-evidence, low-throughput results, as described in detail in Figure 1. We used the resulting list of “gold standard” interactions to investigate into post-transcriptional regulatory mechanisms.

According to this dataset, 3842 human genes are regulated by at least one microRNA, which corresponds to 20% of the estimated 19400 human coding genes (www.genenames.org). We calculated the number of edges for each node (defined as node degree, dc). Concerning the mRNA nodes, the smallest class, with dc=1, contains 1964 presences, i.e. more than 1900 mRNAs interact with one microRNA; while the largest class with dc=108 is present once (and corresponds to PTEN phosphatase). The abundance decreases quickly between 5 to 10 degrees, so that a minority of nodes have more than 10 interactors (Figure 2A, Supplementary Table 1). This scenario is consistent with the abundance distribution best fitting a power law, where large values are exponentially less likely, and a large proportion of occurrences is held by a small number of entities (Pareto distribution). Only 2% of the 3842 mRNAs regulated by microRNAs are controlled by more than 20 microRNAs, and we defined them as network hubs for the following analyses.

The miRNA dataset lists 962 microRNA nodes, which have an average higher degree in comparison to mRNAs. The smallest class with dc=1 shows 239 presences and 14% of interacting microRNAs regulate more than 20 genes. Notably, thirteen microRNAs target more than 100 genes each (Figure 2B, Supplementary Table 1), and this number is destined to increase as new data will be published.

3.2. mRNA Highly Regulated Are Regulators

To gain more insights into microRNA-mediated post-transcriptional regulation, we investigated the molecular functions of the highly regulated genes, to estimate if specific typologies of genes have a higher propensity to be controlled by multiple miRNAs. The Gene Ontology is the most used bioinformatic resource to unify and interpret gene functions and attributes for biomedical purposes [39]. We calculated the enrichment in genes performing specific functions among those regulated by more than 20 microRNAs, in respect to the whole set of human genes, using three different tools: Gprofiler [31] VLAD [32] and BinGO [33]. These three methods gave similar results, showing a very high enrichment for genes related to transcription regulatory activity. Figure 3 shows Ontology Terms related to the transcription and transcription regulation, highly enriched among “hub” genes in respect to the others. In Supplementary table 2 all the enriched terms retrieved with BinGO are listed with their frequency in hub genes and in the whole genome. Reactome analysis tool [34] also returns among the top-ranking pathways “Generic Transcription Pathway” and “RNA Polymerase II Transcription” and “Gene Expression”, all with FDR 1,05E-14.

A similar result had been obtained in [40] from the analysis of hub genes predicted with Targetscan [41].

Figure 3. Transcription regulation function in genes controlled by more than 20 microRNAs.

Figure 3. Transcription regulation function in genes controlled by more than 20 microRNAs.

A) The histogram presents the molecular-function GO terms statistically overrepresented (significance level 0,0005) in the hub genes compared to the whole genome, clearly showing an enrichment for genes related to transcription regulation and signal transduction. The complete list of molecular functions enriched in hub genes and their frequencies are listed in Supplementary Table 2. B) Ontology graph drawn with BinGO in Cytoscape [33], showing the GO hierarchy of enriched child terms related to transcription. The area of a node is proportional to the number of genes annotated to the corresponding GO category among the “hubs” and the colours indicate the statistical significance, higher in dark orange nodes.

On the other hand, not all the transcription factors are regulated by microRNAs. In a recent paper, members of the GO Consortium have produced a curated catalogue of 1456 human DNA binding transcription factors [24]. We compared the distribution of the interacting microRNAs of these 1456 genes with the whole set of human mRNAs, and we found a very similar distribution (Figure 4). For 72% of the DNA-binding TF listed in [24] no interacting microRNA has been reported, while 28% have at least one. Only 25 transcription factors are regulated by more than 20 microRNA, corresponding to 1.7%. The rank abundance distribution of microRNAs binding to transcription factors fits a Pareto distribution as for the other mRNAs, with a minimum value of AIC = 174 as compared to AIC= 401 for the other genes, indicating that they are represented by a power law distribution.

To better estimate if there is any relationship between the microRNA regulation and the transcriptional regulation we performed a Mutual Information analysis, considering, for each transcription factor, two variables: the number of interacting microRNAs and the number of genes regulated by. To estimate the number of genes regulated by each transcription factor we referred to TFlink [25]. The MI value (MI=0.1161408) and the rMI value (rMI=0.02567784), both fall in a range within the CIs with a highly significant p-value, indicating a non-random co-occurrence. However, the low rMI value indicates that other variables influence the numbers. From these observations we suggest that, although genes controlled by multiple microRNAs are highly enriched for transcription regulation functions, not all the transcription factors are regulated by microRNAs. As shown in Supplementary Table 3 however, the most interconnected factors are highly controlled by microRNAs.

Genes involved in receptor signalling and kinase binding are also enriched among microRNA regulated genes (Figure 3A), and the proteins they code for are highly interconnected in PPI networks. For example, considering the term “protein kinase activity”, the genes annotated for this term code for the proteins gsk3b cdkn1a cdk6 ccnd1 erbb2 akt1 met tgfbr1 egfr igf1r tgfbr2, all interacting with hundreds of proteins, according to IntAct database [17]. Interestingly, the mRNA which interacts with the highest number of microRNAs, more than one hundred, encodes for the Phosphatase and Tensin Homolog PTEN , a well-known tumour suppressor gene, which inhibits the PI3K/AKT growth signalling. PTEN is downregulated by several “oncogenic” microRNAs (see figure 5) which exert their tumorigenic function through this key phosphatase [42].

Supplementary Figure 1 shows the enrichment diagram obtained by VLAD [32], where p values for enriched terms are reported.

3.3. microRNAs Regulate Multiple Targets

It has been proposed that one of the main functions of microRNAs is to modulate gene expression to finely coordinate the level of proteins involved in the same biological process or pathway. According to the dataset presented in this work, 590 microRNAs control more than 2 targets and 13, listed in Table 1, regulate more than 100 transcripts each (Supplementary Table 1). Many microRNAs are associated with cancer or other pathologies and their expression is deregulated in the affected cells [43]. Table 1 shows that miRNA controlling one hundred or more of genes are “oncomiR” or “tumour suppressor”, which means that they promote or repress cancer phenotype respectively. Despite their opposite effects, these microRNAs have some common targets: the Venn diagram in Figure 5 shows that 23 targets are in common between miR-155 and miR-21 (both described as oncomiR) and 12 genes are regulated by miR-155 (oncomiR) and miR-34a (tumour suppressor).

Figure 5. OncomiR and suppressor microRNAs share common targets. A) Venn diagram showing the overlapping genes regulated by two oncomiR: miR-155 and miR-21. The 23 common targets are: dock1 pten myd88 nef pik3r1 apaf1 egfr msh6 socs6 smarca4 bcl2 bcl6 thrb cebpb bcl10 icam1 jun foxo1 socs1 vhl msh2 casp3 sox6. B) Venn diagram showing the overlapping genes regulated by the tumour suppressor miR-34a and the oncomiR miR-155. The 12 common targets are: xiap, nfe2l2, myb, spi1, rad51, csf1r, bcl2, agtr1, cebpb, bcl6, myc, ccnd1.

Figure 5. OncomiR and suppressor microRNAs share common targets. A) Venn diagram showing the overlapping genes regulated by two oncomiR: miR-155 and miR-21. The 23 common targets are: dock1 pten myd88 nef pik3r1 apaf1 egfr msh6 socs6 smarca4 bcl2 bcl6 thrb cebpb bcl10 icam1 jun foxo1 socs1 vhl msh2 casp3 sox6. B) Venn diagram showing the overlapping genes regulated by the tumour suppressor miR-34a and the oncomiR miR-155. The 12 common targets are: xiap, nfe2l2, myb, spi1, rad51, csf1r, bcl2, agtr1, cebpb, bcl6, myc, ccnd1.

Interestingly, miR-155-5p is the most connected microRNA in the network, and it regulates approximately 260 genes. It is considered an oncomiR as it was found overexpressed in lymphoma, leukaemia and several solid tumours. In normal tissues, it is expressed at a very low level with the exception of hematopoietic cells, lymph nodes and spleen, which suggest its implication in inflammatory response. Indeed, pathway analysis tool on its targets shows that the list of regulated genes is enriched in proteins involved in: Interleukin-4 and Interleukin-13 signalling (15 genes), signalling by Interleukins (33 proteins), Cytokine Signalling in Immune system (39 proteins) and Immune System (68 proteins) (FDR 1.54E-11 , 3.51E-10, 1.66E-7 and 5,5 x E-4 respectively). It is also significantly enriched in both TLR3 and TLR4 cascade with 11 and 13 proteins respectively.

It has been suggested that genes regulated by the same microRNA may be involved in the same pathway to orchestrate the coregulation [48,49]. Consistently, we observed that several microRNAs controlling multiple genes, target more entities of a specific pathway. For example, miR206 regulates 6 genes (out of 15) involved in the Reactome “NFE2L2 regulates pentose phosphate pathway “ (R-HSA-9818028, FDR 6.18E-8), or miR 33a regulates 7 genes involved in Regulation of lipid metabolism by PPARalpha (R-HSA-400206, FDR 5.3E-9) as shown in Figure 6A. Interestingly, miR-105-5p downregulates 3 genes involved in the pathway “Regulation of lipid metabolism by PPARalpha”, while miR-105-3p, transcribed from the complementary strand, downregulates the hub PTEN, involved in the same pathway (Figure 6B).

3.4. The Length of the Untranslated Regions Do Not Correlate with the Abundance of Interactors

Most human protein-coding genes contain regions complementary to seed sequences of the microRNAs, that potentially mediate their binding. The seeds are short nucleotide sequences (approximately 6-7 nucleotides), so the complementary stretches may occur by chance in the transcripts without mirroring any regulative function. For example, 6 nucleotides may occur accidentally every 4096 nucleotides, and the longer is the transcript sequence, the more likely it will contain potential microRNA binders. The number of potential binding sites is even larger due to the fact that the base pairing rarely needs a perfect match in animals but few mismatching are also allowed [6,7]. Interestingly, the average length of UTR has been correlated with species complexity, even within vertebrates [50,51]. Moreover, it has been demonstrated that different transcript isoforms, bearing longer/shorter UTRs, may be differentially regulated by microRNAs [52,53].

We performed a simple test, comparing the 3’ UTR length of a sample of transcripts with the number of interactors experimentally verified and the number of predicted binding sites, to verify if more regulated mRNAs bear longer untranslated regions. The sample included highly connected mRNAs as well as transcripts with no verified interactors (Supplementary Table 4). To estimate the predicted binding sites, we selected two high-performance prediction tools: MiRDB [54] and MirWalk [55], both using additional criteria besides the matching sequence, spanning from deep learning from high-throughput data, or expression profiling data, to energetic considerations.

As expected, the two predictors returned some common as well as some different microRNAs and in both cases we found a significant positive correlation between the UTR lengths and the binding sites predicted (r=0.654 for miRDB and 0.737 for MirWalk). On the contrary, when comparing the length of the UTR regions with the number of microRNAs that were demonstrated to bind them, we found very low correlation (r=0.312) suggesting that the true binding sites are selected and maintained and do not depend on the length of the sequence. Surprisingly, several of the sequences predicted to have hundreds of potential interactors, have not been demonstrated for any of them, generating an unexpected discrepancy between the number of predicted versus verified binders (Supplementary Table 4).

4. Discussion

Since the relevance of the microRNA-mediated, post-transcriptional regulation is supported by evidence from thousands of studies, in normal cells and pathologies, a paradox must be fixed, as the amount of protein or mRNA that decreases in response to the expression of a microRNA is quite limited, and the rate of repression is usually in the range of 2-3 folds, making it difficult to explain the phenotypes observed in microRNA knockdowns [10,56]. This points to the possibility that microRNAs exert their most evident functions regulating key genes such as master transcription factors, and that they may act on more genes involved in the same process or pathway, thus potentiating the final effect. This is consistent with previous observations, obtained on a predicted microRNA-mRNA network, where extensive interactions between microRNAs and transcription factors were highlighted [40]. Notably, feed-forward loops were also shown, with transcription factors regulating microRNAs which downregulate them, or both co-regulating a third protein [57]. For example, Myc activates the expression of the the transcription factor E2F1 and of the miR-17 cluster, composed of six miRNAs, comprising miR-17, miR18a, miR19a, miR20, miR19b, and miR92 [58], which, in turn, control E2F1 level [59]. Genome-wide profiling of the transcriptome will help in elucidating which regulative proteins control microRNA expression, better clarifying the extension of the cross-talk [60]

In this study, the analysis of a collection of trustworthy targets have confirmed that microRNAs finely control driver transcription factors and other key regulatory proteins. In fact, we show a remarkable enrichment of transcription factors among highly regulated mRNAs, yet limited to a defined subgroup of genes, possibly those determining the cell fate. For example, Zeb, Twist1 and Snai1 have been recognised as EMT drivers [61], and all of them are controlled by more than 20 miRNAs (Supplementary Table 3).

Regulatory circuitries are commonly tampered in cancer cells and in other diseases, so they are likely to represent optimal therapeutic targets. Indeed, many cancer associated genetic mutations bear down on genes and DNA regions involved in the transcriptional and post-transcriptional control, as transcription factors, signalling proteins, chromatin modifying enzymes or ncRNAs. Inhibitors of oncogenic transcription factors and oncomiRs have recently entered in clinical trials [45], so that elucidating the interplay between different levels of regulation is crucial to predict the effect of interfering drugs and molecules. The presence of multiple overlapping targets for microRNAs candidated for therapies (Figure 5, Table 2), suggests the urgency of an in-depth understanding of the full picture of relationships.

The purpose to complete a trustworthy microRNA-mRNA interaction network is still far to be achieved: the intersection between low scale luciferase assay network (this study) and the network obtained with a high throughput CLASH approach [16], corresponds to less than 1% of the interactions (data not shown), pointing to incompleteness on both sides, similarly to the first attempts to deciphering PPI networks [62]

On the other hand, information coming from high-throughput experiments and predictions weaken the reliability of the result interpretation [63]. A similar problem, concerning the reliability of protein-protein interactions (PPI), was handled in the last 20 years by the IMEx consortium [21], by defining guidelines as what to consider experimental proof of binding and by adding a score of reliability to each interaction [64], which consider, among the others, how many times it has been described in the literature. In the IntAct database, the same score was applied to evaluate microRNA interactions, but the amount of data concerning microRNAs is still too limited to allow for filtering [13]. Several PPI databases follow common principles in order to offer highly reliable data to the user [21]. On the contrary, ncRNAs are mainly collected in independent resources, following different criteria not always clear to the user [65]. Moreover, the integration of gene regulatory networks (GRN) with post-transcriptional regulation is still at its infancy [66]. As the relevance of this intricate and dynamic interplay which determines the cell fate is clearly emerging, an in-depth understanding of the global and local architecture of the regulatory system is certainly needed.

5. Conclusions

Countless studies have contributed to the elucidation of the effects of microRNAs on specific genes, but very little is still certain about the extent of regulation mediated by microRNAs and how they cooperate with the other layers of regulation. Transcriptional and post-transcriptional effects certainly contribute to the final expression level, but the knowledge in the two fields is still modestly interconnected. Moreover, the excess of potential binding sites, largely exceeding the microRNA abundance, imply that the regulation may be sensitive to microRNA levels [7], which also should be taken into account when drawing models.

Despite the lacking knowledge on the biology of ncRNAs, from the regulation of their expression to their function, and even on their coherent classification, there is no doubt as to these molecules have made a revolution in our understanding of the cell behaviour and must, from now on, be taken in consideration on any perturbation prediction.