Profiling the population of microRNAs in hypoxic MCF-7 cells by next-generation sequencing

The total number of sequencing reads obtained was in the range between 32 and 35 million reads per sample. After filtering, sequencing reads with a length between 18 and 26 nucleotides represented 75.2-87.5% of total number of reads obtained initially (Additional file 1: Table S1). On average, 99% of these reads (98.9-99.5%) were mapped to the human genome using Bowtie, from which 75-85.9% had a unique mapping position (Additional file 1: Table S1). Most of the reads mapped to annotated regions (98-98.9%), and as expected the vast majority of them matched genomic regions allocating microRNAs (94.5-96.1% of reads mapped to annotated regions). The residual fraction of mapped reads not corresponding to microRNAs was distributed among a miscellanea of annotated regions including other RNA and small RNAs species, protein coding exons, pseudogenes and repeat molecules (Additional file 1: Table S1).

We were able to identify 502 microRNAs consistently across all conditions in MCF-7 cells (Additional file 1: Table S2). Among them, 412 microRNAs had strand (-5p/-3p) annotations with 53.2% of 5p and 46.8% of 3p-microRNAs. We found that 184 microRNAs were from only one strand (either-5p or-3p), so there were no detectable traces of the other hairpin arm and 114 microRNAs have both the 5p and 3p strands expressed.

The sequencing data also allowed us to identify precursors on different genomic loci that generate the same mature microRNA. For six mature microRNAs (hsa-let-7f-5p, hsa-miR-16-5p, hsa-miR-181a-5p, hsa-miR-24-3p, hsa-miR-29b-3p and hsa-miR-30c-5p) we were able to identify microRNAs from the other strand derived from the two different precursors. These findings provide evidence that both loci of these particular microRNAs would be expressed in MCF-7 cells. Compared to a similar sequencing based study performed using human umbilical vein endothelial cells (HUVEC) [17], there was an overlap of 80 microRNA sequences.

Differential expression of microRNAs under hypoxia in MCF-7 cells

When looking at differential expression between hypoxia and normoxia, we found that the number of differentially expressed microRNAs increased along the time of exposure to hypoxia. We identified 1, 49 and 69 microRNAs up-regulated at 16 h, 32 h and 48 h of hypoxia respectively (adj.p-val ≤ 0.05) whereas 4, 37 and 50 microRNAs were detected as down-regulated at 16 h, 32 h and 48 h of hypoxia respectively (adj.p-val ≤ 0.05) (see Figure 1A). Among each group of significant microRNAs we found microRNAs from both 5p and 3p strands.

We have previously demonstrated the up-regulation of hsa-miR-210-3p (previously annotated as hsa-miR-210) across a hypoxia time course in MCF-7 cells [27], but this in depth study has revealed a number of hypoxia regulated microRNAs that have not been reported in previous studies [11–20] (highlighted in bold in Tables 1 and 2).

Agilent microRNA microarrays (3 biological replicates per condition) were also used to complement and validate the sequencing data in a broad manner, as the two technologies have their own strengths and weaknesses in microRNA profiling [28]. We found 228 microRNAs commonly detected between sequencing and microarrays (Additional file 1: Table S3), with a significant correlation between the number of reads and array intensity (Spearman: 0.57 (average), p-val ≤ 0.05). When comparing the fold changes obtained for each microRNA across the time course in each platform there is overall a significant correlation (Pearson: 0.30 (average), p-val ≤ 0.05). Most of the compared microRNAs show a clear trend towards the same direction in both platforms despite some punctual disagreements (see Additional file 2: Figure S1). Moreover around half of the microRNAs found significantly regulated through next-generation sequencing are also found significantly regulated in the microarray experiment (Additional file 2: Figure S1 and Additional file 1: Table S3).

