MDA-MB468 and DLD-1 cells were cultured in DMEM/HamsF12 with L-glutamine (PAA Laboratories; GE Healthcare Bio-Sciences Austria GmbH); LNCaP cells were cultured in RPMI-1640 medium with L-glutamine (Invitrogen; Thermo Fisher Scientific, Inc.). All cells were supplemented with 10% Fetal Bovine Serum Gold (PAA Laboratories; GE Healthcare Bio-Sciences Austria GmbH). Hera cell 240 (MultiTemp Scientific AG), IG150 (Jouan; Thermo Fisher Scientific, Inc.) and CB 53 (Binder GmbH) incubators were used for cell culture. For normoxic conditions, cells were cultured in room air with 5% CO₂, whereas hypoxic incubation was conducted with 1% O₂ and 5% CO₂; both at 37°C in H₂O-saturated atmosphere. MDA-MB468, DLD-1 and LNCaP cells were authenticated by SNP typing (Multiplexion GmbH) and approved for mycoplasma negativity using the VenorGeM Mycoplasma PCR detection kit (Minerva Biolabs GmbH) prior to experiments. Transient MB-knockdown cells for all three cell lines were generated using small interfering (si)RNA. Cells at 50-60% confluency were seeded in 24-well plates. After 24 h of incubation, 20 pmol MB siRNA (HS_MB_6, 2607575; anti-sense sequence, 5′-UGA UUA AUC AGA CAA UUG CTA-3′; Qiagen GmbH) or scrambled (scr) RNA (cat. no. 1022076; Qiagen GmbH) were diluted in Opti-MEM (cat. no. 31985062; Thermo Fisher Scientific, Inc.) for each well. Lipofectamine^® 2000 (1 µl; Thermo Fisher Scientific, Inc.) was added to 50 µl Opti-MEM and incubated for 5 min at room temperature. The two prepared solutions were then gently mixed together and incubated for 20 min at room temperature. The oligomer-Lipofectamine 2000 complexes were then added to the cells. After 24 h of incubation, the culture medium was changed.

Following 72 h of incubation under 1% O₂ or normoxia, RNA was extracted from cells using the RNeasy Mini kit (Qiagen GmbH), including the RNase-Free DNase I treatment (Qiagen GmbH) for DNA digestion. RNA integrity was confirmed using the Agilent RNA 6000 Nano kit on a 2100 Bioanalyzer (Agilent Technologies GmbH). RNA samples from four replicate cell passages were pooled for each condition and cell line. Illumina TruSeq Total RNA transcriptome libraries (Illumina, Inc.) were generated for the pools of the siMB cell types (3 cell lines and 2 O₂ conditions) and the respective MB⁺ controls (StarSEQ GmbH). The 50 bp single-end Illumina sequencing was performed by the next generation sequencing core facility of the Biology Department of Johannes Gutenberg University (Mainz, Germany) on an Illumina HiSeq2500 sequencing platform.

Using the CLC Genomics Workbench 8.0.1 sequence trimmer (Qiagen GmbH), 12 terminal nucleotides from the 5′ end, remaining Illumina adapters and low-quality sequences (below Phred 13) were removed from each read. The minimum sequence length was set to 15 bp, allowing ≤2 ambiguous nucleotides per read. The reads were then mapped to the annotated human genome hg19 (ftp.ensembl.org/pub/release-75/gtf/homo_sapiens/Homo_sapiens.GRCh37.75. gtf.gz) with the CLC Genomics Workbench 8.0.1 RNA-Seq tool. Mapping parameters were set as default, but included intergenic regions and allowed ≤10 hits per read. Paired two-group comparisons were conducted between the mapping results of siMB and MB⁺ cells for the three cell lines with or without hypoxia using the CLC Genomics Workbench 8.0.1 with the original reads per kilobase of transcript per million mapped reads (RPKM) values (30). Proportion-based statistics were calculated with Kal's Z-tests [Bonferroni- and false discovery rate (FDR)-corrected] to identify statistically significant differentially expressed genes at [−log₁₀(p)]≥0.5 (31).

Low frequency variant detection was performed with the CLC Genomics Workbench 8.5.1 at default parameters using a minimum read count of 6, a minimum coverage of six reads and including a read direction filter of 1%. Only variants of QUAL ≥100 (probability that a particular variant exists ≥0.9999999999) were listed and annotated with human genome hg19 gene track overlap information to infer potential non-synonymous amino acid changes.

To access the transcriptome-wide effects of MB expression in epithelial cancer cells, all genes that were differentially expressed between siMB and MB⁺ cells were identified (Table I). A list of genes that were differentially expressed depending on the MB status in at least four of the six gene lists was compiled to filter general MB functions common to the majority of the investigated cell types. This merged gene list was used for gene set enrichment analysis (hypergeometric tests; FDR-corrected P<0.05) using the Cytoscape 3 application BiNGO 3.0.3 (32) to generate a map of significantly overrepresented Gene Ontology (GO)-categories.

