The HT29 (human colonic) cell line was obtained from ATCC 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% CO₂ – 95% air gas mixture). Under these conditions, cells form polarized monolayers with transepithelial resistances of >1000 Ω·cm². 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 . 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[32]. We used Onto-Express (available at , last accessed in March 2007) [33] 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[34], 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[35].

Data analysis was performed using three replicates of array data per group (within group Mean r² = 0.996, range 0.004–0.998). Genes were considered differentially expressed between normoxic and hypoxic conditions if absolute expression changes were 1.5-fold or greater, a lower than 0.1% false discovery rate (FDR) was observed, and the p-value from a Welch two-sample t-test was less than 0.001. Over 8% of all probes satisfied these criteria (1999 probes up-regulated and 2099 down-regulated, see Figure 1). These probe sets map on to 1476 up- and 1649 down-regulated genes with a HUGO gene symbol, along with 208 probe sets that did not map to a HUGO identifier. Table 1 provides examples of selected genes known to be regulated by hypoxia with a broad biological significance, which are also predicted to contain target sites (in the 3'UTR) for miRNAs, according to the miRanda(microrna.org) target prediction database (for complete dataset of all genes altered by oxygen restriction in our studies, see GEO series GSE9234). Interestingly, the miRanda(mirorna.org) Target Database predicted that every 3'UTR in this set had target sites for between 4 and 69 miRNAs, although only 30% of the genome has been suggested to be regulated by miRNA-based mechanisms[1,4,5].

We next analyzed our mRNA array data to identify various cellular processes affected by hypoxia according to Gene Ontology (GO) annotation[36]. GO is organized into three partially overlapping categories that consider three different aspects of each gene: biological processes, molecular functions, and cellular component. To investigate whether any GO terms were significantly over-represented among the differentially expressed genes, we used the Onto-Express tool[33] to calculate statistical significance values for each category. Categories with the most significant corrected p-value are shown in Table 2. The highest fraction of differentially expressed genes among the biological processes class included genes linked to cholesterol biosynthesis, protein metabolism, or ribosome biogenesis. Among the molecular function class, highly represented genes included those involved in L-ascorbic acid binding, ATP-dependent RNA helicase activity, oxidoreductase activity, NAD binding, and tRNA binding. A high fraction of mitochondrial inner membrane presequence translocase complex proteins was found among the cellular components class. Only a small portion of GO terms in each category changed significantly (less than 1% of entire terms, p < 0.05), suggesting the hypoxia-response may be highly specialized. Taken together, all three categories indicate considerable effects on glycolysis, translation and protein metabolism, and RNA processing. Of note was one group of proteins particularly affected by hypoxia. At least 20 members of the DEAD box family were significantly down-regulated, 2 to 33 fold, and comprised individual gene products in molecular function and biological process categories important in cellular pathways involving RNA (Table 3) [37,38]. These results suggest a substantial impact on RNA metabolism and activity by hypoxia.

Unlike well-characterized hypoxia-mediated transcriptional activation[24], the general mechanisms underlying gene repression due to hypoxia are not well understood. We hypothesized that miRNAs may play an important role in the down-regulation of gene expression by hypoxia. We addressed this hypothesis by combining a comprehensive mRNA expression array with miRNA bioarray to evaluate not only the potential for a specific, hypoxic stress-induced miRNA profile, but also to correlate the expression of specific miRNAs with their algorithm-predicted targets. The mirVana miRNA Bioarrays V2(Ambion, Inc., TX) used in this study contain probes for mouse, rat, and human miRNAs in miRBase. A hierarchical cluster analysis, using average linkage and Pearson's correlation as the weight function, of all the significant miRNAs evaluated in three hypoxic and three normoxic conditions is shown in Figure 2, illustrative of miRNAs differentially expressed by hypoxia. Approximately 8% (53 miRNAs) of all human miRNAs tested (640 unique miRNAs) were significantly different in the hypoxic group compared to the normoxic group. Twenty-eight miRNAs (4%) were up-regulated in a statistically significant fashion (p < 0.05) by hypoxia and twenty-five miRNAs (4%) were significantly down-regulated (p < 0.05) (Tables 4 and 5). No correlation of target mRNAs as predicted by miRanda (microrna.org) with any particular GO term was identified, supporting previous reports that miRNAs regulate a wide variety of mRNAs and their action is not restricted to specific groups of genes[5,39].

While it is well established that miRNAs play an important role governing gene expression, the correlation of mRNA and miRNA profiles across the entire genome and in response to an acute cellular perturbation (such as hypoxia) have not been thoroughly assessed. Accordingly, we next evaluated the extent to which the profile of miRNAs observed during hypoxia was consonant with mRNA expression (i.e., the transcriptome response). The distribution of t-statistics for all probes on mRNA expression arrays was used as a reference distribution (dotted line), as shown in Figure 3. In comparison to the reference distribution, the thirteen most significant differentially expressed miRNAs (Tables 4 and 5) and twenty randomly selected non-significant miRNAs were analyzed. The distribution of t-statistics for mRNA probes that are predicted targets of the given miRNAs were generated (solid lines, Figures 3a and 3b). Figure 3 depicts the distribution of mRNAs with one predicted target site for a given miRNA (microrna.org); a deviation from the reference curve (dotted lines) would be interpreted as a significant relationship between the mRNA and miRNA arrays, as mediated via the target predictions. This methodological approach to investigating 'cause and effect' allows a visualization of mRNA/miRNA correlations. No consistent deviations, however, were observed in the combined t-statistics. As a further test of this finding, evaluations of individual miRNAs were performed (representative examples in Figures 3c and 3d). No consistent deviations were noted in t-statistics for any single miRNA in either the significant or the non-significant miRNA group. Figure 3 illustrates the distribution for miRanda-predicted targets, while Figure 4 through Figure 6 contain the results for PicTar, TargetScanS, and miRBase, respectively. When all miRNAs that changed significantly during oxygen restriction were compared against mRNAs containing one target prediction site in the 3'UTR, none of the target prediction programs indicated a significant relationship between data from the mRNA and miRNA expression arrays. When the findings were evaluated using the most recent TargetScanS predictive method (i.e. a more stringent threshold [40] that includes only the highest ranked targets with a context score percentile of 85%), no consistent deviations were observed in the combined t-statistics (Figure 5).

Similar comparisons were made for mRNAs containing at least three target sites within each 3'UTR for the same miRNA (Figure 7). As in Figure 3, the 13 most significant differentially expressed miRNAs (Figure 7a) and 20 randomly selected, non-significant miRNAs (Figure 7b) were investigated. Certain miRNAs did exhibit some relationship with their predicted targets (indicated by arrows); however, miRNAs without statistically significant differential expression also demonstrated a roughly comparable relationship (arrows, Figures 7c and 7d), suggesting that this observation occurred by chance.

Combinatorial regulation by groups of diverse miRNA species binding to different numbers of target sites within a single 3' UTR of a given gene has been hypothesized as a mechanism underlying miRNA-mediated gene repression[4]. Therefore, we next compared expression under hypoxia or normoxia of all miRNAs predicted to target a specific gene (termed a gene-specific miRNA group, Figure 8). This comparison evaluates whether specific groups of miRNAs (namely those predicted to target a specific gene) are significantly up- or down-regulated as a cluster. For each gene, a regression line (solid line, estimated best fit using the specific group of miRNAs indicated by red dots) was plotted and an ANOVA test performed to determine whether the regression line differed significantly from equality (dotted line). Figure 8a shows a histogram of the resulting p-values, indicating that a higher-than-expected number (expected result shown by dotted line) of gene-specific miRNA groups were regulated coordinately in association with a particular mRNA target.

A representative sample of statistically significant gene-specific miRNA groups is shown in Figures 8b–d. We chose to analyze three genes; CFTR, KIAA2026, and C16orf73. Selection of CFTR was based on our laboratory interest in regulation of that particular gene product. The other two genes shown in the figure (whose functions are not known) were selected because of highly significant p-values, indicating very strong co-regulation by gene-specific miRNA groups. Scatterplots were identical by this method (black dots representing all miRNAs), indicating the group of miRNAs predicted to regulate expression of a particular gene of interest. Significant deviation from the line of equality (dotted line) indicates co-regulation of a cohort of miRNAs. The results presented here therefore provide some of the first evidence to suggest that miRNAs may be coordinately regulated in groups relevant to specific 3' UTRs. For example, Figure 8b shows predicted regulation of CFTR mRNA expression by CFTR specific miRNAs (red dots, 9 predicted by miRanda(microrna.org)). The regression line indicates considerable deviation, suggesting that coordinated regulation of CFTR mRNA may exist by these specific miRNAs. CFTR mRNA levels were among the most significantly decreased under hypoxia based on the mRNA array (Table 1). Figure 8 shows results for miRanda-predicted targets. Results for PicTar, TargetScanS, and miRanda(miRBase) are provided in the Figures 9, 10, 11. All four algorithms indicate gene-specific groups of miRNAs reacted coordinately to hypoxia, as shown by the histograms of p-values. However, different algorithms identified somewhat different gene-specific groups as being most strongly affected by hypoxia, as shown by the scatterplots (Figures 8b–d; and Figures 9, 10, 11).

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.