Combinatorial network of primary and secondary microRNA-driven regulatory mechanisms.

Recent miRNA transfection experiments show strong evidence that miRNAs influence not only their target but also non-target genes; the precise mechanism of the extended regulatory effects of miRNAs remains to be elucidated. A hypothetical two-layer regulatory network in which transcription factors (TFs) function as important mediators of miRNA-initiated regulatory effects was envisioned, and a comprehensive strategy was developed to map such miRNA-centered regulatory cascades. Given gene expression profiles after miRNA-perturbation, along with putative miRNA-gene and TF-gene regulatory relationships, highly likely degraded targets were fetched by a non-parametric statistical test; miRNA-regulated TFs and their downstream targets were mined out through linear regression modeling. When applied to 53 expression datasets, this strategy discovered combinatorial regulatory networks centered around 19 miRNAs. A tumor-related regulatory network was diagrammed as an example, with the important tumor-related regulators TP53 and MYC playing hub connector roles. A web server is provided for query and analysis of all reported data in this article. Our results reinforce the growing awareness that non-coding RNAs may play key roles in the transcription regulatory network. Our strategy could be applied to reveal conditional regulatory pathways in many more cellular contexts.

INTRODUCTION

As the major endogenous triggers for posttranscriptional silencing, microRNAs (miRNAs) are predicted to target over one-third of human genes and regulate a variety of biological processes, including development, cell death, cell proliferation, hematopoiesis, and nervous system patterning (1). Complementarity between the 5′ end of miRNA (positions 2–7) and potential targets, also known as seed-region matching, is critical for miRNA-based regulation, but perfect seed pairing alone is insufficient to predict functional targets (2). Computational algorithms focused primarily on sequence complementarity have been devised to predict miRNA targets; such algorithms include PicTar (3), TargetScan (4,5) and miRanda (6). Because these algorithms are thought to generate a significant proportion of false positives, independent experimental validation is required for refined miRNA target screening.

In general, miRNA silencing of target gene expression is achieved by mRNA degradation or translational inhibition (1). Recently, to help identify targets of miRNA-induced degradation (‘degraded targets’), miRNA-perturbed gene expression (MPGE) experiments have been performed, in which a particular miRNA of interest is over-expressed or knocked-down, and then large-scale mRNA expression levels are measured. (Note: In this article, we will use the term ‘MPGE dataset’ to denote datasets generated by over-expression of miRNA, as all datasets involved in this study are of this type.) Most typically, a large set of putative miRNA targets is first identified using complementarity-based target-prediction algorithms; from the putative targets, actual degraded targets are confirmed based on significant down-regulation of mRNA levels in the MPGE dataset (7–9). This experimental design, however, can be error-prone, due to the complex nature of miRNA-triggered gene regulation. In a MPGE dataset, only a small fraction of the down-regulated genes are direct miRNA targets; non-target genes may also have decreased mRNA levels because of a repressed activator. In addition, mRNA levels of both target and non-target genes may be up-regulated because of a repressed inhibitor. Such situations were indeed discovered in our investigation of some MPGE datasets (Figure 1A–C). In general, we propose that the mRNA level changes observed in MPGE datasets are combinatory effects resulting from both direct and indirect regulation triggered by miRNAs; thus, to resolve miRNA-induced gene regulation, a combinatorial network analysis is needed.

Figure 1.

(A) mRNA level changes of miR-124’s putative target genes, measured at various time points after transfection of miR-124 into the HepG2 cell line (GDS2657). The putative targets are a combined set of prediction results from three predicting algorithms: PicTar4way, TargetScanS and miRandaXL. The x-axis indicates time after transfection in hours, and the y-axis indicates the log2-transformed ratio of gene expression between treatment (after transfection) and control (before transfection). (B) mRNA level changes of miR-124’s putative non-target genes in GDS2657. Putative non-target genes are the whole set of genes found in GDS2657 minus the putative target genes. (C) Numbers of up- (up panel) and down-regulated (below panel) genes at different time points in GDS2657. The up- or down-regulated genes are further divided into three groups based on our target-identification work: direct targets (first targets), secondary targets (second targets), and the rest. The exact numbers are included in Supplementary Table 8.

Given that miRNAs exert direct as well as indirect regulatory effects, an important question arises: among direct miRNA targets, which serve as important mediators of the regulatory signal? Previous works have found that transcription factors (TFs) prevail among miRNA targets in plants (10) and insects (6); we also discovered an enrichment of TFs among the potential targets of human miRNAs (Supplementary Figure 1), as well as a preference for TFs with more targets over those with fewer targets (one-sided Wilcox test P < 0.05). Given these observations, we envisioned that miRNAs trigger a two-layer regulatory network in which TFs play an important role in propagating miRNA-initiated regulatory effects. In this kind of network, mRNAs directly regulated by miRNAs via degradation or translational inhibition are at the primary level of the regulatory cascade, whereas target genes of miRNA-regulated TFs are secondary targets, subject to indirect regulation by miRNAs (Figure 2A). This two-layer regulatory mechanism explains how non-target mRNAs can also be influenced by miRNA, as well as how miRNA targets can be up-regulated in MPGE datasets.

Figure 2.