For a more detailed analysis of the impact of MB in each cell line under each O₂ condition, the six lists of genes differentially expressed between siMB and MB⁺ cells in matching conditions were further analyzed for GO term enrichment (P<0.05) and gene enrichment in KEGG and Biocarta pathways using the DAVID Bioinformatics Resources 6.7 functional annotation tool (33,34), choosing either the upregulated or downregulated genes as inputs. Genes differentially expressed between siMB and MB⁺ cells were further annotated and interpreted by Core Ingenuity Pathway analysis (IPA) (http://www.ingenuity.com). Analyses were conducted using the fold-changes determined by CLC two-group-comparisons, including the IPA a priori knowledge of direct and indirect relationships between genes observed in all human tissues. For visualization, a list of significantly active upstream regulators in each condition was compiled based on the direction of regulation of their target genes.

To investigate the function of endogenously expressed MB in epithelial cancer cells, siRNA was used to generate MB-knockdown breast (MDA-MB468), prostate (LNCaP) and colon (DLD-1) cancer cells, and their transcriptome profiles were compared to those of matching MB-wild-type cells treated with scr siRNA as a control. Different tumor types were selected for studying MB expression to discriminate tumor-specific effects [e.g., ER status (27)] from common changes that may be associated with MB expression throughout different tumor types of epithelial origin. As MB can be either oxygenated or deoxygenated, experiments for all three cell lines were conducted in room air (normoxia) and 1% O₂ (hypoxia), the latter causing a fractional MB O₂ desaturation of ~42% (35). To specifically study the impact of MB in cells adapted to long-term hypoxia, mimicking tumors, the cells were incubated for 72 h at hypoxic vs. normoxic conditions; previous experiments on MDA-MB468 siRNA MB-knockdown cells demonstrated strong phenotypic effects at 72 h (28). Using Illumina transcriptome sequencing and read mapping to the annotated human genome, gene expression profiles were generated for each cell line and condition. The numbers of successfully mapped reads ranged between 23.7 and 53.3 Mio (Table I).

Following RNAi, MB was downregulated 6- and 8-fold in DLD-1, 11- and 16-fold in LNCaP and 26- and 38-fold in MDA-MB468 cells compared with the control cells in normoxic and hypoxic conditions, respectively. In addition, between 160 and 857 genes were differentially expressed between siMB and MB⁺ control cells under equivalent oxygen conditions (Tables I and SI). This divergence implied differences in the impact of MB on specific epithelial cell lines and encouraged the investigation of the generic and cell-type specific effects of MB expression in epithelial cancer cells. Differentially expressed genes were subjected to GO term analysis, and overrepresentations of GO terms within gene lists were statistically evaluated.

The effects of MB expression common among the three epithelial cancer cell lines were summarized by BiNGO analysis of a subset of 145 genes, which were consistently differentially expressed in at least four of the six experimental groups (Table I). A number of overrepresented GO terms, presented as colored nodes in Fig. 1, were associated with metabolic shifts (i.e., 'glycolysis', 'fatty acid biosynthetic process', 'oxidative phosphorylation' and 'generation of precursor metabolites and energy') and anti-apoptotic cell signaling (i.e., 'anti-apoptosis', 'negative regulation of programmed cell death' and 'regulation of DNA damage response', signaling transduction by p53 class mediator'). These terms were linked to the overrepresented nodes 'cell cycle' and 'protein ubiquiti-nation' in the BiNGO graph. In addition, overrepresentation of the GO categories 'response to hypoxia', 'response to oxidative stress' and 'response to ROS' was observed, suggesting that MB also participates in cellular response to oxidative stress.

The present study further aimed to identify individual cell type- and O₂ level-specific characteristics in the three cancer cell lines (Table II). GO term and KEGG and Biocarta pathway enrichment analyses were performed separately for each subset of significantly upregulated and downregulated genes to compare siMB and MB⁺ cells (Table SI). Lists of enriched GO terms revealed a number of categories that repeatedly appeared across all studied cell types. These functional terms were related to those identified before in the global BiNGO analysis. In addition, the terms were linked to a range of known molecular functions exerted by MB. Several GO categories unique either to a certain cell type or to the O₂ level were also identified (Table II).

In prostate and breast cancer cells, the GO categories 'response to O₂ levels' and 'response to hypoxia' were overrepresented in the set of downregulated genes in siMB cells independent of the O₂ conditions (Table II). Since these genes included a number of known hypoxia marker genes (Table SI), these results suggested that MB may modulate the cellular hypoxia response in breast and prostate cancer cell lines.

GO enrichment analysis further indicated a metabolic shift towards a reduced rate of glycolysis in MB-knockdown prostate and breast cancer cells under the two O₂ conditions. This was evidenced by the downregulation of the hypoxia-regulated lactate dehydrogenase A (LDHA) gene and the glycolytic genes triosephospate isomerase 1 (TPI1), GAPDH, protein kinase CGMP-dependent 1 (PGK1) and glycerol-3-phosphate acyltransferase, mitochondrial (GPAM1). By contrast, MB-knockdown in hypoxic colon cancer cells resulted in an increase in glycolytic flux, as LDHA and the glycolysis genes TPI1, PGK1, GPAM1 and enolase 1 were upregulated (Tables SI and SII).

In normoxic MB-knockdown breast cancer cells, the GO categories 'unsaturated fatty acid metabolic process' and 'unsaturated fatty acid biosynthetic process' were overrepresented in the set of downregulated genes (Table II), including fatty acid desaturase 3 and δ4-desaturase, sphingolipid 1 (Tables SI and SIII). By contrast, genes associated with sterol, cholesterol and steroid metabolic processes were upregulated in MB-knockdown breast cancer cells, such as sterol regulatory element binding transcription factor 1 (SREBF1), phosphomevalonate kinase, 24-dehydrocholesterol reductase (DHCR24) and cytochrome P450 1B1 (Tables SI and SIII). These results suggested that oxy-MB may stimulate the metabolism of unsaturated fatty acids in oxygenated breast cancer cells. However, the globin may also interfere with the synthesis of sterols.

In prostate cancer cells, siMB treatment resulted in an increased ability to metabolize or synthetize fatty acids, sterols and eicosanoids, indicated by the upregulation of hydroxyacyl-CoA dehydrogenase trifunctional multienzyme complex β, DHCR24, sterol-C5-desaturase, 15-hydroxypros-taglandin dehydrogenase (HPGD), elongation of very long chain fatty acids protein 5 (ELOVL5) and acyl-CoA synthetase long chain family member 3 in normoxic siMB cells, and by upregulation of ELOVL5, HPGD, fatty acid synthase and phospholipase A2 group V in siMB hypoxic cells (Tables SI and SIII). By contrast, in hypoxic MB-knockdown colon cancer cells, the expression of genes associated with cholesterol and sterol metabolism and biosynthesis [DHCR7, 3-hydroxy-3-methylglutaryl-CoA synthase 1 (HMGCS1), farnesyl diphosphate synthase, farnesyl-diphosphate farne-syltransferase 1 and SREBF1] was decreased (Tables SI and SIII). In summary, the results of the transcriptome analysis suggested that MB may regulate fatty acid metabolism in a cell-type specific manner.

The GO categories 'induction of programmed cell death' and 'induction of apoptosis' were overrepresented in the set of downregulated genes in hypoxic MB-knockdown prostate cancer cells, suggesting that MB positivity may be associated with the stimulation of apoptotic processes in these cells. One option of maintaining apoptosis is by triggering cell cycle checkpoint signaling (36). Accordingly, in both normoxic and hypoxic siMB LNCaP cells, the GO categories 'regulation of cell cycle' and 'positive regulation of cell proliferation' were overrepresented among the downregulated genes, whereas the category 'negative regulation of cell proliferation' was enriched in the upregulated genes. In addition, the G₁/S cell cycle checkpoint regulator cyclin-dependent kinase inhibitor 1a (CDKN1A, also termed TP21) was downregulated in siMB LNCaP cells. These results indicated direct or indirect involvement of MB in TP21 expression and a possible induction of cell cycle arrest in the presence of MB. However, this scenario is most probably of a more complex nature, since a subset of cell cycle progression-associated genes such as cyclin B1 (CCNB1) and CCND1 were also downregulated in siMB LNCaP cells (Tables SI and SIII). Enrichment analyses also revealed that the GO categories 'cell motion' and 'cell migration' were significantly overrepresented among the upregulated genes of siMB LNCaP cells (Table II), including activated leukocyte cell adhesion molecule, CD9, RAC1, laminin subunit γ1 and thrombospondin (THBS1) under both O₂ conditions (Tables SI and SIII). These observations strengthened our hypothesis that MB induced apoptosis and decreased the migratory and motion capacity of prostate cancer cells.

In breast cancer cells, genes associated with apoptosis were overrepresented in the sets of up- and downregulated MDA-MB468 cells under both O₂ conditions (Table II), which made it difficult to discriminate the impact of MB on promoting or inhibiting apoptosis. The GO terms 'mitotic cell cycle', 'cell cycle process' and 'positive regulation of cell proliferation' were overrepresented in the downregulated genes of siMB cells, whereas only a limited number of negative regulators of cell proliferation were upregulated (including retinoic acid receptor responder 3, prostaglandin E synthase and CDKN2C; Tables SI and SIII). These results may imply increased proliferation rates in MB-positive breast cancer cells. In addition, the GO terms 'cell motion' and 'cell migration' were significantly overrepresented in the downregulated genes of siMB cells (Table II), as evidenced by the lower expression of THBS1, TNF receptor superfamily member 12A, vinculin, moesin and tensin 3 (Tables SI and SIII). Thus, the predictions made about the impact of MB on cell motion in hypoxic and normoxic breast cancer cells contradicted the results observed in prostate cancer cells.

In colon cancer cells, genes associated with apoptosis were also overrepresented in both up- and downregulated genes (Table II). The term 'cell cycle process' was overrepresented in the upregulated genes of hypoxic and normoxic siMB DLD-1 cells, including TPX2 microtubule nucleation factor, CDC20, non-SMC condensin I complex subunit D2, family with sequence similarity 83 member D, CCNB1 and CCNB2. In addition, the term 'cell cycle arrest' was overrepresented in the downregulated genes of hypoxic siMB samples (Table II) due to the downregulation of genes such as CDKN1C, CDKN1A, protein phosphatase 1 regulatory subunit 15A and DNA damage inducible transcript 3 (Tables SI and SIII). Therefore, the DLD-1 data suggested that MB expression may mediate the suppression of cell cycle progression in colon cancer cells, possibly via apoptotic pathways.

In line with the results of the prostate cancer cell line, the GO term 'cell motion' was also overrepresented in the down-regulated genes of normoxic siMB DLD-1 cells.

In hypoxic and normoxic MB-knockdown breast cancer cells, the GO term 'regulation of NO^· biosynthetic process' and other GO categories referring to the biosynthetic process of N-compounds were overrepresented in the downregulated genes, suggesting an increased production of NO^· in hypoxic and normoxic MB-positive breast cancer cells (Table SIII). This dependency between MB expression and NO^·-related GO terms was also observed in prostate and colon cancer cells, which suggested involvement of MB in NO^· metabolism in all three cell types. All siMB cells also exhibited enrichment of the GO category 'response to ROS' in the downregulated genes irrespective of the O₂ culture conditions. Therefore, MB may be involved in increasing ROS levels in all tested cell lines.

To predict the differential activation of key transcriptional regulators in siMB vs. MB⁺ cells, the mRNA expression of their target genes was examined using bioinformatic IPA upstream regulator analysis (Table III). A subset of 40 transcription factors was suggested to be differentially active in at least one cell line at hypoxia or normoxia, with high Z-scores indicating strong activation, and low Z-scores indicating transcription factor inhibition. All siMB cells subjected to hypoxia and normoxic MDA-MB468 cells were predicted to exhibit decreased HIF1α activity compared with the MB⁺ controls. This was supported by the downregulation of HIF1α target genes in the siMB cell groups, including TP21 and nucleophosmin 1 in all cell lines, LDHA in LNCaP and MDA-MB468 cells, aldolase fructose-bisphosphate A and interleukin 8 in MDA-MB468 and DLD-1 cells, vascular endothelial growth factor A in LNCaP cells and HIF1A, solute carrier family 2 member 1 and carbonic anhydrase 9 in MDA-MB468 cells (Table SI). Thus, MB may activate HIF1α in breast cancer cells and, to a lesser extent, in prostate and colon cancer cells.

IPA upstream activator analysis also suggested the involvement of p53 signaling in the three cancer cell lines (Table III). Of note, in our set up, only LNCaP cells encoded the TP53 wild-type-like Pro72Arg variant, whereas the MDA-MB468 and DLD-1 cells carried TP53 mutations Arg273His and Ser241Phe, respectively (Table SII), which affected the key residues of DNA-binding (37). RNA-Seq analysis of LNCaP cells demonstrated that MB⁺ was associated with the upregulation of p53 target genes including TP21, BCL2-interacting killer and heme oxygenase 1 (HMOX1) (Table SI), which suggested an activation of p53 signaling by MB in the prostate cancer model. In MDA-MB468 and DLD-1 cells expressing MB, even the mutant p53 was active, which was evidenced by the decreased expression of p63 target genes (Table III), possibly due to lowered p63 activity by the binding of mutated p53 (38,39). Accordingly, IPA analysis results suggested activation of p63 in MB⁺ hypoxic LNCaP cells expressing the wild-type-like TP53 mutation Pro72Arg (Table III). Another fitting observation noted for the breast cancer cell line was enhanced transcriptional activity predicted for the oncogenic factor RELA (Table III), which is known to be activated by mutant p53 (40).

In line with the observed metabolic shifts in MB⁺ cancer cells, the gluconeogenesis and adipogenesis regulator fork-head box O1 was another transcription factor predicted to be dysregulated in the three cell models (Table II).