Documented changes in levels of microRNAs (miRNA) in a variety of diseases including cancer are leading to their development as early indicators of disease, and as a potential new class of therapeutic agents. A significant hurdle to the rational application of miRNAs as therapeutics is our current inability to reliably predict the range of molecular and cellular consequences of perturbations in the levels of specific miRNAs on targeted cells. While the direct gene (mRNA) targets of individual miRNAs can be computationally predicted with reasonable degrees of accuracy, reliable predictions of the indirect molecular effects of perturbations in miRNA levels remain a major challenge in molecular systems biology.
Changes in gene (mRNA) and miRNA expression levels between normal precursor and ovarian cancer cells isolated from patient tissue samples were measured by microarray. Expression of 31 miRNAs was significantly elevated in the cancer samples. Consistent with previous reports, the expected decrease in expression of the mRNA targets of upregulated miRNAs was observed in only 20-30% of the cancer samples. We present and provide experimental support for a network model (The Transcriptional Override Model; TOM) to account for the unexpected regulatory consequences of modulations in the expression of miRNAs on expression levels of their target mRNAs in ovarian cancer.
The direct and indirect regulatory effects of changes in miRNA expression levels in vivo are interactive and complex but amenable to systems level modeling. Although TOM has been developed and validated within the context of ovarian cancer, it may be applicable in other biological contexts as well, including of potential future use in the rational design of miRNA-based strategies for the treatment of cancers and other diseases.
Cancer systems biologyFeed-forward loopsGene regulationmiRNAsOvarian cancer
Human miRNAs regulate gene expression post-transcriptionally by degrading target mRNAs and/or blocking their translation . As a consequence, mRNA expression changes are expected to be inversely correlated (IC) with changes in levels of their targeting miRNAs. Although this expectation has been validated in studies of individual miRNAs and specific mRNA targets, the expected inverse relationship is often not observed in global transcriptome level studies [2–4]. While these unexpected findings may, in some instances, be attributed to inaccuracies in miRNA target prediction algorithms , recent evidence suggests that many of the unexpected regulatory effects may be the result of feed-back or feed-forward loops and/or other system level complexities [3, 6].
To systematically address the relationship between miRNAs and their regulated mRNA targets in the same cellular context, we employed microarray gene expression profiling to compare differences in expression levels of mRNAs and miRNAs in ovarian surface epithelial cells (OSE) vs. serous papillary ovarian cancer epithelial cells (CEPI) isolated from patient tissues (Additional file 1) by laser capture microdissection. Gene expression profiling identified 5910 significantly differentially expressed genes (mRNAs) between OSE and CEPI (Additional file 2). Of these, 2232 (38%) were upregulated and 3678 (62%) downregulated in CEPI. MiRNA expression profiling identified 31 significantly differentially expressed miRNAs between OSE and CEPI. All of these miRNAs were upregulated in CEPI (Table 1).
MiRNA expression profiling identified 31 significantly differentially expressed miRNAs between OSE and CEPI
Expression profiles for miRNAs in OSE and CEPI samples were determined using the miRChip system (Asuragen Inc.). Human probe sets with >65% present calls in either of the two groups (OSE and CEPI) were selected for analysis. All reported microarray data are described in accordance with MIAME guidelines. The processed and raw data files for the samples used in this study have been deposited in the Gene Expression Omnibus (GEO) as SuperSeries GSE42460.
Employing three commonly used miRNA target prediction algorithms (miRanda-mirSVR , TargetScan , SVMicrO ), we identified putative mRNA targets of these 31 miRNAs to determine if differences in their levels of expression between the OSE and CEPI samples were IC, positively correlated (PC) or unchanged (NC). Based on the established molecular mechanism of miRNA regulation, levels of miRNAs are expected to be IC with levels of their target mRNAs. Contrary to this expectation and consistent with previous findings , we observed a consistently low percentage (23-31%) of target mRNAs displaying expression level changes IC with their regulating miRNAs (Table 2). Since predictions from the most commonly employed mirSVR algorithm were representative of the results from all tested algorithms, mirSVR was employed in subsequent analyses.
Percentage of miRNA target genes displaying changes in expression between OSE and CEPI
No change detected
Predicted gene targets of miRNAs were computed independently using three commonly used prediction algorithms: miRanda-mirSVR , TargetScan  and SVMicrO . For each prediction algorithm, the frequency of target genes displaying changes in gene expression IC, PC or NC with changes in expression of their regulating miRNAs are shown. The number of predicted target genes falling into each class is presented in parentheses.
The transcriptional override model (TOM) postulates that downregulation of target genes induced by elevated levels of regulating miRNAs may be masked (NC) or overridden (PC) by increases in expression mediated by the downregulation of repressor genes that are themselves targets of upregulated miRNAs (Figure 1A). The possibility of such feed-forward loops was prompted by the fact that several of the predicted mRNA targets of the 31 overexpressed miRNAs encode documented repressors of gene expression (e.g. ZNF24, YY1, SPEN, BACH1). MiRNA-mediated downregulation of these repressor genes would be expected to result in the derepression of their respective target genes and a consequent increase in levels of expression. If these derepressed gene targets were also the targets of upregulated miRNAs, the expected downregulation of these genes by the miRNAs (IC) could be masked (NC) or overridden (PC). For example (Figure 1B), one of the predicted targets of ten of the 31 miRNAs upregulated in cancer is the well-documented repressor gene ZNF24. Consistent with the fact that ZNF24 is targeted by upregulated miRNAs, its expression in CEPI is significantly reduced. An experimentally validated target of ZNF24 is VEGFA. Despite the fact that VEGFA is itself directly targeted by 11 upregulated miRNAs (including five of those targeting ZNF24), its level of expression is significantly increased (PC) in CEPI. These results are consistent with the hypothesis that ZNF24-mediated derepression is overriding the expected downregulatory effects of the upregulated miRNAs on VEGFA expression.
Although many of the genes falling within the NC category could be the result of partial transcriptional override, they might also simply be the result of no or slight miRNA regulatory effects. Since we cannot experimentally distinguish between these two possibilities, we will operationally only consider PC differences in expression as being inconsistent with the expected IC differences.
To further evaluate TOM, we identified 105 genes that are 1) targets of one-or more of the 31 upregulated miRNAs, 2) significantly downregulated in our cancer samples and 3) previously characterized as transcriptional repressors  (Additional file 3). The targets of ten (Table 3) of these 105 genes have been previously identified in the Transcription Factor Binding Site database (TRANSFAC) . This resulted in 843 genes (Additional file 4) predicted to be directly targeted by both the ten downregulated repressors and one-or-more of the 31 upregulated miRNAs. From the perspective of miRNA regulation, all 843 of these target genes are expected to be downregulated while from the perspective of the downregulated repressor genes all the targets are expected to be upregulated. The observed reality lies somewhere between these two expectations (Figure 2). TOM predicts that the response of any particular target gene will be determined by the relative strengths of these two opposing regulatory controls.
Ten genes characterized as validated repressors and predicted targets of one or more of the 31 miRNAs upregulated in CEPI
The number of miRNAs targeting individual human genes (mRNAs) is known to vary from zero to over 100 with an estimated average of 7.1 miRNA targets per gene . Thus, the relative strength of the regulatory effect of miRNAs on target genes might be expected to be a function of the number of miRNAs targeting individual genes. To explore this possibility, we grouped the 843 predicted gene targets of the 31 upregulated miRNAs and the 10 downregulated repressors into 5 groups based upon the number of upregulated miRNAs predicted to target each gene (Figure 3A). A sixth group was comprised of targets of the ten repressor genes (TRANSFAC) but not predicted targets of any of the 31 miRNAs. For each group, we computed the percentage of target genes that were upregulated (PC). In addition, for each group we divided the target genes into those predicted to be regulated by a single repressor vs. those predicted to be regulated by multiple repressors (Figure 3B).
The results demonstrate that as the number of upregulated miRNAs targeting individual genes increases, the percentage of target genes displaying the unexpected PC change in expression decreases. These results are consistent with TOM and indicate that as the relative strength of the miRNA regulatory effect increases, the impact of the opposing derepression effect mediated by the downregulated repressor genes is diminished. However, the fact that ~20% of genes targeted by even large numbers (>15) of upregulated miRNAs continue to display the unexpected PC indicates that, in some cases, the magnitude of derepression is sufficient to completely override miRNA regulation. The results presented in Figure 3B suggest that genes targeted by multiple repressors tend to be associated with a higher percentage of PC genes than those targeted by a single repressor. The effect, however, is not as consistent as observed with increasing numbers of regulating miRNA likely due to the relatively low number of repressor genes in this dataset and the fact that not all repressor genes can be expected to exert the same magnitude of regulatory control.
We were next interested to see if the model’s ability to account for trends observed using the limited dataset described above might also extend more globally. We divided all differentially expressed genes including those that are predicted gene targets of the 31 upregulated miRNAs (4829) and those that are not (1081), into 6 groups based on the number of miRNAs targeting each gene. The results (Figure 4A) demonstrate a clear inverse relationship between the number of miRNAs targeting genes and the percentage of these genes displaying the unexpected PC. Again, however, we found that ~20% of genes targeted by even large numbers (>15) of upregulated miRNAs continue to display the unexpected PC consistent with the hypothesis that the magnitude of repressor gene mediated derepression is, in some instances, sufficient to completely override miRNA regulation.
Testing the model’s ability to globally predict the relative influence of miRNA and repressor gene regulatory controls on target gene expression is problematic for two reasons: first, a compendium of all human repressor genes and their regulatory targets is currently unavailable; second, many regulatory proteins can function as repressors or activators depending on cellular context and protein complex association . One approach taken by systems biologists to model regulatory relationships in complex cellular contexts is to use highly correlated changes in expression patterns among genes as evidence of direct and/or indirect interactions [18, 19]. In our case, we examined variation in gene expression patterns across our OSE samples to identify genes displaying consistent inverse correlations (Pearson correlation coefficient < -0.8) in expression with changes in expression of the 105 repressor genes previously characterized as significantly downregulated in CEPI and regulatory targets of one-or-more of the 31 miRNAs (see above). Genes displaying an inversely correlated pattern of co-expression (1205) were operationally classified as targets of these repressor genes. Genes not displaying this pattern of expression (3624) were classified as non-targets of the designated repressor genes. Having established these classes, it became possible to distinguish between regulatory interactions fulfilling the triangular relationship of TOM from those that do not.
The results presented in Figure 4B again indicate that as the number of upregulated miRNAs targeting genes increases, the percentage of the unexpected PC decreases significantly (chi-square test for trend X2 = 444.6, p < <0.0001). Consistent with the results presented in Figure 4, the overlap of miRNA targets and downregulated genes is highly significant. (Figure 5A; hypergeometric p < 1E-12). Likewise, the overlap of repressor targets and upregulated genes is also highly significant (Figure 5B; hypergeometric p < 1E-12). Collectively, the results indicate that the derepression of target genes mediated by high (6–45) numbers of downregulated repressors is sufficient to nearly or completely override the regulatory controls of even large numbers (>15) of upregulated miRNAs.
Recent studies have clearly established miRNAs as early indicators of disease [20, 21] and as a potential new class of therapeutic agents [22, 23]. Full appreciation of the biological significance of modulations in levels of miRNAs, as well as, the future rational employment of miRNAs as therapeutic agents will require an understanding of both the direct and indirect molecular consequences of changes in the levels of miRNAs on cell function. While the direct gene (mRNA) targets of individual miRNAs can be computationally predicted and experimentally validated with varying degrees of accuracy , reliable predictions of the indirect molecular effects of changes in miRNA levels has remained a major challenge in molecular systems biology [23, 25].
In this paper, we present a regulatory network model (TOM) that explains a significant component of the unexpected low frequency of IC changes in expression levels between mRNAs and their regulating miRNAs. The model postulates that the expected downregulation of target genes induced by elevated levels of regulating miRNAs may be masked or “overridden” by increases in transcriptional initiation mediated by the downregulation of repressor genes that are themselves targets of the same regulating miRNAs (Figure 1A). Depending upon the strength of the transcriptional override (i.e., the relative strengths of miRNA and repressor gene mediated de-repression), TOM predicts that increases in miRNA levels may display no effect (NC) or be positively correlated (PC) with changes in levels of their targeted mRNAs.
It is widely recognized that the operation of regulatory effects mediated by miRNAs in vivo is a complex and interactive process and a number of explanatory models have been previously offered [26–29]. In this paper, we propose an additional model (TOM) that focuses on the feedback interactions that exist between miRNAs, regulatory (repressor) genes and their mutual gene targets. While we have evaluated TOM within the context of its ability to account for global patterns of changes in gene expression, the model also provides a framework for predicting interactions between specific miRNAs and target genes (e.g., Figure 1). Further testing in other cancers (and other biological contexts) will be needed to evaluate the robustness of TOM. Nevertheless, our initial findings in ovarian cancer indicate that interactions between miRNAs and repressor genes may well play a significant role in effecting the unexpected regulatory responses of targeted genes to modulations in levels of their regulatory miRNAs.
It is now widely acknowledged that the complexity of molecular interactions taking place on the cellular level can significantly obscure the expected consequences of molecular processes characterized in vitro [30, 31]. Our findings indicate that the direct and indirect regulatory effects of changes in miRNA expression levels in vivo are interactive and complex but amenable to systems level modeling. We have shown that TOM can account for a major component of the unexpected consequences of changes in miRNA expression levels on their target mRNAs. Although the model has been developed and evaluated within the context of ovarian cancer, we believe it may be applicable in other biological contexts as well including of potential future use in the rational design of miRNA-based strategies for the treatment of cancer and other diseases.
All tissues were collected according to previously published procedures  following approved Institutional Review Board protocols from Northside Hospital (Atlanta) and Georgia Institute of Technology. Informed consent was obtained from all subjects. The histopathology for all cancer patients was serous papillary adenocarcinoma of the ovary and for the control patients the ovaries were considered within normal limits.
mRNA microarray data analysis
Ten OSE (normal) and ten CEPI (cancer) samples were analyzed for mRNA expression using the Affymetrix Gene Chip Operating System (GCOS HG-U133 Plus 2.0). CEL files generated by GCOS were converted to expression values using GCRMA normalization on the arrayanalysis.org  website, which output also included quality control metrics, principal components analysis (PCA) and cluster dendrograms. Present/absent calls were generated from the MAS 5.0 statistical algorithm as implemented in Affymetrix Expression Console. Probe sets with >60% present calls in either of the two groups (OSE and CEPI) were selected for further analysis. After log2 transformation, signal values of those probe sets were submitted to Statistical Analysis of Microarrays (SAM) for multiple testing correction where a 5.5% FDR was applied resulting in 7462 probe sets representing 5910 differentially expressed genes (DEGs). Annotations for probe sets were obtained from Affymetrix . The processed and raw data files for the samples used in this study have been deposited in the Gene Expression Omnibus (GSE52037 with SuperSeries GSE52460).
microRNA microarray data analysis
Expression profiles for microRNAs from three OSE and three CEPI samples were generated by Asuragen (Austin, TX) using Ambion miRChip technology (Life Technologies). Two sets of CEL files, created from 6 biological replicates and two sets of technical replicates were normalized using MAS 5.0 to expression signals, giving 6 values per probe/gene. Probe sets labeled as human (those having an “hsa-” prefix), known to be conserved to mouse, and with at least 65% present calls (calculated by Asuragen) in either of the two groups (OSE and CEPI) were selected for further filtering. Thirty-one differentially expressed microRNAs (fc > 6, p-value < .03) were selected. The repressive potential of all 31 microRNAs was validated by noting that > 65% of the predicted DEG targets of each upregulated microRNA were actually downregulated, while only 44% of DEGs not predicted to be targets of any upregulated microRNA were downregulated. Mean repression over all 31 microRNAs was 71%. The processed and raw data files for the samples used in this study have been deposited in the Gene Expression Omnibus (GSE52459 with SuperSeries GSE52460).
microRNA target prediction
The miRNA target prediction file based on mirSVR was downloaded from microRNA.org (August 2010 release). The mirSVR score refers to targets of microRNAs with scores obtained from their support vector regression algorithm. To reduce the occurrence of false positives, only predicted targets with a mirSVR score less than -.2 were considered. The microRNA target predictions based on TargetScan and SVMicrO were downloaded from http://www.targetscan.org (retrieved 8/2010) and http://www.compgenomics.utsa.edu/Result/Human/hsa_human (retrieved 9/2010), respectively.
Transcriptional repressor selection
Members of the Gene Ontology categories GO:0045892, GO:0000122, GO:0010944, GO:0032088 and GO:0008156, relating to the negative-regulation-of-transcription or its child terms, were downloaded from the European Bioinformatics Institute (EBI) and parsed using UNIX scripts. In that download, we found 439 potential repressor genes. Of those, 109 genes were significantly downregulated according to our microarray analysis and 105 of these genes were also predicted targets of one or more of the 31 upregulated microRNAs. These 105 transcriptional repressor genes formed the basis for microRNA target derepression in our model.
Transcriptional repressor target prediction and experimental validation
To obtain predicted and/or experimentally validated transcription factor binding site data, we downloaded the TRANSFAC data file c3.tft.v3.1.symbols.gmt from GSEA (Gene Set Enrichment Analysis website - http://www.broadinstitute.org/gsea/downloads.jsp). Data files were parsed with UNIX scripts, which extracted pairs of genes consisting of one repressor and one or more binding partners. All repressor-partner pairs under consideration had to be DEGs and predicted targets of at least one of the 31 upregulated microRNAs, and all transcriptional repressors were downregulated in cancer. Further, all repressor-partner pairs were required to show a correlation coefficient of r < -.8 across all normal samples.
Correlation coefficient calculation
For the global analysis of relationships among all 105 transcriptional repressors and their binding partners, Pearson’s correlation coefficient (PCC) was calculated across all ten OSE (normal) samples between all transcriptional repressors and predicted microRNA targets. Specifically, we used the Mathematica  correlation function (n = 10; r < -.8) for a directional significance of (p < .0027). Fold-change from normal to cancer in these genes ranged from -625 to 121.
Availability of supporting data
The processed and raw data files for the samples used in the mRNA and miRNA expression studies have been deposited in the Gene Expression Omnibus (GSE52037 and GSE52459 with SuperSeries GSE52460).
Cancer epithelial cells
Ovarian surface epithelial cells
Funding for this project was provided by Ovarian Cycle, Deborah Nash Endowment Fund, Josephine Robinson Family, and the J.D. Rhodes Trust.
Integrated Cancer Research Center, School of Biology, and Parker H. Petit Institute of Bioengineering and Biosciences, Georgia Institute of Technology
Shahab SW, Matyunina LV, Mezencev R, Walker LD, Bowen NJ, Benigno BB, McDonald JF: Evidence for the complexity of MicroRNA-mediated regulation in ovarian cancer: a systems approach.PLoS One 2011, 6:e22508.PubMed CentralPubMedView Article
Lim LP, Lau NC, Garrett-Engele P, Grimson A, Schelter JM, Castle J, Bartel DP, Linsley PS, Johnson JM: Microarray analysis shows that some microRNAs downregulate large numbers of target mRNAs.Nature 2005, 433:769–773.PubMedView Article
Liu H, Yue D, Chen Y, Gao SJ, Huang Y: Improving performance of mammalian microRNA target prediction.BMC Bioinforma 2010, 11:476.View Article
Jia D, Hasso SM, Chan J, Filingeri D, D’Amore PA, Rice L, Pampo C, Siemann DW, Zurakowski D, Rodig SJ, Moses MA: Transcriptional repression of VEGF by ZNF24: mechanistic studies and vascular consequences in vivo.Blood 2013, 121:707–715.PubMed CentralPubMedView Article
Bonavida B, Baritaki S: The novel role of Yin Yang 1 in the regulation of epithelial to mesenchymal transition in cancer via the dysregulated NF-kappaB/Snail/YY1/RKIP/PTEN Circuitry.Crit Rev Oncog 2011, 16:211–226.PubMedView Article
Ariyoshi M, Schwabe JW: A conserved structural motif reveals the essential transcriptional repression function of Spen proteins and their role in developmental signaling.Genes Dev 2003, 17:1909–1920.PubMed CentralPubMedView Article
Kitamuro T, Takahashi K, Ogawa K, Udono-Fujimori R, Takeda K, Furuyama K, Nakayama M, Sun J, Fujita H, Hida W, Hattori T, Shirato K, Igarashi K, Shibahara S: Bach1 functions as a hypoxia-inducible repressor for the heme oxygenase-1 gene in human cells.J Biol Chem 2003, 278:9125–9133.PubMedView Article
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G: Gene ontology: tool for the unification of biology. The gene ontology consortium.Nat Genet 2000, 25:25–29.PubMed CentralPubMedView Article
Matys V, Kel-Margoulis OV, Fricke E, Liebich I, Land S, Barre-Dirrie A, Reuter I, Chekmenev D, Krull M, Hornischer K, Voss N, Stegmaier P, Lewicki-Potapov B, Saxel H, Kel AE, Wingender E: TRANSFAC and its module TRANSCompel: transcriptional gene regulation in eukaryotes.Nucleic Acids Res 2006, 34:D108-D110.PubMed CentralPubMedView Article
Remenyi A, Scholer HR, Wilmanns M: Combinatorial control of gene expression.Nat Struct Mol Biol 2004, 11:812–815.PubMedView Article
Mansson R, Tsapogas P, Akerlund M, Lagergren A, Gisler R, Sigvardsson M: Pearson correlation analysis of microarray data allows for the identification of genetic targets for early B-cell factor.J Biol Chem 2004, 279:17905–17913.PubMedView Article
Allocco DJ, Kohane IS, Butte AJ: Quantifying the relationship between co-expression, co-regulation and gene function.BMC Bioinforma 2004, 5:18.View Article
Shahab SW, Matyunina LV, Hill CG, Wang L, Mezencev R, Walker LD, McDonald JF: The effects of MicroRNA transfections on global patterns of gene expression in ovarian cancer cells are functionally coordinated.BMC Med Genomics 2012, 5:33.PubMed CentralPubMedView Article
Gurtan AM, Sharp PA: The role of miRNAs in regulating gene expression networks.J Mol Biol 2013, 425:3582–3600.PubMedView Article
Ritchie W, Rasko JE, Flamant S: MicroRNA target prediction and validation.Adv Exp Med Biol 2013, 774:39–53.PubMedView Article
Vera J, Lai X, Schmitz U, Wolkenhauer O: MicroRNA-regulated networks: the perfect storm for classical molecular biology, the ideal scenario for systems biology.Adv Exp Med Biol 2013, 774:55–76.PubMedView Article
Li Y, Goldenberg A, Wong KC, Zhang Z: A probabilistic approach to explore human miRNA targetome by integrating miRNA-overexpression data and sequence information.Bioinformatics 2014, 30:621–628.PubMedView Article
Arvey A, Larsson E, Sander C, Leslie CS, Marks DS: Target mRNA abundance dilutes microRNA and siRNA activity.Mol Syst Biol 2010,2010(6):1–7.
Khan AA, Betel D, Miller ML, Sander C, Leslie CS, Marks DS: Transfection of small RNAs globally perturbs gene regulation by endogenous microRNAs.Nat Biotechnol 2009, 27:549–555.PubMed CentralPubMedView Article
Hornberg JJ, Bruggeman FJ, Westerhoff HV, Lankelma J: Cancer: a systems biology disease.Biosystems 2006, 83:81–90.PubMedView Article
Bowen NJ, Walker LD, Matyunina LV, Logani S, Totten KA, Benigno BB, McDonald JF: Gene expression profiling supports the hypothesis that human ovarian surface epithelia are multipotent and capable of serving as ovarian cancer initiating cells.BMC Med Genomics 2009, 2:71–84.PubMed CentralPubMedView Article
Eijssen LM, Jaillard M, Adriaens ME, Gaj S, de Groot PJ, Muller M, Evelo CT: User-friendly solutions for microarray quality control and pre-processing on ArrayAnalysis.org.Nucleic Acids Res 2013, 41:W71-W76.PubMed CentralPubMedView Article
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.