In order to further validate the sequencing results obtained here, we performed qPCR to test the expression levels of hsa-miR-4521, hsa-miR-145-3p and hsa-miR-222-5p. Indeed, all microRNAs were validated as down-regulated in agreement with the sequencing data (Figure 1C, adj.p-val ≤ 0.05). The expression levels of hsa-miR-145-5p and hsa-miR-222-3p were found either not changed or slightly up by qPCR (Figure 1C). This suggests that one strand for these microRNAs is specifically down-regulated under hypoxic conditions and may have been functional as compared to the other strand.

Role of HIF in microRNA regulation under hypoxia

Direct regulation of microRNA expression by HIF-1α and HIF-2α in MCF-7 was assessed by ChIP-seq, while stabilization of the two proteins in the same cells under hypoxia has been published previously [25, 26]. We used the 500 published high-stringency HIF binding sites [29] in order to look for microRNAs that map to within 50 kb of each binding site. Gene Set Enrichment Analysis (GSEA) against the fold change by hypoxia at each of the three time points shows a weak enrichment for HIF-binding microRNAs amongst those up-regulated by hypoxia, which becomes more apparent at the longer time points (Figure 2A-C). The heat map (Figure 2D) shows the fold change at each time point for each microRNA locus containing a HIF binding site. The most up-regulated microRNA containing a HIF-binding site is indeed hsa-miR-210-3p, but there are clearly other microRNAs (Additional file 1: Table S4). In Figure 2E and F we show the tracks for two loci of particular interest, MIR30D and MIR27A, which are near the top of the list. Indeed, hsa-miR-30d-5p, hsa-miR-30d-3p and hsa-miR-27a-5p are significantly up-regulated at 32 h and 48 h of hypoxia according to our sequencing data. Both microRNA loci have un-annotated, hypoxia inducible host transcripts that appear in directional ribo-depleted RNA-seq data produced by exposing the MCF-7 cells to hypoxia (24 h 1% hypoxia) [30] and strong HIF-binding ChIP-seq peaks (Figure 2E and F).

Hsa-miR-27a is within a close distance to other two microRNAs, hsa-miR-23a and hsa-miR-24-2, whereas hsa-miR-30d is also very close to hsa-miR-30b. They are reported to be clustered microRNAs that seem to be transcribed in the same pri-miRNA, which agrees with our previous observations concerning the host genes for these microRNAs. Therefore, we checked the expression of these microRNAs by qPCR, validating the up-regulations seen in sequencing results for hsa-miR-30d-5p, hsa-miR-30d-3p and hsa-miR-27a-5p (Figure 3A and B, adj.p-val ≤ 0.05). In addition we found hsa-miR-27a-3p significantly up-regulated across the whole hypoxia time course, hsa-miR-23a-3p up-regulated at 48 h and hsa-miR-24-2-5p, hsa-miR-30b-5p and hsa-miR-30b-3p significantly up-regulated at 32 and 48 h (Figure 3A and B, adj-p-val ≤ 0.05). It is worth to mention that hsa-miR-27a-3p was also found significantly up-regulated across all hypoxia time course in Agilent microarray data (Additional file 1: Table S3). All these results would further support the role of HIF on regulating the expression of these microRNA clusters under hypoxic conditions. The results link genome wide HIF occupancy under hypoxic conditions to the up-regulation of pri-miRNAs and mature microRNAs.

Association of hypoxia regulated microRNAs with a hypoxia gene signature in breast cancer samples

microRNAs encoded within host genes