(A) A hypothesized two-layer regulatory mechanism triggered by miRNAs. miRNAs regulate their primary targets, including TFs, via degradation or translational inhibition. Regulated TFs propagate miRNA-initiated regulation to secondary targets, causing changes in their mRNA levels. Because TFs may regulate some primary miRNA targets, some primary targets may be subject to secondary regulation, as well. (B) Flowchart illustrating our approach. Putative miRNA–target relationships predicted by classical sequence-based algorithms and over-expression datasets were analyzed in a hypothesis test model to identify degraded targets. Putative miRNA–target relationships, over-expression datasets and curated TF–target relationships were analyzed by stepwise linear regression, to identify miRNA-regulated TFs. The two outputs (degraded targets and regulated TFs), along with TF targets, were combined to map miRNA-triggered two-layer regulatory networks, upon which function analyses (GO term and KEGG pathway enrichment) were performed to identify biological themes associated with each network.

In this work, we developed a comprehensive strategy to characterize the hypothesized two-layer regulatory network by integrating MPGE datasets with miRNA–gene and TF–gene regulatory relationships (Figure 2B). Given the expression patterns revealed in the MPGE datasets, highly likely degraded targets were extracted from the set of all putative miRNA targets by a non-parametric statistical test; TF mediators of miRNA-triggered regulation and their downstream targets were mined out through linear regression modeling. The identified primary and secondary miRNA targets were mapped to miRNA-centered regulatory cascades, to sketch out the combinatorial regulatory program controlled by miRNAs and TF mediators.

MATERIALS AND METHODS

miRNA–target relationships in miRGen

