Correlation of microRNA levels during hypoxia with predicted target mRNAs through genome-wide microarray analysis
© Guimbellot et al. 2009
Received: 20 March 2008
Accepted: 25 March 2009
Published: 25 March 2009
Skip to main content
© Guimbellot et al. 2009
Received: 20 March 2008
Accepted: 25 March 2009
Published: 25 March 2009
Low levels of oxygen in tissues, seen in situations such as chronic lung disease, necrotic tumors, and high altitude exposures, initiate a signaling pathway that results in active transcription of genes possessing a hypoxia response element (HRE). The aim of this study was to investigate whether a change in miRNA expression following hypoxia could account for changes in the cellular transcriptome based on currently available miRNA target prediction tools.
To identify changes induced by hypoxia, we conducted mRNA- and miRNA-array-based experiments in HT29 cells, and performed comparative analysis of the resulting data sets based on multiple target prediction algorithms. To date, few studies have investigated an environmental perturbation for effects on genome-wide miRNA levels, or their consequent influence on mRNA output.
Comparison of miRNAs with predicted mRNA targets indicated a lower level of concordance than expected. We did, however, find preliminary evidence of combinatorial regulation of mRNA expression by miRNA.
Target prediction programs and expression profiling techniques do not yet adequately represent the complexity of miRNA-mediated gene repression, and new methods may be required to better elucidate these pathways. Our data suggest the physiologic impact of miRNAs on cellular transcription results from a multifaceted network of miRNA and mRNA relationships, working together in an interconnected system and in context of hundreds of RNA species. The methods described here for comparative analysis of cellular miRNA and mRNA will be useful for understanding genome wide regulatory responsiveness and refining miRNA predictive algorithms.
MicroRNAs (miRNA) are approximately 22-nucleotide, non-coding RNA sequences important in the control of gene expression. They are involved in a variety of cellular processes, including development, cell differentiation, signaling, and tumorigenesis, and are believed to represent 1% of the predicted genes in mammalian and nematode genomes[2, 3]. Mammals in general (and primates in particular) appear to have a large number of miRNAs not found in other animal orders, suggesting that many functional miRNAs may have emerged during recent evolutionary periods. According to current functional and predictive models, each miRNA regulates multiple genes during differentiation and/or development at the transcription, translation, and posttranslational levels[1, 4, 5]. However, few of these targets and regulatory pathways have been experimentally validated, and the number of authentic (as opposed to predicted) miRNAs that exist in the mammalian genome as well as the actual number of their targets are not yet known.
Considerable effort has been directed toward understanding which mRNAs within the human genome are subject to regulation by miRNA-mediated repression. The miRGen Targets interface allows users to search either for targets or particular miRNA(s) that influence a particular gene. DIANA-microT, MiRanda, TargetScanS, and PicTar are four genome-wide prediction algorithms whose results are available through miRGEN http://www.diana.pcbi.upenn.edu/miRGen.html[8–11], an integrated database of (i) positional relationships between animal miRNAs and genomic annotation sets, and (ii) animal miRNA targets according to combinations of widely used prediction programs. These algorithms can provide quite a variable picture of miRNA behavior, and it is difficult to assess which in silico predictive method is best for identifying true miRNA targets. It is probable that the use of multiple programs combined with mRNA expression profiling will be necessary to address this question. As a result, we considered four different algorithms (PicTar, TargetScanS, miRanda(microrna.org), and miRanda(miRBase)) in this report when assessing the relationship between mRNA and miRNA expression.
Previous studies have evaluated the influence of a particular miRNA on potential mRNA targets and found a high degree of correlation in specific tissues [13–15]. In general, these reports have examined the regulatory relationships between miRNAs and the genome wide transcriptome, with a focus on pathological conditions (such as cancer) rather than acute perturbations such as hypoxia, although hypoxia in particular has been shown to regulate discrete miRNAs [16–23]. Hypoxia results in a change in expression of a significant portion of the human transcriptome. After oxygen restriction, we observed down-regulation of hundreds of transcripts, including the cystic fibrosis transmembrane conductance regulator (CFTR), in which we could not identify a consensus hypoxia regulatory motif (HRE, hypoxia regulatory element; an indicator of transcriptional regulation by hypoxia inducible factor (HIF)). In the present study, we hypothesized that many of these transcripts may be down-regulated by miRNAs.
Kulshreshtha et al recently demonstrated a functional link between hypoxia and microRNA expression, although the relationship to mRNA expression was not evaluated. We therefore investigated the effects of hypoxia on a model epithelia and found that 3125 unique genes were significantly altered. Of these, approximately 53% were down-regulated, presenting 1649 unique possible targets for miRNA-mediated repression. Expression data from a miRNA Bioarray (see Methods) was compared with the mRNA expression profile, and the strength of correlation against predicted targets with differentially expressed miRNAs was analyzed using computational techniques we developed specifically for this purpose. We found no compelling evidence that miRNA-mediated repression plays a major role in down-regulation of CFTR, and present evidence that the individual miRNA levels do not correlate well with their algorithm-predicted target mRNAs. However, the groups of miRNAs predicted to regulate the same mRNA target were found to be co-regulated, indicating that a level of combinatorial control may exist.
The HT29 (human colonic) cell line was obtained from ATCC http://www.atcc.org and seeded on 12-mm diameter Transwell filters (Corning-Costar, Corning, NY). Cells were cultured in media (HT29: McCoy's 5a medium supplemented with 7% Fetal Bovine Serum (FBS)) for 5–7 days (media bathing both the apical and basolateral compartments) at 37°C (5% CO2 – 95% air gas mixture). Under these conditions, cells form polarized monolayers with transepithelial resistances of >1000 Ω·cm2. In some monolayers, media was removed from the apical side to expose cells to air (A/L or air-liquid interface). Media remained on the apical surface (L/L or liquid-liquid interface) of other monolayers to a depth of one centimeter, a condition that markedly impairs access to ambient oxygen, conferring a hypoxic environment at the cell surface and a glycolytic cytosol due to an impairment of oxygen diffusion through liquid [25–27]. Several studies have established that altered physiology of cells under submerged conditions is primarily the result of a lack of oxygen (with otherwise equal volumes or composition of overlying media) [27–29]. For example, it has been determined that cells cultured under liquid (L/L) are hypoxic as judged by three criteria; (i) HIF protein level is increased in L/L culture, (ii) HRE-driven luciferase activity is increased, and (iii) well known genes activated under hypoxia such as VEGF are dramatically elevated. The assay represents a useful and commonly used test for observing changes in vitro due to hypoxia.
Total RNA was purified from HT29 cells using the mirVana™ miRNA Isolation Kit per manufacturer instructions (Ambion, Inc., Austin, TX), and RNA quality assessed before RNA labeling (2100 Bioanalyzer, Agilent, Palo Alto, CA). Detailed analysis procedures are presented in the Manufacturer's GeneChip Expression Technical Manual (Affymetrix, Santa Clara, CA). Briefly, 2 ug of total RNA from each sample was used to generate double strand cDNA by linear amplification using oligo dT-T7 primer and reverse transcriptase. Subsequently, biotin-labeled cRNA was synthesized by in vitro transcription (IVT) using 3'-Amplification Reagents for IVT labeling (Affymetrix, Santa Clara, CA) followed by cRNA fragmentation. The Affymetrix Human Genome U133 Plus 2.0 Array was used for hybridization. This array contains 54675 probes designed to over 47,000 transcripts and variants. Arrays were hybridized overnight at 45°C, and then washed, stained, and scanned on a GeneChip Scanner 3000 (Affymetrix, Inc., Santa Clara, CA). Gene expression levels were analyzed with GeneChip Operating Software (Affymetrix, Inc., Santa Clara, CA). Raw data were analyzed using Microarray Suite, Version 5.0 software (Affymetrix, Santa Clara, CA). The raw data set is available through Gene Expression Omnibus under accession number GSE9234.
Total RNA was purified from HT29 cells using the mirVana™ miRNA Isolation Kit per manufacturer instructions (Ambion, Inc., Austin, TX) which efficiently purifies RNA as small as 10 nucleotides. Expression profiling was then performed using the mirVana miRNA Bioarrays V2 (Asuragen, Inc., Austin, Texas) which contains probes for all mouse, rat, and human miRNAs (266, 238, 482 confirmed miRNAs, respectively) in miRBase http://microrna.sanger.ac.uk/. Samples for microRNA profiling studies were processed by Asuragen, Inc. (Austin, TX) and the microRNA enriched fraction obtained by passing total RNA through a flash PAGE™ Fractionator apparatus (Ambion, Inc., Austin, TX) and cleaned. The 3' ends of RNA molecules were tailed and labeled using the mirVana™ miRNA Labeling Kit (Ambion, Inc., Austin, TX). Amine-modified nucleotides were incorporated during the poly (A) polymerase mediated tailing reaction, and Cy3 succinimide esters (Amersham Biosciences (GE Healthcare), Piscataway, NJ) were conjugated to amine moieties on microRNAs. Hybridization to the mirVana miRNA Bioarrays (Ambion, Inc., Austin, TX) was performed. The Cy3 fluorescence on the arrays was scanned at an excitation wavelength of 532 nm using a GenePix 4200AL scanner (Molecular Devices, Union City, CA). The fluorescent signal associated with the probes and local background was extracted using GenePix Pro (version 6.0, Molecular Devices, Union City, CA). The raw data set is available through Gene Expression Omnibus under accession number GSE9234.
The raw microarray data obtained from Microarray Suite v5.0 software were analyzed using a two-sided t-test corrected for unequal variances (Welch test) to compare the mean expression level for each gene between the two groups. A Bayesian posterior probability of being a false positive result (expressed as the false discovery rate, FDR) was estimated for each probe set individually, based on the Welch t-test p-values and using a mixture model[30, 31]. We focused on the genes among those most differentially expressed that had corresponding probe sets with a lower than 1% FDR, that is, with a posterior probability of being differentially expressed of 99%. These genes were annotated with their chromosomal locations, UniGene number, LocusLink ID and Gene Ontology (GO) information (genes were grouped according to biological process, cellular component, or molecular function) using NetAffx resources. We used Onto-Express (available at http://vortex.cs.wayne.edu/Projects.html, last accessed in March 2007)  to calculate whether any of the GO terms were significantly over-represented among differentially expressed genes, as determined by a two-sided binomial test. The p-value calculation is only valid if the expression levels of the genes are independent, which is probably not the case in expression studies; thus, the p-values reported for these analyses should only be considered as heuristic ranking statistics. The fold change represents the ratio between microarray probe expression values.
Thresholding and signal scaling were generated using algorithms selected by Asuragen. The background-adjusted fluorescent values generated by GenePix Pro were normalized for each microRNA using a variance stabilization method described by Huber et al, followed by a Welch two-sample t-test carried out for every gene; and a multiplicity correction was conducted to control FDR at 5% using a step-up approach, as described by Benjamini and Hochberg.
Representative list of genes changed between hypoxia and normoxia including number of miRNAs predicted to have target sites in each gene (p ≤ 0.001).
Number of miRNAs as predicted by miRBase
fructose-bisphosphate aldolase C
hypoxia-inducible protein 2
egl nine homolog 3 (C. elegans)
carbonic anhydrase IX
vascular endothelial growth factor
Early growth response 1
DEAD (Asp-Glu-Ala-Asp) box polypeptide 28
NADPH oxidase 1
cystic fibrosis transmembrane conductance regulator, ATP-binding cassette (sub-family C, member 7)
DEAD (Asp-Glu-Ala-Asp) box polypeptide 52
DEAD (Asp-Glu-Ala-Asp) box polypeptide 27
Most significant GO functions in three GO classes.
Most significant GO Biological Processes
Total Genes in Class
Nuclear mRNA splicing, via spliceosome
Regulation of the cyclin dependent protein kinase activity
Amino acid biosynthesis
Regulation of progression through cell cycle
Regulation of translational initiation
Most significant GO Molecular Function.
Total Genes in Class
Unfolded protein binding
L-ascorbic acid binding
ATP-dependent RNA helicase activity
Most significant GO Cellular Component
Total genes in class
mitochondrial inner membrane presequence translocase complex
DEAD and DEAH helicases down-regulated in response to hypoxia.
DEAD (Asp-Glu-Ala-Asp) box polypeptide 1
DEAD (Asp-Glu-Ala-Asp) box polypeptide 5
DEAD (Asp-Glu-Ala-Asp) box polypeptide 10
DEAD (Asp-Glu-Ala-Asp) box polypeptide 17
DEAD (Asp-Glu-Ala-Asp) box polypeptide 18
DEAD (Asp-Glu-Ala-As) box polypeptide 19A
DEAD (Asp-Glu-Ala-Asp) box polypeptide 21
DEAD (Asp-Glu-Ala-Asp) box polypeptide 27
DEAD (Asp-Glu-Ala-Asp) box polypeptide 28
DEAD (Asp-Glu-Ala-Asp) box polypeptide 31
DEAD (Asp-Glu-Ala-Asp) box polypeptide 39
DEAD (Asp-Glu-Ala-Asp) box polypeptide 42
DEAD (Asp-Glu-Ala-Asp) box polypeptide 52
DEAD (Asp-Glu-Ala-Asp) box polypeptide 55
DEAD (Asp-Glu-Ala-Asp) box polypeptide 56
DEAD (Asp-Glu-Ala-Asp) box polypeptide 58
DEAH (Asp-Glu-Ala-His) box polypeptide 29
DEAH (Asp-Glu-Ala-His) box polypeptide 33
DEAH (Asp-Glu-Ala-His) box polypeptide 36
DEAH (Asp-Glu-Ala-His) box polypeptide 9
MicroRNAs up-regulated by hypoxia (p ≤ 0.05).
(Hypoxia vs Normoxia)
MicroRNAs down-regulated by hypoxia (p ≤ 0.05).
(Hypoxia vs Normoxia)
Several mRNA changes were statistically significant due to pronounced differential expression of a single miRNA (for example Figure 10), a finding that may or may not be biologically relevant. In addition, while many mRNAs within a gene-specific miRNA group were found to be coordinately up-regulated, many others predicted to have target sites by miRanda(microrna.org) and other algorithms were not differentially expressed (Figure 8), suggesting that the definition of gene-specific miRNAs requires further refinement. Only a subset of genes with significantly regulated gene-specific miRNA groups were differentially expressed in our data set. The implication of these findings is that although miRNAs may be coordinately regulated, they do not predict expression changes of every predicted target.
Cellular responses to hypoxia can occur through stabilization of HIF, a well-established transcriptional activator, and result in enhanced expression of a variety of hypoxia-related genes. Much less is known regarding hypoxia-dependent transcriptional repression. Our mRNA array data indicate that a large number of transcripts are robustly downregulated following oxygen restriction in human epithelial cells. One goal of the present study was to investigate the extent to which changes in miRNA could account for these variations in the cellular transcriptome. We hypothesized that miRNAs would play a role suppressing certain genes during hypoxia, and tested this by comparing expression data from miRNA and mRNA expression profiles investigated in parallel.
Recent studies describe hypoxia specific miRNA signatures in a variety of cell types [16–23]. A functional link between hypoxia and miRNA expression has therefore been observed by others, although the relationship between mRNA and miRNA from a genome-wide perspective has not been investigated previously. Kulshreshtha and colleagues emphasized that a spectrum of miRNAs can be induced during hypoxia, and at least some of these occur via a HIF-dependent mechanism. Ten miRNAs reported previously as hypoxia-responsive were also identified in our experiments (e.g. miR-23a, -23b, -27b, -30d, -191, -210, -213, -155, -200a, -181b) using a different method of oxygen deprivation (Tables 4 and 5). Interestingly, three miRNAs (miR-155, -200a, -181b) reported to be upregulated in other cell systems (reviewed in [16–23]) were noted to be repressed in colonic epithelia. These differences most likely relate to the various cell-types, growth conditions, or procedural aspects used in earlier studies.
The epithelial model of hypoxia described here represents a well-defined in vitro system for studying subacute (including transcriptional) effects of oxygen restriction [25–29]. We used the model to evaluate miRNA regulation of gene expression. We found changes in the epithelial transcriptome resulting from low oxygen, as well as further evidence for a potential signature of miRNAs induced by hypoxia [16–23]. However, in contrast to several extant models [13–15], we did not observe a significant correlation between mRNA expression levels and miRNAs on a genome-wide scale. Earlier studies have relied primarily on particular tissue types and developmental stages from a variety of organisms, suggesting results most relevant to embryologic gene regulation. The present investigation of mRNA:miRNA association applied a novel analytical approach to widely available data visualization tools, and monitored miRNA and mRNA expression on a genome-wide basis, including the potential role of environmental stressors (found commonly in pathologic conditions) on miRNA-mediated regulation.
Our analysis incorporated four miRNA target predictions programs (MiRanda, PicTar, miRBase, TargetScan). When miRNA targets were compared to mRNA output, the data sets failed to indicate a significant relationship between expression arrays. In a very recent study, Baek et al.  reported that the top third TargetScan predictions (ranked by 'total context score') may correlate best with protein downregulation. In the present experiments, applying this stringent threshold and strict site conservation (after ) did not result in a stronger association. This included use of the most recent TargetScan algorithm (release 4.2; http://www.targetscan.org/) and restriction of targets to a context score of 85% or higher. Our results therefore indicate limitations of the currently available target prediction algorithms. While high stringency methods can be valuable for an individual miRNA , TargetScan/PicTar modifications do not appear to enhance the available algorithms in a broader, genomic context.
The lack of a significant and robust correspondence between mRNA levels and miRNA expression could represent a lack of specificity and/or accuracy of miRanda or other target prediction algorithms. The observed magnitude of miRNA expression changes (Figures 8b–d) in the present experiments is lower than observed for mRNA (Figure 1). In addition, the relatively small sample sizes used in this study could contribute to a lack of information, making it difficult to test the assumptions underlying the statistical method (such as normality), in a fashion that could impact results. MiRanda typically produces more potential targets than other programs, and a large number of false targets would seriously limit the computational methods described here. We also note that the available programs have only partially overlapping predicted targets for the same miRNA and produce smaller data sets than miRanda. Due to the differences among databases and because there are no clearly superior methods, future studies of mRNA and miRNA regulation should consider analysis of multiple predictive algorithms rather than use of a single data analysis tool.
Although miRNAs can act to promote cleavage and subsequent degradation of their mRNA targets, this may not be the only (or even primary) mechanism of miRNA action in mammalian epithelia. A strong consensus is not yet available regarding the predominant pathway(s) that underly miRNA gene repression [41–44]. One explanation for our findings could relate to translational repression as a major action of miRNA in human cells. It has been shown that certain miRNAs bind their targets and prevent adequate translation. However, mRNA levels are not always affected by this process. A quantitative, proteomic approach to evaluate hypoxic protein expression in epithelia followed by in silico statistical correlation would be necessary to investigate this possibility. On the other hand, miRNA levels are also governed by DNA promoter elements, stability of miRNA, degradative pathways related to differential RNA editing, transport into the cytoplasm, and/or deficient processing by Drosha. Alternative transcript splicing and polyadenylation can eliminate miRNA regulatory sites from a message, and miRNA directed repression can be blocked by certain RNA binding proteins. It seems less likely that common promoter element(s) or a single pathway (by itself) could explain the very large number of up- or down-regulated miRNAs noted as a result of oxygen restriction (Figure 1). Moreover, translation of miRNA targets leads to secondary transcriptional and post-transcriptional regulation that contributes to the observed mRNA profile. The diversity of potential regulatory sequences, difficulty predicting biologic regulation based solely on a consensus miRNA binding site, and the increasingly apparent need for confirmation in living cells indicate that additional, cell-based studies should be used in the future to address transcriptome regulation by miRNA.
miRNA expression arrays represent a relatively new technology, and potential issues exist with regard to data acquisition. The correlation of biological replicates in our studies was >0.99, which indicates the technology is precise, although accuracy is undefined. In addition, we randomized the order of miRNA and mRNA extraction to minimize non-biological, confounding variables. The goal of identifing a method to predict levels of mRNAs based on miRNA profiling, regardless of the underlying regulatory mechanism, was strengthened by correlation against predicted mRNA targets across the entire transcriptome. While previous studies have evaluated effects of a single miRNA after high level recombinant overexpression, the present experiments allowed us to study the dynamics of miRNA and mRNA regulation in parallel with a common physiologic insult (oxygen deprivation). This approach avoided potential variables introduced by overexpression of foreign DNA elements or otherwise manipulating the cellular genome.
The present findings suggest that correlation between miRNAs and their predicted targets based primarily on the number of consensus sites in the 3'UTR may be overly simplistic. Combinatorial analysis reveals much more significant agreement between specific genes and their predicted miRNA regulators as a group; however, this too may reflect a one-dimensional view of miRNA activity. Based on evidence presented here that entire (GO) functional categories of mRNAs are regulated in parallel by hypoxia (Table 2), higher order miRNA groupings may exist along functional or developmental lines that respond as networks. In either case, the present experiments provide a means by which other predicted target lists – either currently available or under development – may be optimized to yield a better correlation between miRNA levels and gene expression.
The observation that a gene-specific group of miRNAs may work in concert to repress CFTR mRNA during hypoxia also points to a novel mechanism of regulation. Previous experiments have failed to establish a direct role for HIF during the pronounced inhibition of CFTR that occurs during oxygen deprivation. Moreover, very few gene products are believed to be down-regulated in a direct fashion by HIF. If a cohort of miRNAs target CFTR and coordinately suppress its message, this could represent an important and novel example of miRNA based repression following an environmental stress. The findings may also help explain in vivo suppression of CFTR mRNA during low oxygen exposure[45, 46], and suggest a role for miRNAs governing levels of hundreds of gene products following hypoxic insult (Figures 2 and 8).
Our results suggest that the expected inverse relationship between miRNA and target mRNA may be a rare event. Several previous studies [13–15] have indicated a clear correlation between a specific miRNA and suppression of a target mRNA. These earlier studies in some cases were based on marked overexpression of a particular miRNA, followed by expression studies of the mRNAs of interest. However, our experiments suggest that under physiological conditions in human epithelium, miRNA acts in a more subtle fashion distinct from that of marked overexpression. In addition, the physiologic impact of miRNAs on cellular transcription appears to result from a multifaceted network of miRNA and mRNA relationships, working together in an interconnected system and in context of hundreds of other RNA species. It may be that target prediction algorithms and expression profiling techniques do not yet adequately represent the complexity of miRNA-mediated gene repression, and new methods may be required to truly understand these systemic aspects.
This work was supported by Cystic Fibrosis Foundation (R464) and U.S. National Institutes of Health (P30DK072482, F30ES014987, U54 AT100949, 5T32HL072757-04 to SWE).
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 cited.