Some of the microRNAs found up- or down-regulated under hypoxia are located within intronic regions of protein coding genes. Therefore, we generated mRNA expression profiles of the same time course samples using Illumina arrays. If only considering microRNAs with a unique genomic position, we found 185 microRNAs encoded within 120 host genes that were detected in our microarray experiment (Additional file 1: Table S5). No correlation was found between microRNA expression and host gene expression at any hypoxia time point (Spearman rank correlation: 0.137 (P.val = 0.0613) at 16 h; -0.0949 (P.val = 0.1961) at 32 h; -0.0127 (P.val = 0.8621) at 48 h). Looking at the expression of genes containing microRNAs, we observed that they do not show a lot of variation. When dividing them in 3 groups according to whether they were significantly up-regulated, significantly down-regulated or without changes in hypoxia, the average fold-change of the up- or down-regulated groups is not very different to the group without changes (Figure 4A), suggesting that these genes are very lightly regulated by hypoxia at transcriptional level. For each of the group of genes, we represented boxplots with the expression of the microRNAs they contained (Figure 4B), showing the lack of correlation between the average fold-change of the microRNAs and the average fold-change of the corresponding host genes. Similarly from the microRNA point of view, we divided the microRNAs located within a host gene in three groups depending on their fold-change expression at 32 h and 48 h of hypoxia (up-regulated, down-regulated or no change, Figure 4C). Then for each group of microRNAs, we represented the fold-changes of the corresponding host genes (Figure 4D). While up-regulated and down-regulated microRNAs showed an average fold-change of 2 and 0.5 respectively (Figure 4C), the average fold-change of the corresponding host genes was very similar and close to 1 for both groups (Figure 4D). The picture is the same when Agilent microarray data are used to determine microRNA expression fold changes instead of small RNA-seq (Additional file 3: Figure S2). Detailed examples of these relationships can be seen in Additional file 1: Table S6. In general we observed around a 50% agreement between microRNA and host gene expression. In order to gain more information, we focused on cases where there is no correlation between microRNA and host gene expression and performed additional validations. The first example is hsa-miR-942-5p, found significantly up-regulated at 32 and 48 h, whereas its hosting gene TTF2 was significantly down-regulated in array data (Additional file 1: Tables S5 and S6). Indeed TTF2 was found significantly down-regulated in hypoxia by qPCR (adj.p-val ≤ 0.05, data not shown). Also hsa-miR-3140-3p is up-regulated in sequencing data whereas the host gene FBXW7 is down-regulated in array data (Additional file 1: Tables S5 and S6). FBXW7 has four isoforms and only two would contain MIR3140 within its sequence (Figure 5A). We designed a pair of primers (Figure 5A, F1 and R1) that would amplify three isoforms, including two that do not include the MIR3140 locus and a pair of primers specific for the two isoforms containing MIR3140 (Figure 5A, F2 and R2). All products were found significantly down-regulated in hypoxia by qPCR (Figure 5A, adj.p-val ≤ 0.05).

Another interesting example is the cluster composed by hsa-miR-106b and hsa-miR-93 and their host gene MCM7. The microRNA sequencing data showed a significant up-regulation of hsa-miR-106b-3p and down-regulation of hsa-miR-93-3p at 32 h and 48 h of hypoxia, whereas MCM7 was significantly down-regulated according to array data (Additional file 1: Tables S5 and S6). MCM7 has two isoforms and was confirmed as significantly down-regulated in hypoxia by qPCR when primers able to amplify both isoforms were used (Figure 5B, primers F0 and R0). However when using specific primers for testing the expression of each isoform separately, we found one isoform significantly down-regulated at all hypoxia times (Figure 5B, primers F1 and R1, adj.p-val ≤ 0.05) whereas the other isoform was significantly up-regulated at 32 and 48 h (Figure 5B, primers F2 and R2, adj.p-val ≤ 0.05). This isoform specific regulation would match the regulation of the microRNAs provided they are processed from preferred isoforms.

Changes in expression of microRNA processing genes due to hypoxia