Candidate miRNA–target relationships were downloaded from miRGen 3.0 (11) (http://www.diana.pcbi.upenn.edu/miRGen), which collects miRNA target results from widely used target-prediction programs. We obtained the union of prediction results from three algorithms: PicTar (3), TargetScan (4) and miRanda (6); these results included 118 408 relationships among 276 human miRNAs and 10 255 targets (for detailed statistics, please see Supplementary Table 2).

TF–target relationships developed from UCSC and TRED data

Two files, TFbsConFactors.txt and TFbsConsSites.txt, were downloaded from UCSC hg18 (http://genome.ucsc.edu/). TFbsConsSites gives predicted chromosomal coordinates of TF binding sites (TFBSs) on human, mouse and rat genes, while TFbsConfactors.txt links the internal TF accessions to SWISS-PROT IDs. These SWISS-PROT IDs were further converted to NCBI gene IDs via BioMart (http://www.ebi.ac.uk/biomart/), and NCBI’s homologene.data file was used to find the human homologs of mouse and rat TFs, enabling us to compile an enlarged set of human TF–TFBS relationships.

Also from UCSC hg18, we downloaded gene coordinate information (refGene.txt file), which specifies the chromosomal locations of 18 620 human genes. The promoter region of each gene [from 1 kb upstream of the transcription start site (TSS) to 0.5 kb downstream of the TSS] was scanned for the TFBSs identified in the above TF–TFBS relationships. If an occurrence of a certain TFBS was found, the corresponding TF was linked with that gene. In this way, we developed a set of TF–gene regulatory relationships.

In addition, we retrieved TF target information from another source, the TRED database (12) (http://rulai.cshl.edu/TRED/), which collects mammalian cis- and trans-regulatory elements, accompanied by experimental evidence. Because the TF target information contained in TRED was of better quality and was not completely covered by the above results, we supplemented the above set of TF–gene regulatory relationships with that retrieved from TRED. The final TF–gene set included 130 338 relationships between 214 human TFs and 16 534 targets. For more information, see Supplementary Table 3.

MPGE datasets

Five groups of MPGE datasets (GDS1858, GDS2657, GSE6474, GSE6838 and GSE7864) were downloaded from the GEO database; these groups include 53 individual datasets involving 19 miRNAs. The GDS1858 dataset group (2) includes data on HeLa cells at 12 or 24 h after transfection with wild type or mutant miR-1, miR-124, or miR-373. GDS2657 (13) includes gene expression profiles at seven time points (4, 8, 16, 24, 32, 72 and 120 h) after overexpression of miR-124. GSE6474 (14) includes four replicated measurements of gene expression changes after overexpression of let-7a. GSE6838 (15) includes gene expression data over a time course (6, 10, 14 and 24 h) after overexpression of 12 different miRNAs from the miR-16 family. GSE7864 (16) reports gene expression changes at 24 h after overexpression of miR-34a, miR-34b and miR-34c, in five different cell lines. For more information, see Supplementary Table 4.

In each MPGE dataset, every gene is associated with a fold-change value in mRNA level triggered by miRNA perturbation; this value was log2 transformed and is referred to henceforth as mRNA change (MC).

miRNA families

Information on miRNA families was downloaded from miRBase (http://miRNA.sanger.ac.uk). All 19 miRNAs involved in this study were assigned to one of nine miRNA families (see Supplementary Table 5).

Determining miRNA’s degradation-inducing ability and its degraded targets

Classical miRNA target-prediction algorithms have predicted miRNA–gene regulatory relationships exclusively on the basis of seed-pairing and inter-species conservation. While the putative targets include a high percentage of false positives, non-target predictions are much more likely to be correct. Indeed, the MCs of putative non-targets have no evident bias generally: up- and down-regulations are observed almost equally (Figure 1B). Their MC distribution, caused mainly by secondary regulation by miRNA, was termed D0. By contrast, potential targets comprise a sizeable portion of false-positive predictions, and true targets therein are further divided into degraded targets and translationally inhibited ones. As the degraded targets tend to be down-regulated and their MCs have a mean value that is significantly lower than that of the non-targets, the distribution of their MCs was termed D1. Like non-targets, false-positive predictions and translationally inhibited targets are subject to potential secondary regulation by miRNA, with an MC expected to follow D0. In summary, non-target MCs were expected to follow D0; MCs of putative targets were expected to follow the combinatory distribution D2 = D0 + D1. The degradation-inducing ability of miRNA was judged by comparing D0 and D2, with a one-sided Kolmogorov–Smirnov (K–S) test.

To answer whether a putative miRNA target, gene i, is a truely degraded targets, we performed a non-parametric statistical test on a null hypothesis that MC is sampled from D0. As we argued above, MCs of translationally inhibited targets and false-positive predictions should follow the same distribution as those of putative non-targets, so the MCs of putative non-targets were used to approximate the null distribution D0, and the significance of gene i being a true degraded targets (P-value) was calculated as the proportion of genes in non-targets that have a lower MC than gene i. The null hypothesis and alternative hypothesis are as follows: The P-values out of this non-parametric statistical test were then adjusted to false discovery rates (FDRs) using the BH method (17). With a FDR threshold, α1, we picked out the degraded targets from all putative targets.

H0 : MC follows D0, i.e. target i is a translationally inhibited target or a falsely predicted target.

H1 : MC does not follow D0, but follows D1, i.e. target i is a degraded target.

Identifying TF mediators that propagate miRNA-initiated regulation to secondary targets

It is assumed that, if a TF is repressed by a miRNA (through either mRNA degradation or translational inhibition), the regulatory effect would be systematically transmitted to the TF’s target genes (i.e. miRNA’s secondary targets). Thus, by comparing the MCs of a TF’s targets and those of its non-targets, we could infer whether a particular TF was inhibited by the miRNA or not.

Given the miRNA–gene and TF–gene regulatory relationships developed as described above (see ‘Materials and Methods’ section), we assumed that a gene's MC could be regarded as a linear combination of regulatory effects from miRNA and TFs, which could be estimated using conventional linear regression. Theoretically, all TFs identified as a miRNA’s putative targets should be included in the linear combination; in our situation, however, involving too many candidate TFs would render intolerable time-complexity. Therefore, a feature-selection procedure was performed first to evaluate the power of each candidate TF in explaining a gene's MC. The MC of a gene (g) was preliminarily modeled as a combination of the primary regulatory effect from the miRNA and one secondary regulatory effect, contributed by a TF, as follows: Here, MC denotes g’s MC; m refers to the particular miRNA perturbed in the MPGE experiment; a quantifies this miRNA’s primary regulatory strength; TF denotes a TF putatively regulated by m (i.e. m-TF appeared in the miRNA–gene regulatory relationship set); a quantifies this TF’s regulatory strength; b and b, with values zero or one, encode the putative regulatory relationships between g and miRNA m, and g and TF, respectively; interception is a constant across all genes capturing a baseline MC; and err is the residual error, assumed to follow normal distribution. In Equation (1), MC, b, and b are known, and the unknown coefficients a and a need to be estimated.

All TFs predicted to be regulated by miRNA m underwent the above feature-selection procedure, and those with P-values smaller than α2 were further filtered by Equation (2), the formal linear regression model of g’s MC. MC, b and b are exactly the same as in Equation (1), whereas regulatory strengths of the miRNA and TFs (i.e. a and a) need to be re-estimated.

We used stepwise linear regression to optimize Equation (2) according to BIC score (18), ultimately yielding a model best interpreting the variation in gene MCs. Significant TFs surviving in the final model were regarded as TF mediators that delivered miRNA-triggered regulatory effects to secondary targets. In other words, if the coefficient for TF (i.e. a) turned out to be significantly different from zero, it suggested a significant difference between MCs of TF’s target genes and those of its non-target genes; hence, TF was inferred to be inhibited by the miRNA of interest. If TF itself did not show a significantly down-regulated MC, its inhibition was very likely achieved by translational inhibition rather than mRNA degradation; this might be worthy of validation with carefully designed experiments.

The fitness of the regressed model to the observed data was evaluated by BIC score during standard implementation of stepwise linear regression. In our situation, however, we were afraid that the significance of the model could be attributed largely to the effect of miRNA-induced degradation, without adequately reflecting TF-mediated secondary regulatory effects. To rule out this possibility, we additionally shuffled the TF–gene regulatory relationships randomly 100 times and defined a ‘shuffling P-value’ as the proportion of regressed models in the 100 TF-target shuffling experiments that have an identical or greater BIC score than that of the regressed Equation (2). This additional shuffling-based evaluation was designed to examine particularly the contribution of the TFs in the final model; if the shuffling P-value is very low, that is, an equally valid model is rarely observed in the random shuffling experiments, we could safely conclude that the TFs included in the model do make a significant contribution in explaining the variation of gene MCs.

Building miRNA-centered two-layer regulatory networks

Combining primary and secondary miRNA targets identified as described above, we mapped combinatorial regulatory networks in which miRNAs act as triggering regulators. To achieve a concise and refined representation of these networks, we reserved only differentially expressed genes with MC values in the upper 5% or lower 5% ranges of the overall MC distribution.

In addition, we performed gene set enrichment analyses with GO terms and KEGG pathways, so as to summarize the biological themes of each miRNA-triggered regulatory network. For the genes included in a regulatory network, we counted the genes annotated to a GO term or a KEGG pathway, and performed a one-sided Fisher's exact test to compare the percentages of annotated genes inside the network and outside the network. A BH multiple test correction (17) was applied on the P-values to estimate the corresponding FDRs.

The methods proposed above to identify primary and secondary targets of miRNA-initiated regulation and to build miRNA-centered regulatory networks have been implemented as a web application at http://www.biosino.org/de.cleaved.repressed. Researchers can use this application to build miRNA-centered two-layer regulatory networks with their own MPGE datasets.

RESULTS

Predicted miRNA targets of high reliability

Setting the FDR threshold α1 = 0.25 for degradation target identification and α2 = 0.05 for pre-screening of TF mediators, we implemented the approaches described above on five public dataset groups encompassing 53 MPGE datasets and 19 distinct miRNAs. Results associated with individual MPGE datasets are included in the ‘wrapped results’ available at http://www.biosino.org/~kanghu/DCR/ supplementary file1.zip, and full lists of degraded targets are available in Supplementary Table 6. Predicted direct targets are summarized for each MPGE dataset in Table 1, and for each miRNA in Supplementary Table 7.

Table 1.

Targets of miRNA-induced degradation and TF mediators of miRNA-triggered regulation, summarized for each of 53 MPGE datasets

Dataset group	miRNA	Cell line	Time point	K–S test P-value	De-graded targets	TF mediators	Shuffling P-value
GDS1858	miR-1	HeLa	12	1.2e–15	91	ETS1, CREB1, YY1	0
			24	4.7e–15	107	TFAP2A, CREB1, YY1, SREBF1	0
	miR-124		12	1.8e–12	109	GLI3*	0.04
			24	6.5e–22	132	MLLT7, NKX6.1	0.03
	miR-373		12	9.7e–05	17	MYCN	0.03
			24	4.0e–4	12	NFYA, TAL1, TFAP4*, KLF12	0
GDS2657	miR-124	HepG2	4	4.1e–20	0	AHR, CREB1, SP1, ETS1, EGR1, RELA, KLF12, RREB1, RFX1, NR3C1, BACH2, STAT3	0
			8	3.5e–29	89	AHR*, RELA, RREB1, MEIS1	0
			16	4.7e–49	283	AHR*, CREB1, SP1*, KLF12, RREB1, NR3C1*, BACH2, IRF1	0
			24	1.3e–73	366	AHR*, RREB1	0
			32	6.2e–64	329	AHR*, SP1*, EGR1, RELA*, RREB1, NR3C1*, SP2	0
			72	2.2e–59	292	CREB1, SP1, ETS1, MLLT7, SP2	0
			120	1.1e–19	144	AHR*, SP1*, MLLT7	0
GSE6474	let-7a3	A549	Not known	1.1e–2	1	PAX3, HOXA1, BACH2, EGR3, MYC	0.02
GSE6838	let-7c	HCT116 Dicer−/− #2	24	1.5e–54	211	MYC	0.05
	miR-103		10	5.7e–08	82	MEF2A	0.08
			24	4.4e–22	77	NFATC3, MEF2A	0.04
	miR-106		6	3.4e–55	234	–	–
			10	8.0e–34	158	FOXJ2*	0.05
			24	6.0e–51	246	EGR2	0.07
	miR-107		10	3.2e–06	1	FOXJ2*	0.04
			24	7.8e–15	0	–	–
	miR-15a		6	1.3e–13	0	HOXC8, TBP, POU3F2, FOXC1	0.01
			10	8.3e–51	224	–	–
			14	7.2e–25	0	POU3F2	0.02
			24	3.9e–36	78	WT1	0.06
	miR-15b		10	2.0e–26	0	FOXC1	0.06
			24	4.6e–37	69	FOXC1	0.07
	miR-16		6	2.3e–23	85	–	–
			10	2.6e–29	181	–	–
			14	5.9e–40	0	BACH2	0.11
			24	1.1e–20	46	–	–
	miR-17-5p		24	1.1e–60	215	E2F1*, BCL6, STAT3	0
	miR-192		24	7.8e–07	29	–	–
	miR-195		10	4.3e–30	137	SMAD7, NFATC3	0.05
			24	7.4e–31	98	FOXC1	0.08
	miR-20		24	9.3e–30	59	–	–
	miR-215		24	2.6e–11	38	–	–
GSE7864	miR-34a	A549 H-1 term	24	1.0e–29	112	E2F5*, YY1	0.03
		HCT116 Dicer −/− #2	24	1.1e–29	132	E2F3, YY1, NFE2L1	0.02
		TOV21G H1-term	24	1.4e–23	70	E2F5, BACH2	0.02
		DLD Dicer −/− #2	24	1.6e–33	119	YY1	0.04
		HeLa	24	1.1e–24	128	YY1, BACH2	0.03
	miR-34b	A549 H-1 term	24	1.7e–09	21	E2F5*, MYB*	0.04
		HCT116 Dicer −/− #2	24	4.7e–13	26	E2F5*, MYC, ETS1	0
		TOV21G H1-term	24	3.9e–07	17	MYC, E2F5*, ETS1	0.01
		DLD Dicer −/− #2	24	2.3e–08	18	E2F5*	0.07
		HeLa	24	1.4e–10	40	MYC, BACH2, CREB1, HOXC8	0
	miR-34c	A549 H-1 term	24	1.4e–23	131	E2F3*, MYC, HOXC8	0
		HCT116 Dicer −/− #2	24	1.4e–21	81	E2F3*, MYC, MLLT7, ETS1	0.01
		TOV21G H1-term	24	7.0e–20	37	MYC, TBP, HOXC8, NFE2L1, KLF12, MYB*	0
		DLD Dicer −/− #2	24	4.4e–18	77	–	–
		HeLa	24	7.2e–27	149	MYC, MLLT7, MYB*	0.03

‘K–S test P-value’ is the P-value of a one-sided two-sample Kolmogorov–Smirnov test to assess the degradation-inducing ability of each miRNA. ‘Degraded targets’ gives the total number of degraded targets for each miRNA in each dataset, identified through analysis of the dataset. ‘TF mediators’ lists the TF mediators of miRNA-triggered regulation, identified based on linear regression as shown in Equation (2). The asterisk indicates that the TF is a degraded target. ‘Shuffling P-value’ denotes the proportion of regressed models in the 100 TF-target shuffling experiments that have an identical or greater BIC score than that of the regressed Equation (2).

Many of our predictions are supported by independent experimental studies. For example, our analyses of two MPGE datasets for miR-1 (GDS1858) predicted 130 primary targets; 50 (38.4%) of these targets appeared in TarBase (19), a database collecting experimentally validated miRNA targets. Some miRNAs, like let-7c, miR-16 and miR-17-5p, are associated with only a few targets in TarBase; for the most part, these few targets recurred in our predicted lists. Overall, our screening of putative miRNA–gene relationships resulted in a refined set of mappings with a significantly higher percentage found to be TarBase-validated: 4.3% (121 TarBase-validated mappings out of 2828 relationships in our refined set) versus 1.3% (157 TarBase-validated mappings out of 11 466 putative miRNA–target relationships originally culled from PicTar, TargetScan and miRanda) (see Supplementary Table 7).

Apart from TarBase evidence, our predictions include several that agree with recent literature reports, such as the prediction of MYC as a target of both let-7a-3 and let-7c (14,20).

miRNA’s degradation-inducing ability in human cells

Translational inhibition has long been considered the exclusive mechanism of miRNA-mediated regulation in mammals (21,22). This view, however, may need to be adjusted, especially given recent reports showing that some mammalian miRNAs can down-regulate large numbers of target mRNAs (2,21). Because our work involves a large collection of human MPGE datasets, we hypothesized that we might be able to test the degradation-inducing ability of a number of miRNAs in many different cellular contexts. If putative miRNA targets include a significant portion of degraded targets, whose transcript levels are down-regulated, there should be a significant shift between the MC distribution of putative targets (D2) and that of non-targets (D0.). By contrast, if a miRNA performs its function primarily through translational inhibition, its potential targets should be composed primarily of translationally inhibited targets plus false-positives; the MC distribution of putative targets, D2, degenerates to D0, and hence no significant difference in the two distributions should be found. For this reason, the P-values of the K–S tests of D0 and D2 (Table 1) were taken as an assessment of an miRNA’s degradation-inducing ability. In addition, percent–percent curves similar to quantile–quantile plots were used to demonstrate the difference in D0 and D2 visually (Please find these plots in the ‘wrapped results’ available at http://www.biosino.org/kanghu/~DCR/ supplementary file1.zip).

With the exception of let-7a-3 in the A549 cell line (GSE6474), all miRNAs were found to have degradation-inducing ability in all surveyed situations, as the K–S test P-values were exclusively <0.001 (Table 1). Some miRNAs, such as miR-124, let-7c and miR-106, were associated with hundreds of degraded targets; other miRNAs, such as miR-373, miR-192 and miR-34b, had tens of targets. In some MPGE datasets, although significant degradation-inducing abilities were confirmed for miR-124, miR-107, miR-15a, miR-15b and miR-16 (see K–S test P-values in Table 1), no degraded targets survived the FDR threshold of 0.25 (see the ‘degraded targets’ column). In these cases, small peaks of large negative MCs in the non-target MC distribution (D0, see ‘Materials and Methods’ section) were often observed; their occurrence obscured the detection of degraded targets. Despite these few exceptions, our results suggest that inducing mRNA degradation is a widespread mechanism in human miRNAs.

miRNA-targeted transcription factors and secondary regulation

In our work, a model taking both direct regulation and TF-mediated secondary regulation of one miRNA into account [Equation (2)] was regressed to explain the mRNA level changes upon the miRNA perturbation. Random shuffling of TF–gene regulatory relationships was performed repeatedly to evaluate the significance of the TF-mediated secondary regulations in particular (see ‘Materials and Methods’ section). It turned out that TF-mediated secondary regulations were generally significant in the surveyed MPGE datasets, as an equally valid model could rarely be obtained in the random shuffling experiments (the shuffling P-values were mostly <0.1, see Table 1). The number of identified TF mediators turned out to be far fewer than that of degraded targets; in a network predicted from a MPGE dataset, on average, less than five TF mediators were found responsible for delivering miRNA-initiated regulatory signals to hundreds of secondary targets (Table 1). We identified TF mediators of the miRNA regulatory signal based on a significant difference between MCs of a TF’s target genes and those of its non-target genes. Some of these TF mediators are degraded targets of the miRNA of interest (Table 1); others are not, with trivial or zero MC values. These latter TF mediators are likely to be translationally inhibited, a hypothesis in need of future confirmation by miRNA-perturbed protein-level experiments.

With the transduction of miRNA-initiated regulatory effects through a few TF mediators, a greater number of secondary targets, in addition to primary targets, can be associated with miRNA activity. These secondary targets account for a significant proportion of the observed mRNA level changes in an MPGE dataset (Figure 1C). For example, at the 32-h time point after miR-124 overexpression (dataset from the GDS2657 group), miRNA’s direct regulation could explain decreased MCs of only 181 genes; with our predicted two-layer regulatory model, the decreased MCs of an additional 98 genes and increased MCs of another 42 genes were attributed to miRNA’s indirect regulation, raising the proportion of explainable MCs from 27.8 to 47.7%. The classifications of regulated genes at all time-points in GDS2657 are shown in Supplementary Table 8, where a general trend is evident that the direct regulation decreases rapidly while the secondary regulation maintains at a considerable multitude, resulting in an actual dilution of direct regulation in the total regulations as time passes. This trend is in accordance with the data releasers’ observation (13) and our expectation. In terms of the number of regulated secondary targets, TF-mediated effects reached their summit at around 16-h post transfection. To further address the timing of secondary regulatory effects mediated by TFs, we performed a hypergeometric test of the consistency between early weak up-regulations and later strong up-regulations. The results, provided in Supplementary Text 9, indicate that the shift from random noise to genuine secondary up-regulations probably takes place between 8 and 16 h in MPGE datasets.

Two-layer regulatory networks centered on miRNA and mediated by TFs

Figure 3 depicts a typical two-layer regulatory network, mined from an MPGE dataset measured at the 12-h time point after overexpression of miR-1 (dataset from the GDS1858 group). In addition to directly down-regulating 91 degraded targets (blue arrows), miR-1 overexpression causes expression changes in more than 100 non-target genes, possibly through translationally inhibiting three TFs (ETS1, CREB1 and YY1; red arrows).

Figure 3.

A two-layer regulatory network centered around miR-1, mapped based on analysis of the 12-h MPGE dataset in the GDS1858 dataset group. The red circle indicates miR-1; red arrows represent regulatory effects of miRNAs on TFs; blue solid arrows indicate regulation of target genes by miRNA-induced degradation of corresponding mRNAs; black dashed arrows represent regulatory effects of TFs on their target genes; light green circles indicate down-regulated genes; pink circles indicate up-regulated genes; and gray circles indicate TFs translationally inhibited by miRNA. This network was found to be most significantly associated with three biological themes: multicellular organismal development, cell development and cell–cell signaling.

According to enrichment analyses of GO and KEGG terms (see ‘Materials and Methods’ section), the miR-1 network is associated with 196 biological themes, such as cell development (FDR = 0.12) and cell–cell signaling (FDR = 0.12). Notably, six GO terms on muscle development or regulation of muscle contraction were enriched in this network (GO:0045987, GO:000693, GO:004593, GO:000694, GO:005114, GO:000752; all FDRs < 0.29), consistent with the known fact that miR-1 is expressed selectively in heart and skeletal muscle. Similar analyses were performed on the miR-124 network (GDS2657, 32 h), resulting in the identification of 129 significant biological themes (FDR < 0.25), among which neuron apoptosis (GO:0051402) is in accordance with miR-124’s proven role in development of the nervous system (23). Detailed lists of biological themes significantly associated with each network are included in the ‘wrapped results’ (http://www.biosino.org/~kanghu/DCR/ supplementary file1.zip).

Other miRNA-centered networks also were found to be significantly enriched with certain GO terms and KEGG pathways. Some networks, such as those triggered by miR-17-5p and miR-34 family members, were enriched with quite a lot of biological themes, indicative of their active participation in a broad scope of cellular processes (Supplementary Table 10); certain biological themes (such as cell cycle, regulation of kinase activity and cell proliferation) appeared repeatedly in many miRNA networks, indicative of the complex nature of their regulation by a number of miRNAs and TFs (Supplementary Table 10). A global view of the biological themes regulated by miRNAs in different cellular contexts is shown in a heatmap (Supplementary Figure 11).

Finally, individual networks were integrated into a global two-layer combinatorial network, involving 19 miRNAs, 45 TFs, 1,685 direct target genes, 2755 degradation-based direct regulation effects, and 53 translational-inhibition-based direct regulation effects. A topological graph of the network was developed using Cytoscape (24) (see Supplementary Figure 12). In agreement with similar analysis in Caenorhabditis elegans (25), the degrees of target genes in the global network manifested an exponential distribution with a slope of about −0.38 and R2 = ∼0.9978 (Supplementary Figure 13), indicating that most vertices (96.8%) were regulated by very few miRNAs (up to four), while few targets (3.2%) had many regulators (five to eight). Our web service allows online query and analysis of this global regulatory network (http://www.biosino.org/de.cleaved.repressed).

A tumor-related, miRNA/TF-mediated regulatory network

miRNAs are speculated to have an intrinsic function in tumor suppression, as a global decrease in miRNA levels is often observed in human cancers (26). Recent studies have related several miRNAs with the well-known oncogene MYC and tumor suppressor gene TP53 (14,27), highlighting an important role for miRNA in the regulation of tumorigenesis. According to the Tumor Suppressor and Oncogene Directory (http://embryology.med.unsw.edu.au/DNA/DNA10.htm) and a manual literature search, 12 of the 45 total TFs mapped to the global combinatorial network described above are either oncogenes or tumor suppressors, providing further evidence of miRNA’s close relationship with tumorigenesis. Some miRNAs known to suppress tumorigenesis, such as the miR-34 family (16) and let-7 family (14,27), were found to inhibit MYC, MYB and E2F family oncogenes (Table 1). In contrast, the miR-15 family was found to inhibit two tumor suppressors, FOXC1 and WT1 (Table 1), suggesting a possible tumor-induction effect.

Based on existing knowledge and our novel discoveries, we propose that the regulatory cascades triggered by the miR-34, let-7, miR-17 and miR-1 families may be linked together to form a tumor-related regulatory network (Figure 4). In addition to confirming the recent discovery that the miR-34 family are direct transcriptional targets of p53 (16), our results showed further that miR-34 family miRNAs regulate 13 downstream TFs. Five of these 13 TFs are oncogenes: E2F3, E2F5, MYB, MYC and ETS1. Previous reports have experimentally confirmed E2F family members as targets of miR-34a (28), but our study is the first to identify MYC, MYB and ETS1 as targets of the miR-34 family. MYC bridges the miR-34 sub-network and the let-7 sub-network, as a target for 7a-3 and let-7c, as well as for miR-34. In Figure 4, the bidirectional solid line between MYC and let-7 represents a reciprocal regulation mechanism: MYC promotes let-7 transcription (14), and let-7 inhibits MYC in return (14,27). The integration of the oncogene-enriched miR-17 sub-network with the tumor-related network also occurs via MYC, which has been reported to induce miR-17 transcription (29). All three TF targets of miR-17 (BCL6, E2F1 and STAT3) are oncogenes, adding significantly to the tumor-suppressing effect of the composite network. Note that inhibition of E2F1 by miR-17 has been reported previously (29). Finally, a predicted target of miR-1, YY1, has been shown to down-regulate TP53 (30), a regulator of the miR-34 family. This connects the miR-1 sub-network with the tumor-related network; miR-1 further connects to the tumor-related network via the predicted oncogene target ETS1, a target shared with the miR-34 family. Overall, the composite network mapped here connects two TFs (TP53 and MYC) to three miRNA families (miR-34, miR-17 and let-7) via three transcriptional regulations, and four miRNA (miR-34, miR-17, miR-1 and let-7) families to 22 TFs via 27 post-transcriptional regulations. This network has a strong association with the regulation of tumorigenesis, because nine of the 22 total TFs are oncogenes or tumor suppressors, with the important tumor-related regulators TP53 and MYC playing hub connector roles. Our results reinforce the growing awareness that non-coding RNAs may play key roles in tumor development, and might bring new insights to the mechanisms of p53 and MYC.

Figure 4.

A miRNA/TF-mediated regulatory pathway associated with tumorigenesis. miRNAs (miR-34, miR-17, let-7 and miR-1) are shown as small pink circles. Oncogene TFs (MYC, MYB, E2F1, E2F3, E2F5, BCL6, ETS1 and STAT3) are shown as red circles. The tumor suppressor gene TP53 is shown as a green circle. Other miRNA-regulated TFs are shown as light blue circles. Red and green arrows indicate activation and inhibition, respectively. Predicted regulations are shown as dashed arrows, while validated regulations are shown as solid arrows.

DISCUSSION

In this work, we developed algorithms for mining miRNA-initiated primary and secondary regulations, through the analysis of 53 gene expression datasets involving overexpression of 19 miRNAs in four types of human cell lines. As well as illustrating that miRNA-induced degradation of targets is widespread in human cell lines, these algorithms allow for refined mining of such targets, which will help increase the reliability of miRNA target prediction. The mining of secondary regulations for the first time in MPGE datasets sheds light on the scope of the extended regulatory effects of miRNA through important mediators: transcription factors. Combining primary and secondary regulations allowed for depiction of a two-layer network including all 19 miRNAs analyzed. In the future, combinatory networks centered around individual miRNAs or specific biological themes (e.g. tumor-related) could be analyzed in greater detail. A web server is provided for query and analysis of all reported data in this article.

miRNAs are conjectured to silence gene expression mainly through mRNA degradation if they base pair perfectly with target mRNA sequences, or translational repression if they display imperfect sequence complementarity with target mRNAs. MPGE experiments have been performed to help identify targets of miRNA-induced degradation. In this regard, MPGE datasets are beneficial experimental supplements to conventional complementarity-based target prediction. Many computational studies have been carried out to utilize MPGE datasets to improve the prediction of miRNA targets (7–9). In our work, a novel strategy was proposed to identify true degraded miRNA targets based on seed-region match and expression change patterns revealed in miRNA overexpression datasets. A non-parametric statistical test of mRNA level change against a null distribution, in which the existence of true degraded targets was presumably denied, was shown to be an effective way of distinguishing highly likely degraded targets from a set of all putative targets predicted by seed-region match-based algorithms. Our refined set of degraded targets is more reliable than the full set of putative targets, based on validation by TarBase records. Moreover, we proposed to determine miRNA’s overall degradation-inducing ability by discriminating the mRNA level change distribution of putative targets and that of putative non-targets. Results of our statistical tests suggested that induction of mRNA degradation is a ubiquitous mechanism of human miRNAs in observed cell lines.

In general, presence (or absence) of a miRNA seed-match sequence in mRNA 3′ UTRs seems to be correlated with relative down-regulation (or up-regulation) (31), but many mRNA level changes remain to be illuminated; examples include the up-regulation of seed-matched targets and the regulation of non-targets, such that the percentage of predicted targets among down-regulated genes decreases rapidly with time from miRNA perturbation (13). In this work, we attempted to address this issue by considering miRNA’s secondary regulatory effects transmitted by target TFs. Similar to a forgoing work (31), we tried to model mRNA level changes with a linear regression, but we made an innovation by adopting TFs, in addition to miRNAs of interest, in explanation of mRNA level changes. With this innovation, the number of explainable mRNA level changes in MPGE datasets increased significantly compared to previous reports. By adopting TF targets as mediators of miRNA secondary regulation, it is also possible to mine translationally inhibited TF targets, based on systematic mRNA changes in secondary target genes from miRNA overexpression datasets; this has never been done before.

One miRNA can target hundreds of genes. The mechanism of such extended effects has long been a hot topic of research. Increasing evidence links miRNAs to transcription factors or complexes. For example, in regulation of stress-dependent cardiac growth and gene expression by miR-208, strong experimental evidence suggests that miR-208 acts partly by repressing expression of THRAP1, which can exert positive and negative effects on transcription (32). Based on a similar perspective, we envisioned transcription factors to be important mediators of secondary regulation by miRNAs in our study of MPGE datasets. MPGE datasets, however, have never been used to mine secondary regulations. Secondary regulations are usually considered to happen later than 24 h after miRNA transfection. This is true if a stringent cutoff threshold is applied to mine secondary up-regulated targets that cannot be explained by primary miRNA regulation. Contrary to the threshold-based dichotomous view, in our work, an opposite perspective was adopted through use of a linear regression algorithm [Equations (1) and (2)], which takes advantage of the log ratios (relative expression levels) themselves that reflect the continuous shifting in relative expression level from non-regulated extreme to regulated extreme. Specifically, a TF is identified as a potential secondary effect mediator if a significant difference between the log ratios of its target genes and those of its non-target genes is observed. The difference may be small, but considering the large sample size, it is reasonable that a small but systematic difference could be detected. As discussed in the ‘Results’ section, we found that MPGE datasets after 8 h could be used to discriminate weak secondary effects from random noise. Thus, the utility of MPGE datasets could be extended to studying miRNA–TF–gene cascades or feed-forward loops and eventually constructing a two-layer regulatory network.

In this report, novel methodologies were proposed to infer miRNA-triggered cellular regulatory networks through integrated analyses of MPGE datasets and putative miRNA–TF and TF–gene regulatory relationships. We are confident in the reliability of our results, because of the observations concordant with independent studies and the recapitulation of expected functional roles of perturbed miRNA (Supplementary Table 7 and Supplementary Figure 11). TFs repressed by miRNA were found to be important mediators that propagate miRNA’s regulatory effects to downstream targets. In most cases, each miRNA affected fewer than five TFs, which then influenced hundreds of secondary targets, making up a miRNA–TF combinatorial regulatory network. Recent studies have indicated that tissue-specific miRNAs may function at multiple hierarchical levels of gene regulatory networks, from targeting hundreds of effector genes incompatible with the differentiated state to controlling the levels of global regulators of transcription. This multilevel regulation may allow individual miRNAs to profoundly affect the gene expression program of differentiated cells (33). Our work reported here represents both a good example and an important approach in such an exciting research field.

As an initiative to uncover miRNA-triggered regulatory cascades, our work has limitations. One miRNA could potentially target hundreds of genes, but the exact silencing mechanism is not static and the regulated targets are not constant, as regulation could depend on other RNA-binding proteins or specific tissue/cell types. Our work is limited to TFs and performed on public MPGE datasets from, primarily, established cancer cell lines. In addition, because we lack protein-level validation data, our identification of potential translationally inhibited TF targets should be considered as hypotheses. With these caveats in mind, we would encourage the preparation of more appropriate datasets to ensure optimal analysis results. With an anticipated expansion of MPGE datasets and proteomics data in the near future, our strategy could be applied to reveal conditional regulatory pathways in many more cellular contexts, whose integration may lead to elucidation of a reference regulatory network in human cells (27).

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.

FUNDING

National High-Tech R&D Program (863: 2006AA02Z334, 2007AA02Z304, 2006AA02Z330, 2009AA02Z304); International collaboration program (2007DFA31040); Key Research Program of CAS (KSCX2-YW-R-112); China National Key Projects for Infectious Disease (2008ZX10002-021); Shanghai Natural Science Foundation (08ZR1415800), Science and Technology Commission of Shanghai Municipality (2008PJ14084), and the National Natural Science Foundation of China (30770497). Funding for open access charge: National key basic research program (973:2006CB910700) from Ministry of Science and Technology of China.

Conflict of interest statement. None declared.