Overall there is a certain degree of coordination between host gene and microRNA expression under hypoxia but also discrepancy. This suggests that other factors may be involved in the changes observed in mature microRNA expression under hypoxia. We used the gene expression profiles to specifically monitor any change in the expression of genes coding for proteins involved in the microRNA maturation pathway in hypoxic MCF-7 cells. Indeed, significant changes (see Additional file 1: Table S7, adj.p-val < 0.05) were detected in the expression of several genes, including DDX5 (encoding P68 which forms part of the protein complex lead by DROSHA), XPO5 and RAN (encoding EXPORTIN5 and RAN, respectively, both involved in microRNA transport from the nucleus to the cytoplasm), DICER and EIF2C2 and EIF2C4 (encoding AGO2 and AGO4, respectively, both associated to RISC complex). We also found that some of these changes were HIF-1α and HIF-2α dependent based on previous data, in which gene expression regulations dependent on these transcription factors were studied through microarray profiles after short interfering (si) RNA-based suppression of HIF subunits [26]. In order to validate these results, we performed qPCR using the MCF-7 hypoxia time course, and two HIF-1α and HIF-2α siRNA sample sets (Figure 6). Indeed, we could confirm a significant down-regulation of DDX5, XPO5, RAN, EIF2C2 and DICER consistently across the hypoxia time course in MCF-7 cells, as well as a significant up-regulation of EIF2C4 (Figure 6). We also confirmed by qPCR a significant up-regulation of both genes DICER and EIF2C2 when HIF-1α and HIF-2α are silenced in MCF-7 compared to the control, in both independent sample sets (Figure 6). Concerning EIF2C4, its expression was found significantly down-regulated upon HIF-1α silencing (Figure 6). Therefore the down-regulation of DICER and EIF2C2 under hypoxia is dependent on both HIF-1α and HIF-2α, whereas the up-regulation of EIF2C4 seems to be only dependent on HIF-1α. Changes observed in DDX5, XPO5 and RAN in hypoxic MCF-7 seem to be independent of HIF-1α or HIF-2α, as their expression did not change in the HIF siRNA experiments (Figure 6).

Despite the fact that some of the changes observed under hypoxia for these genes are HIF dependent, we could not find any high stringency HIF binding site close to any of these genes in our datasets. Since the down-regulation of DICER in hypoxic MCF-7 cells is the most strongly associated to HIF among the genes involved in microRNA processing, we decided to test whether this change was also observed in vivo. We used the gene expression data generated for the same breast cancer series we have referred in the previous section (GSE22220) [32] and investigated the correlation of DICER expression with the hypoxia score based on the expression of 99 genes (univariate analysis) [31]. Interestingly, we found that the expression of DICER was significantly inverse correlated to hypoxia in breast cancer (univariate analysis: rho = -0.37, pval = 3.e-08), meaning that the expression of DICER is lower as hypoxia increases. Therefore we show that DICER is also down-regulated under hypoxia in vivo in breast cancer and that this regulation is possibly dependent on HIF although not in a direct manner.

Cell culture and hypoxia exposure

Breast adenocarcinoma cell line MCF-7 cells were grown in Dulbecco’s Modified Eagle Medium supplemented with 10% fetal bovine serum (Sigma-Aldrich, St Loius, Missouri, USA) and 2 mmol/L of L-glutamine. Exposure of cell cultures to hypoxia (1% Oxygen) was undertaken in an In vivo2 Hypoxia Work Station (Ruskin Technologies, Pencoed, UK) in parallel with cells maintained in normoxic conditions (21% oxygen). Cells were recovered after 16 h, 32 h and 48 h of exposure to hypoxia. All experiments were done in triplicate from independent cell cultures.

RNA extraction

RNA was extracted from cells using the miRVana miRNA Isolation Kit (Ambion, Life Technologies, Carlsbad, California, USA). RNA quality and abundance were determined after extraction using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA) and a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA), respectively.

Preparation of microRNA libraries for next-generation sequencing

For each experimental condition (normoxia, 16 h hypoxia, 32 h hypoxia, 48 h hypoxia), small RNA libraries for next-generation sequencing were prepared from RNA samples obtained from two biological replicates. Each library was prepared from the 20 to 30 nucleotide RNA fraction. This was isolated from 2 μg of total RNA after being run in a 15% urea-TBE gel (Invitrogen, Life technologies, Carlsbad, California, USA) for 1 h at 200 V. The RNA contained in the excised band was eluted in 300 μl of 0.3 M NaCl by 4 h of constant rotation at room temperature. The elute was separated from the gel debris using a Spin-X-column (Thermo Fisher Scientific, Waltham, Massachusetts, USA) and RNA was precipitated by adding 750 μl of 100% ethanol and 3 μl of glycogen (1 mg/ml) and incubating for 30 minutes at -80°C. The precipitated RNA was centrifuged at 14,000 rpm for 25 minutes at 4°C, washed with 75% ethanol, air dried and resuspended in 5 μl of RNAse-free water. This material was used to prepare sequencing libraries using the DGE-small RNA sample preparation kit (v1.5) from Illumina (Illumina inc., San Diego, CA, USA ) according to the protocol provided by the manufacturer. The abundance of the libraries was measured using Qubit dsDNA HS Assay kit and a Qubit 2.0 fluorometer (both from Life technologies, Carlsbad, California, USA). The size of the products contained in the libraries was assessed with a DNA 1000 chip and an Agilent 2100 Bioanalyser (both from Agilent Technologies, Santa Clara, California, USA). The libraries were sequenced using an Illumina Hiseq sequencer according to the manufacture’s standard protocols at 36 bp read length.

Alignment and analysis of microRNA next-generation sequencing data

Small RNA sequencing data were processed from raw FASTQ files and analysed using the Kraken pipeline (http://www.ebi.ac.uk/research/enright/software/kraken) [49]. The 3’ adaptors were stripped using the criteria: 12-nt alignment stretch with no more than 2 mismatches and no gaps. The reads were then filtered for low-complexity regions and the remaining reads were size-selected for 18-26nt, resulting in inclusion of more than 75% of reads on average for each sample (Additional file 1: Table S1). A file with unique reads and their corresponding counts was generated for each sample. The final processed unique reads were mapped to the human genome (hg19 based on Ensembl v64 [50]) using Bowtie 0.12.7 [51] allowing for two mismatches and using the best alignment stratum option. Reads mapped to more than 20 loci were discarded. Mapped reads were classified and counted based on genomic annotations (RNAs, coding genes, pseudogenes and repeats) obtained from Ensembl v64 (Additional file 1: Table S1). For reads mapping to multiple loci, counts were divided by total number of reads. The sense and antisense distribution of mapped reads across all chromosomes was also analysed.

For differential expression of microRNAs, the overlap between mapped reads and human mature microRNAs (based on miRBase v17 [52]) was found using the function findOverlaps available in the Bioconductor GenomicRanges package (http://bioconductor.org/packages/2.6/bioc/html/GenomicRanges.html). The normalization and differential expression was performed with Bioconductor edgeR package v2.3.57 [53] using both common and tagwise dispersion. The significant differentially expressed microRNAs were determined by an adjusted p value lower than 0.05 based on the Benjamini and Hochberg multiple testing correction [54]. Data has been submitted to GEO (reference series: GSE47534; miRNA-seq: GSE47602).

HIF-1α and HIF-2α ChIP-seq analysis

We used the 500 published high-stringency HIF binding sites [29] in order to look for microRNAs that map to within 50 kb of each binding site. Gene Set Enrichment Analysis (GSEA) was carried out as described previously [29].

microRNA microarrays

microRNA expression profiles were performed using the human microRNA microarray V2 (Agilent Technologies, Santa Clara, California, USA) which contains 723 human microRNAs (Sanger database V.10.1) and 76 human viral microRNAs. Three biological replicate samples for each experimental condition (normoxia, 16 h hypoxia, 32 h hypoxia and 48 h hypoxia) were analysed. Samples were prepared with the Mirna Labelling Reagent and Hybridization kit (Agilent Technologies, Santa Clara, California, USA) following the instructions supplied by the manufacturer. The microarray hybridization was performed at 55°C for 20 h and at a constant speed rotation of 20 rpm. The hybridization signals were detected by the Agilent Microarray scanner G2565BA and quantified using the Agilent Feature extraction software version 9.5.1 (both from Agilent Technologies, Santa Clara, California, USA). A GeneView file including the Total Gene Signal (raw data) and the Total Gene Error for each of the microRNAs interrogated on the microarray was generated for all the samples and used for further analysis. One sample corresponding to MCF-7 cells exposed to hypoxia for 16 h was affected by a random technical problem and had to be removed due to the bad quality of the scan.

microRNA microarray data analysis

The Total Gene Signal was normalised using a variance stabilization method described by Huber and colleagues [55]. Furthermore a detection factor was calculated for each microRNA in every sample as the following subtraction: Total Gene Signal - Total Gene Error. A microRNA was considered detected when the value of this factor was higher than zero. The normalised data was then filtered based on the detection factor and on the microRNA class. The minimal requirements for a microRNA to be selected for further analysis were: (i) detection in at least two biological replicates of any experimental condition and (ii) be included as human microRNA in miRBase (http://www.mirbase.org). Limma analysis was then performed to assess independently the effects of 16 h, 32 h and 48 h of hypoxia exposure compared to normoxic control in microRNA expression levels in MCF-7 cells. Benjamini and Hochberg method was used to correct for multiple testing [54]. microRNAs with adjusted p values lower than 0.05 were considered significantly regulated. Data has been submitted to GEO (reference series: GSE47534; miRNA: GSE47532).

Gene expression microarrays

mRNA expression was measured using Illumina HumanWG-6 v3.0 Expression Bead Chip arrays (Illumina inc., San Diego, CA, USA). The Illumina platform was chosen as it performed well in large comparative studies [56]. Three biological replicate samples for each experimental condition (normoxia, 16 h hypoxia, 32 h hypoxia and 48 h hypoxia) were analysed. RNA was amplified using the Illumina Total Prep RNA Amplification Kit (Ambion, Life Technologies, Carlsbad, California, USA). Amplified RNA product (850 ng) was hybridised to the Illumina microarrays using single chamber hybridization cartridges. Washing, staining and scanning were carried out as specified in the Illumina Whole Genome Expression Manual v.1.

Gene expression microarray data analysis

Average signal was background subtracted with local background subtraction (BeadStudio) and consolidated per gene. Data was quantile normalized in Bioconductor (http://www.bioconductor.org) and filtered based on gene detection level using the pvalDet parameter. Only genes with pvalDet lower than 0.05 in at least 2 replicates of any of the experimental condition were considered for further analysis. Limma analysis was then performed to assess independently the effects of 16 h, 32 h and 48 h of hypoxia exposure compared to normoxic control in gene expression levels in MCF-7 cells. Benjamini and Hochberg method was used to correct for multiple testing [54]. Genes with adjusted p values lower than 0.05 were considered significantly regulated. Data has been submitted to GEO (reference series: GSE47534; mRNA: GSE47533).

Datasets from published studies (presented in Additional file 1: Table S10) have been analysed as described above.

Real time PCR (qPCR) analysis

MicroRNA expression was assessed by qPCR with TaqMan microRNA assay protocol (Applied Biosystems, Life Technologies, Carlsbad, California, USA) using 5 ng of total RNA per microRNA as indicated by the manufacturer. QPCR was done in a CFX96 real-time PCR detection system (BioRad, Hercules, California, US). Cycling conditions included a pre-incubation step (95°C for 3 minutes) and 40 cycles of amplification (95°C for 30 seconds, 60°C for 1 minute). For gene expression, total RNA (2 μg) was reverse transcribed using the SuperScript II Reverse Transcription kit, 50 ng of random primers (both from Invitrogen, Life technologies, Carlsbad, California, USA) and 0.5 mM dNTP mixture (Bioline, London, UK). QPCR was done in a CFX96 real-time PCR detection system using 30 ng of cDNA, 1X iQ SYBR-green supermix (Biorad, Hercules, California, US) and primers at a final concentration of 2 mM. Cycling conditions included a pre-incubation step (95°C for 3 minutes), 40 cycles of amplification (95°C for 30 seconds, annealing temperature for 30 seconds and 72°C for 30 seconds) and a melting curve ( 95°C for 1 minute, 55°C for 1 minute, denaturation from 55°C to 95°C at 0.5°C/10 seconds increments). Each reaction was done in triplicate. All primers used are listed in Additional file 1: Table S8.

Expression values were normalised to the geometrical mean of housekeeping genes RPL11, RPL30 and RPS6 and fold-changes between treatments and controls were determined by the 2^-ΔΔCt method [57] as implemented in the SLqPCR Bioconductor package (http://www.bioconductor.org/packages/2.12/bioc/html/SLqPCR.html). The variance between sample groups was assessed through the Barlett test for variance. If variances could be assumed as equal between groups, significant differences were established using ANOVA followed by pairwise t-test. If variances could not be assumed as equal, significant differences were then assessed using a one way test for equal means followed by pairwise t-test. After multiple test correction by Benjamini and Hochberg method [54], results were considered significant at two levels: with adjusted p value lower or equal to 0.05 and with adjusted p value lower or equal to 0.08. All analysis was done using R v9 (http://www.r-project.org).

Protein immunoblotting

Cells were lysed in urea/SDS buffer (6.7 M urea, 10 mM Tris-Cl pH 7.4, 10% glycerol, 1% SDS) supplemented with 1 mM dithiothreitol. Extracts were resolved by SDS-PAGE, electroblotted onto PVDF membranes (Millipore) and probed with primary antibodies followed by HRP-conjugated secondary antibodies. SuperSignal Chemiluminescent Substrates (Pierce) were used to visualize immunoreactive species. Antibodies used: anti-HIF1α 610959 from BD Biosciences, anti-HIF2α 190b from Peter Ratcliffe Lab, and anti-β-Actin/HPR ab49900 from Abcam.

Short interfering RNA treatment of MCF-7 cells

RNA oligonucleotides (Ambion, Life technologies, Carlsbad, California, USA) were used for siRNA suppression of individual HIF family members (HIF-1α: sense, CUGAUGACCAGCAACUUGAtt; antisense, UCAAGUUGCUGGUCAUCAGtt. HIF-2α; sense, CAGCAUCUUUGAUAGCAGUtt; antisense, ACUGCUAUCAAAGAUGCUGtt). A set of oligonucleotides (Ambion, Life technologies, Carlsbad, California, USA) against a sequence from the Drosophila melanogaster HIF sequence, which lacks any substantial sequence similarity with human HIF or other genes, was used as a negative control (D-HIF: sense, CCUACAUCCCGAUCGAUGAtt; antisense, UCAUCGAUCGGGAUGUAGGtt) in one of the experiments. MCF-7 cells were seeded at 30% confluence and grown in normoxic conditions. Cells were transfected twice with 20 nmol/L of siRNA at 24 and 48 hours using OligofectAMINE (Invitrogen, Life technologies, Carlsbad, California, USA) according to the instructions of the manufacturer. At 55 hours, cells were exposed to 1% oxygen, and after an additional 16 hours, RNA was extracted. Each siRNA treatment was done in triplicate as well as the parallel control using OligofectAMINE reagent alone. The same protocol was used for the generation of protein extracts. We ensured that specific and substantial knockdown was occurring in MCF-7 cells under hypoxia by assaying HIF-α levels by immunoblotting (as previously described [26]).

Statistical methods

A signature of 99 genes up-regulated under hypoxic conditions has been derived in head and neck tumors in vivo and it has been assessed for its prognostic value in a head and neck and breast cancer data set [31]. Here, this signature was used to calculate a hypoxic score for a breast cancer series using the gene expression profiles (GSE22220) published in [32]. The score is a measure of the level of expression of genes in the hypoxia signature: a high positive score reflects a higher than average expression of the genes in the hypoxia signature and indicates that the sample is hypoxic, a low negative score reflects a lower than average expression of the genes in the hypoxia signature and indicates that the sample is normoxic [31]. Correlation of microRNA expression from the same breast cancer series (GSE22220) and hypoxia score was assessed using Spearman’s rank tests. Statistical analyses were done using R.