A Systems Biology Approach to Understanding the Pathophysiology of High-Grade Serous Ovarian Cancer: Focus on Iron and Fatty Acid Metabolism.

Ovarian cancer (OVC) is the most lethal of the gynecological malignancies, with diagnosis often occurring during advanced stages of the disease. Moreover, a majority of cases become refractory to chemotherapeutic approaches. Therefore, it is important to improve our understanding of the molecular dependencies underlying the disease to identify novel diagnostic and precision therapeutics for OVC. Cancer cells are known to sequester iron, which can potentiate cancer progression through mechanisms that have not yet been completely elucidated. We developed an algorithm to identify novel links between iron and pathways implicated in high-grade serous ovarian cancer (HGSOC), the most common and deadliest subtype of OVC, using microarray gene expression data from both clinical sources and an experimental model. Using our approach, we identified several links between fatty acid (FA) and iron metabolism, and subsequently developed a network for iron involvement in FA metabolism in HGSOC. FA import and synthesis pathways are upregulated in HGSOC and other cancers, but a link between these processes and iron-related genes has not yet been identified. We used the network to derive hypotheses of specific mechanisms by which iron and iron-related genes impact and interact with FA metabolic pathways to promote tumorigenesis. These results suggest a novel mechanism by which iron sequestration by cancer cells can potentiate cancer progression, and may provide novel targets for use in diagnosis and/or treatment of HGSOC.

Introduction

Ovarian cancer (OVC) is the sixth most common cancer worldwide, and will account for ∼14,000 deaths in the United States in 2017 (Siegel et al., 2017). OVC is the most lethal (i.e., has the lowest 5-year survival rate, ∼45%) of all gynecological cancers. High-grade serous ovarian cancer (HGSOC) is the most common and deadliest subtype, with a 5-year survival rate of ∼35% (Sankaranarayanan and Ferlay, 2006). The distal fallopian tube is believed to be the site of origin for HGSOC (Vang et al., 2009).

High levels of systemic iron are associated with increased risk for a number of cancers. With respect to dietary iron, a meta-analysis of 59 epidemiological studies showed increased relative risk for colorectal, colon, lung, and to a smaller extent breast cancers with increases in heme iron intake (Fonseca-Nunes et al., 2014). Elevated systemic levels of iron can occur independent of diet, such as with the genetic disorder hereditary hemochromatosis (HH). Patients with HH are at especially high risk for liver cancer (Elmberg et al., 2003; Hsing et al., 1995), as well as other nonhepatic cancers (Hsing et al., 1995; Osborne et al., 2010; Shaheen et al., 2003).

Cancer cells themselves often exhibit a phenotype of iron retention, with increased levels of proteins involved in iron import and decreased levels of proteins involved in iron export when compared to normal cells (Torti and Torti, 2013). HGSOC conforms to this paradigm (Basuli et al., 2017). Moreover, iron derived from hemolysis of erythrocytes, which are present due to menstruation, may directly contribute to the initiation and progression of OVC (this is referred to as the incessant menstruation hypothesis) (Vercellini et al., 2011). In addition to contributing to cellular oxidative stress, iron is required for deoxyribonucleic acid replication and is implicated in several oncogenic processes and signaling cascades (Xue and Shah, 2013).

However, the precise role of iron in OVC remains unclear. We were thus motivated to consider how iron can impact the progression of OVC from a systems viewpoint. To do so, we analyzed clinical and cell culture gene expression data from multiple sources and searched for perturbed pathways that showed involvement of iron-related genes in a majority of the datasets. We considered that a consistency in results among different datasets indicated a robust mechanism.

Using our analysis, we identified an involvement of iron-related genes in pathways involved in fatty acid (FA) import and synthesis. It has been established that cancer cells increase FA synthesis and import to aid in generation of phospholipids for cell division and lipid-mediated signaling (Currie et al., 2013). We propose a novel link between increased FA synthesis and import and the cancer-associated increase in intracellular iron in HGSOC.

While a link between iron and FA synthesis has not been investigated in the context of cancer, the effects of iron on FA metabolism in the normal liver have been studied, although frequently with contradictory results (Ahmed et al., 2012). For example, studies have shown that iron deficiency can either increase [e.g., (Sherman, 1978)] or decrease [e.g., (Stangl and Kirchgessner, 1998)] hepatic lipogenesis, with the diverse methodology used to collect data cited as a potential cause of the discrepancy (Ahmed et al., 2012).

Nevertheless, it is agreed that iron can have a major impact on FA metabolism through several mechanisms: one is by directly oxidizing FAs by the Fenton reaction (Bacon and Britton, 1990), the second is by contributing to the activity of heme- and iron-binding enzymes involved in FA metabolism (Stangl and Kirchgessner, 1998), and the third is indirectly by the downstream products of iron deficiency or overload (Davis et al., 2012). Such actions of iron are not liver specific and can transfer to the setting of high iron load in HGSOC and other cancers, as we show.

To establish a systems-level understanding of the relationship between increased iron utilization and increased FA import and synthesis, we use the results of our analysis in tandem with the literature to create a network of cancer-associated FA import and synthesis pathways that are dependent on the activity and expression of iron-related genes. From this network, we derive hypotheses regarding how the action of iron and iron-related genes can promote tumor growth through their activity in subcomponents of the FA metabolic network.

Materials and Methods

Data sources

Dataset [1] comes from The Cancer Genome Atlas Research Network (TCGA) (Cancer Genome Atlas Research Network, 2011) and dataset [2] from the study by Tothill et al. (2008), both which include mRNA expression profiling data on HGSOC (Stages II–IV, Grade 2–3) biopsy samples. Since all cell types are included in the biopsy samples, the data cannot give a clear understanding of which cell types express which genes. However, they can provide an overview of gene expression in the tumor microenvironment and may identify important paracrine interactions between the distinct resident cell types.

The next dataset [3] is from Tone et al. (2008) where fallopian epithelial cells were extracted from subjects who did (FTPb) and did not (FTPn) have known BRCA mutations. Epithelial cells were also extracted from patients with high-grade serous fallopian and ovarian cancer, and it was shown that gene expression between the two cancer types was indistinguishable, and thus we group the two cancer datasets into the HGSOC set. The Tone et al. data represent epithelial (stem and differentiated) cells from normal and cancerous tissue.

Finally, we use expression data (dataset [4]) from primary follicular stem cells, which are thought to be the cells of origin of HGSOC, and model cancer stem cells derived from them (Yamamoto et al., 2016). Yamamoto et al. (2016) obtained the cells through primary culture of fallopian tube stem cells (FTstem) followed by immortalization and transformation with h-TERT, SV40 large T antigen, and c-Myc (FTtr), and further characterized them by gene expression profiling at the time of isolation. The FT cells represent an in vitro model system of normal and cancer ovarian stem cells, and can thus be considered a model subset of the Tone et al. data (Fig. 1). The data sources are summarized in Table 1.

Overview of the datasets used with respect to the tissues (normal ovary/fallopian tube and HGSOC) and cell types (stem, nonstem, and microenvironment) that they represent. Numbers in parentheses correspond to the numbering of the datasets in the text. [1] TCGA (Cancer Genome Atlas Research Network, 2011) and [2] Tothill et al. (2008) datasets come from all cells derived from HGSOC tissue biopsies, [3] Tone et al. (2008) includes samples from nonmalignant epithelial fallopian and HGSOC cells, and [4] are model fallopian stem cells, the predicted precursors of HGSOC, which have been transformed to cancer stem cells (Yamamoto et al., 2016). TCGA, The Cancer Genome Atlas Research Network.

1.

Datasets Used for Analysis, Total Number of Samples (

Dataset	[1] TCGA	[2] Tothill et al.	[3] Tone et al.	[4] FT model cell
n	489	218	37	6
Cell type	Biopsy, heterogeneous	Biopsy, heterogeneous	Biopsy, epithelial cells	Model cell line
Platform	hgu133a	hgu133plus2	hgu133plus2	hgu133plus2
Comparison	Mes. (n = 108) v. Diff. (n = 134) subtypes	CLOVAR Mes. (n = 80) v. Diff. (n = 35) subtypes	HGSOC samples (n = 13) v. normal (n = 12) fallopian	HGSOC cancer stem (FTtr) (n = 2) v. normal fallopian stem (Ftstem) (n = 2) model cell line
Public access	TCGA Data Portal (2014)	GEO acc: GSE9899	GEO acc: GSE10971	GEO acc: GSE69428
Reference	Cancer Genome Atlas Research Network (2011)	Tothill et al. (2008)	Tone et al. (2008)	Yamamoto et al. (2016)

TCGA, The Cancer Genome Atlas Research Network.

Data processing

Affymetrix HGU-133A CEL files (1) from TCGA (Cancer Genome Atlas Research Network, 2011) and HGU-U133plus2 CEL files from Tothill et al. (2008) (GSE9891), Tone et al. (2008) (GSE10971), and FT cell lines (Yamamoto et al., 2016) were robust multiarray average preprocessed and filtered for markers (without the use of variance-based methods) in R (v.3.4.0) using the rma and nsFilter functions, respectively, from the simpleaffy package (v.2.52.0) (Wilson and Miller, 2005). The probe set with the highest average expression per gene in each dataset was selected to represent the expression of that gene.

The comparisons, as described below, for each dataset were between the following subsets (Table 2): for [1] TCGA, Mesenchymal (Mes) (n = 108) were compared against Differentiated (Diff) (n = 134) samples, for [2] Tothill et al. (2008), CLOVAR Mes (n = 80) samples were compared against CLOVAR Diff (n = 35) samples, for [3] Tone et al. (2008), normal fallopian tube epithelium (FTEn) samples (n = 12) were compared against high-grade serous ovarian cancer (HGSOC) samples (n = 13), and for [4] FT cells, c-Myc-transformed FTtr samples (n = 2) were compared against normal FTstem samples (n = 2). Subtype classification for datasets [1] and [2] was obtained from Supplementary Table S1 in (Verhaak et al., 2013).

2.

The Top Three Mutually Perturbed Pathways for the Three Datasets Indicated

	KEGG Pathway	FDR p value	% iron DE genes	Iron DE genes
[1] TCGA Cancer Genome Atlas Research Network (2011)	Mineral absorption	9.33E-22	44.44	FTL, HMOX1, SLC11A2, HEPH
	PPAR signaling pathway	3.10E-19	9.09	SCD, FADS2
	Metabolism of xenobiotics by cytochrome P450	1.46E-05	11.76	CYP1A2, CYP1B1
[2] Tothill et al. (2008)	PPAR signaling pathway	4.72E-32	13.04	CYP27A1, SCD, FADS2
	Metabolism of xenobiotics by cytochrome P450	1.71E-22	23.81	CYP1A2, CYP1B1, CYP2A13, CYP3A5, CYP2S1
	Mineral absorption	2.50E-22	40.00	FTL, HMOX1, CYBRD1, HEPH
[3] Tone et al. (2008)	Mineral absorption	3.42E-06	31.58	SLC46A1, FTH1, FTL, STEAP2, SLC40A1, HMOX1
	PPAR signaling pathway	3.37E-03	20.00	CYP27A1, SCD, SCD5, FADS2
	Metabolism of xenobiotics by cytochrome P450	8.08E-03	16.67	CYP1B1, CYP2E1

“Iron DE genes” refers to iron-related genes that are differentially expressed and “% iron DE genes” refers to the percent of iron DE genes with respect to all DE genes in that pathway.

KEGG, Kyoto Encyclopedia of Genes and Genomes; PPAR, peroxisome proliferator-activated receptors.

Genes that were differentially expressed within each dataset were found using the linear modeling paradigm limma in R/Bioconductor (Ritchie et al., 2015; Smyth, 2004). Significantly perturbed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways within each dataset were found using the gene set analysis method Generally Applicable Gene-set Enrichment (GAGE) (Luo et al., 2009), using the bidirectional option (same.dir = F) and unpaired sample setting (compare = “unpaired”).

Comparison of subtypes within the datasets

The decision to compare the Mes and Diff subtypes in dataset [1] and CLOVAR Mes and Diff subtypes in dataset [2] was made based on the following rationale: while none of the original molecular subtypes showed significant differences in survival in TCGA (Cancer Genome Atlas Research Network, 2011), Verhaak et al. (2013) showed that an optimized classification algorithm using the same training data and original four profiles as in the original study gave four modified profiles (termed CLOVAR Diff, Pro, Imr, and Mes) that were different in terms of survival. In particular, when applied to test data of 800 samples, the CLOVAR Mes subtype had significantly lower survival than the CLOVAR Diff subtype.

The classification by TCGA was based on the expression values of 1,500 of the most variable genes from an analysis of multiple micorarrays from TCGA, whereas the CLOVAR classification was based on the 100 genes most correlated or anticorrelated with survival of the TCGA patients (Cancer Genome Atlas Research Network, 2011; Verhaak et al., 2013). We used both the [1] Mes v. Diff and [2] CLOVAR Mes v. Diff subsets to search for impact of iron-related genes on pathways perturbed in the Mes v. Diff samples since the original classification could potentially more accurately identify molecular subtypes, whereas the CLOVAR-based classification was part of the test set that contributed to significant differences in survival between CLOVAR subtypes.

Iron-related genes and mutually perturbed pathway analysis

Genes that bind iron or are influential in iron-related homeostasis were considered part of a set of iron-related genes. This list of 274 iron-related genes was compiled manually from the literature and the gene ontology (GO) category of iron ion binding (GO:0005506) (Ashburner et al., 2000; Thomas, 2017). The list includes genes involved in iron metabolism itself (i.e., iron uptake, storage, and efflux), genes encoding proteins that use iron, heme, or iron-sulfur clusters as cofactors, as well as genes that regulate proteins of iron metabolism, and is provided as Supplementary Table S1.

For each KEGG pathway identified by GAGE to be significantly perturbed (FDR adjusted p value <0.05) per dataset, the number of iron-related genes that were differentially expressed (as indicated by limma) and had log2-fold change ≥|0.2| was counted, and the percentage of genes that were iron-related genes in the full list of differentially expressed genes per pathway was calculated. A pathway was considered mutually perturbed if it was found to be perturbed in datasets [1]-[3], and had ≥2 differentially expressed iron-related genes. A summary of all mutually perturbed pathways identified is found in Supplementary Table S2, and an ordering of the pathways in ascending order of cumulative rank is found in Supplementary Table S3. The full results of the GAGE pathway analysis for each dataset, along with differentially expressed iron-related genes for each pathway, are found in Supplementary Tables S4–S7.

Overrepresentation analysis (ORA) of iron-related genes in KEGG pathways by the hypergeometric test was performed using the ConsensusPathDB database (Kamburov et al., 2011).

Results

Overview of data used for analysis

We utilize a variety of datasets to investigate the impact of iron on HGSOC (Table 1, Fig. 1). The datasets consist of microarray expression data from biopsy samples of HGSOC patients (Stages II–IV, Grade 2–3) collected by [1] TCGA (Cancer Genome Atlas Research Network, 2011) and [2] Tothill et al. (2008), [3] laser capture microdissected fallopian epithelial cells from healthy controls and patients with high-grade serous fallopian and ovarian cancer (Tone et al., 2008), and [4] primary follicular stem cells FTstem and model OVC stem cells (FTtr, “tr” stands for transformed) derived from them (Yamamoto et al., 2016).

Dataset [2] does not contain samples from nonmalignant tissue, and dataset [1] contains eight such samples when queried in June, 2017 (Broad Institute TCGA Genome Data Analysis Center, 2016), which we deemed were not an appropriate sample size to compare against the qualitatively larger number of tumor samples in the database. Thus, in these datasets, we compared molecular subtypes that differ in survival. Specifically, HGSOC samples have been classified into four molecular profiles: Differentiated (Diff), Proliferative (Pro), Immunoreactive (Imr), and Mesenchymal (Mes) (Cancer Genome Atlas Research Network, 2011).

A simplified classification algorithm was developed by (Verhaak et al., 2013) and termed CLOVAR. Using dataset [2] as part of a test set, Verhaak et al. (2013) showed that the CLOVAR Mes subtype showed significantly lower survival than the CLOVAR Diff subtype. We therefore chose to compare the Mes v. Diff subtypes as a means of comparing subsets with higher (Diff) v. lower (Mes) survival. Datasets [3] and [4] do contain normal, or model normal, samples, and thus a comparison between HGSOC and normal samples was performed for these two datasets (see Table 1 for a summary of comparisons for all datasets). See the Materials and Methods section for more details on the datasets and comparisons.

Mutually perturbed pathways with differentially expressed iron-related genes

We define iron-related genes as genes with a GO annotation of “iron-ion binding,” as well as genes known to function in intracellular iron regulation (transferrin receptor [TFRC], ferroportin, etc.) (assembly of the list of iron-related genes is further described in the Materials and Methods section). To investigate how iron-related genes can contribute to HGSOC development and progression, we identified significantly perturbed KEGG pathways in each dataset using GAGE, a gene set analysis method for pathway analysis (Khatri et al., 2012), since it is optimized for use with both small and large datasets (Luo et al., 2009). GAGE, like Gene-Set Enrichment Analysis by which it was inspired, is a pathway analysis method in the functional class scoring (FCS) category where perturbed pathways can be identified from cumulative changes in gene expression between genes in a pathway, even when individual genes may not be identified as differentially expressed.

In parallel, we identified significantly differentially expressed genes between subtypes of the datasets as shown in Table 1 using the Bioconductor/R limma package for differential gene expression, which fits linear models to gene expression data and uses an empirical Bayes method to estimate the variance in gene expression by “borrowing information” across genes (Ritchie et al., 2015; Smyth, 2004). We analyzed which pathways were significantly perturbed and included several differentially expressed iron-related genes, and how consistent this observation was among the datasets (Fig. 2). We overlaid the perturbed pathways identified using GAGE analysis with differentially expressed iron-related genes as a conservative method to identify perturbed pathways in which iron-related genes make a contribution to the alteration in pathway activity.

Summary of the algorithm to identify mutually perturbed pathways with substantial iron-associated gene involvement. DE, differentially expressed; FDR, false discovery rate; lfc, log-fold change.

We assayed manually curated pathways specific to humans, each with its own pathway identifier (e.g., hsa04930) from the KEGG (Kanehisa et al., 2017). As described in the Materials and Methods section, a pathway was considered to be mutually perturbed if it was perturbed in datasets [1]–[3], and had ≥2 differentially expressed iron-related genes with a log2-fold change of ≥|0.2|. This threshold, which keeps genes with fold change ≤0.87 and ≥1.15, removes genes with very low fold changes in expression, while still allowing the capture of genes with potentially low fold changes in expression, but that may correspond to a pattern of changes in gene expression within the pathway, which is in the spirit of the FCS pathway analysis methods.

Four such pathways were identified (Supplementary Table S2). Three pathways with the lowest cumulative ranking (Supplementary Table S3) were considered for further analysis. These are the peroxisome proliferator-activated receptors (PPAR) signaling pathway (hsa03320), mineral absorption (hsa04978), and metabolism of xenobiotics by cytochrome P450 (hsa00980). Possibly due to the small sample size, GAGE analysis of dataset [4] only gave 11 significantly perturbed pathways, none of which satisfied our mutually perturbed pathway criteria. Therefore, to find mutually perturbed pathways, we used only datasets [1]–[3].

The top three mutually perturbed pathways with iron-related differentially expressed genes are in Table 2 (a full list of significantly perturbed KEGG pathways, along with differentially expressed iron-related genes for datasets [1]–[3] are in Supplementary Tables S4–S7, respectively). Confirming the dependence of OVC on iron, we first note that mineral absorption, the KEGG pathway that contains the preponderance of genes of iron metabolism, is a mutually perturbed pathway.

Furthermore, for each dataset in Table 2, several genes involved in intracellular iron regulation from the mineral absorption pathway are differentially expressed, including the iron storage gene ferritin (FTH1, FTL), iron exporters hephaestin (HEPH) and ferroportin (SLC40A1), iron importer DMT1 (SLC11A2), iron-import stimulator CYBRD1, and heme oxygenase 1 (HMOX1), the inducible member of heme oxygenase that catabolyzes heme to biliverdin, CO, and ferrous iron (Fe2+). While the mineral absorption pathway is not significantly perturbed in dataset [4], the genes FTL, FTH1, HMOX1, and SLC11A2 are all differentially expressed. Moreover, Basuli et al. (2017) showed that FTtr displays an iron retention phenotype in comparison to FTstem.

We also note that TFRC, which is not in the mineral absorption pathway, is significantly increased in the cancer v. nonmalignant cell-type comparisons in datasets [3] and [4] (with fold change 2.52 and 1.82, respectively). A full list of fold changes in the iron-related genes in the mineral absorption pathway that are significantly differentially expressed in each dataset are found in Supplementary Table S8. These results are consistent with increased intracellular iron in these HGSOC datasets.

Pathways involving FA import and synthesis are perturbed in HGSOC and include differentially expressed iron-related genes

As discussed in the Introduction, cancer cells increase synthesis (and, in some cases import) of FAs for a variety of purposes, including membrane biosynthesis and generation of signaling molecules. Increasingly, FA synthesis pathways are being considered targets for novel cancer therapeutics (Currie et al., 2013; Röhrig and Schulze, 2016). Thus, the presence of the PPAR signaling pathway, which integrates a number of important FA metabolic and signaling processes in the list of mutually perturbed pathways, indicates that there may be an important, and as yet unappreciated, connection between HGSOC-mediated changes in FA metabolism and iron-related genes.

PPARs are nuclear receptors that are activated by by-products of lipid metabolism and act as transcription factors for, and regulators of, key enzymes and proteins for lipid and carbohydrate metabolism, among other pathways. The three subtypes of PPARs are PPARα, PPARγ, and PPARβ/δ, which are encoded by three separate genes (Schoonjans et al., 1996; Tyagi et al., 2011). The KEGG PPAR signaling pathway (hsa03320) mainly includes the transcriptional targets of PPAR proteins, which are genes involved in various aspects of FA metabolism. Targets include SCD, SCD5, FADS1, and FADS2, which are iron-related genes that participate in FA synthesis and metabolism pathways by catalyzing desaturation of FAs to monounsaturated (SCD and SCD5) and polyunsaturated (FADS1and FADS2) FAs.

A further analysis of the datasets (and, the set of iron-related genes) provides additional evidence that FA pathways are altered in HGSOC, and that iron and iron-related genes may play an important role in these alterations. Critical genes involved in de novo FA synthesis, desaturation, and uptake, as well as arachidonic acid (AA) metabolism, the pentose phosphate pathway, the NADPH oxidase (NOX) family, and glycerolipid metabolism were found to be differentially expressed in all the datasets (Table 3), indicating that a perturbation in FA metabolism, with concomitant involvement of iron-related genes, had occurred in all the datasets.

3.

Fold Changes in Datasets [1]–[4] for Significantly Differentially Expressed (FDR Adjusted

	Dataset					Dataset
	[1]	[2]	[3]	[4]		[1]	[2]	[3]	[4]
FA synthesis and regulation					Pentose phosphate pathway
PPARG	1.15	1.33	1.78	3.12	G6PD				2.36
SREBF1				1.49	PGD	1.36		2.79	2.13
ACLY	1.23	1.41		2.21
FASN				1.83
ACSLx	‡‡	‡‡	‡‡	‡‡
FA desaturation					NOX proteins
FADS2•	1.22	1.34	1.48		CYBB• (NOX2)	1.95	1.66	1.82
SCD•	1.54	1.72	6.85	1.34	NOX4•	2.10	3.09	7.49	0.66
ELOVL2				2.65	NOX5•				2.01
ELOVL5	1.38	1.56			DUOX1•	0.80		0.12	0.15
					DUOX2•				0.51
AA metabolism					Glycerolipid metabolism
PLA2G4A			0.09	3.58	AGPATx	‡	‡	‡	‡‡
CYPx•	‡	‡‡	‡‡	‡‡	DGAT1			1.42	0.41
ALOXx•	‡		‡	‡‡	LPIN1	0.85		3.11	6.15
PTGS2•(COX-2)	1.26	1.63	0.10	2.39	LPIN2		1.20	2.10	3.33
					GK			0.42	3.74
					GPAM		1.56	2.92	1.85
FA uptake
LPL	2.00	2.08		20.41
CD36	2.80	4.62
FABP4	6.50	9.81
APOC2				0.75

Iron-related genes are indicated by a () following their name, genes that are part of a family are labeled with an “x” at the end (for example ACSL1-6 comprise ACSLx), and if one gene of that family is differentially expressed in a dataset, then one ‡ is used, if ≥2 genes of that family are differentially expressed in that dataset, then two ‡‡ are used. Unmarked cells indicate there were no significant differences in gene expression in that dataset.

AA, arachidonic acid; FA, fatty acid; NOX, NADPH oxidase.

In addition, as previously discussed, GAGE analysis of perturbed networks in dataset [4] generated far fewer perturbed pathways, but glycerolipid metabolism (hsa:00561) was found to be perturbed near the significance level cutoff (adjusted p value of 0.071 and rank order 13 of 222 human KEGG pathways). This is a pathway that links FA uptake and synthesis with synthesis of phospholipids (Prentki and Madiraju, 2008). Also, several critical genes involved in FA uptake (LPL), synthesis (FASN, SREBP1, and ACLY), desaturation (SCD and ELOVL2), and AA metabolism (PLA2G4A, PTGS2) were also found to be upregulated in the FTtr v. FTstem cells in [4], indicating that FA metabolism is perturbed in the model HGSOC stem cells (Table 3).

Finally, an ORA for iron-related genes in KEGG pathways showed that iron-related genes are significantly overrepresented in several FA metabolic and signaling pathways. ORA uses a statistical test, based most commonly on a hypergeometric distribution, to ascertain whether a set of genes is significantly overrepresented in a select pathway or pathways of interest (Khatri et al., 2012; Zhang et al., 2005a).

Notably, AA metabolism is ranked first in this analysis (iron-binding CYP and ALOX proteins are involved in AA metabolism and oxidation), see Supplementary Table S9 for full results. We also note that the KEGG pathways for mineral absorption, metabolism of xenobiotics by cytochrome P450, and lysine degradation are also observed in the ORA analysis, suggesting further that all these pathways have a potentially therapeutically tractable dependence on iron through the activity of iron-related genes.

A network for iron involvement in FA metabolism in OVC

To synthesize these results into actionable hypotheses, we constructed a network for involvement of iron-related genes in FA metabolism in OVC (Fig. 3). The network uses information from the preceding analysis by incorporating the key players in FA synthesis and import that have been especially implicated in cancer, several of which are iron-related genes and/or were differentially expressed in the datasets, and many of which are strongly linked to the PPARγ signaling pathway. The network is also deeply rooted in the known biochemistry of FA metabolism, and all links are supported by the literature (Supplementary Table S10). The value of the network lies in the ability to derive novel, experimentally verifiable hypotheses from it, which we do in the following subsections.

The network shows how both iron-binding (red outline) and other enzymes (orange fill) act on exogenous and de novo synthesized FAs in HGSOC. Cellular processes are labeled in green, the dashed line separates the exterior of the cell (top) from the interior (bottom). Relevant connections to specific hypotheses deriving from this network are shaded blue, red, and purple. FA, fatty acid; TAG, triacylglycerol; LPL, lipoprotein lipase; LA, linoleic acid; AA, arachidonic acid; MUFA, monounsaturated FA; PL, phospholipid; LPA, lysophosphatidic acid; FA-coA, fatty acyl-CoA; Lipid-OOH, lipid peroxide.

PPARγ- mediated promotion of de novo FA synthesis

It has been previously observed that expression and nuclear localization of the ligand-activated nuclear hormone receptor PPARγ is increased with increasing tumor grade in ovarian carcinoma (Zhang et al., 2005b). It is upregulated in all the datasets analyzed (Table 3), whereas PPARD and PPARA are not consistently significantly differentially expressed in more than any one dataset (not shown).

PPARγ upregulates lipogenesis and FA uptake, as well as inducing differentiation and restricting cellular progression through altering (by expression or activity) proteins involved in the cell cycle and inflammation (Michalik et al., 2004). Notably, despite the promotion by PPARγ of de novo FA synthesis, activation of PPARγ by agonist drugs can initiate cell apoptosis and other antitumorigenic processes (Vignati et al., 2006). In prostate cancer cells, it has been shown that the paradox may be due to a dose dependency: at low doses of an agonist drug, PPARγ activation results in upregulation of several lipogenic genes, whereas at higher doses, this effect disappears with a corresponding increase in cell apoptosis (Mansour et al., 2011).

In addition, the cumulative effect of PPARγ (protumorigenic vs. antitumorigenic) also depends on the molecular milieu of the given cell or tumor type, since PPARγ activity is regulated by several ligands resulting from the metabolism of AA and linoleic acid (LA), as well as kinase activity by MAPK, EGF, and PDGF, whose proportional concentrations are often altered in cancerous cells, and in turn regulates pathways such as TP53, Wnt, AKT/PTEN, and others, which may have mutated, upregulated, or downregulated components in a given cellular context (Dong, 2013; Robbins and Nie, 2012; Varga et al., 2011).

We consider that in HGSOC, PPARG is upregulated and activated to support FA import and synthesis, but strong pharmaceutical or endogenous activation will lead to PPARγ-dependent activation of proapoptotic and antitumorigenic pathways. The latter process is shown in Figure 3 using a thicker arrow connecting PPARγ and antitumorigenic processes. We observe that there exists a positive feedback cycle involving PPARγ: PPARG upregulation and activation lead to FA import and synthesis, which increase AA production and metabolism. The metabolites of AA are PPARγ agonists (Krey et al., 1997), thus completing the loop. We therefore hypothesize that iron-related genes involved in AA metabolism promote PPARγ activation, and thus FA synthesis and import processes in HGSOC cells (Fig. 3, relevant connections in blue).

Indeed, iron-related genes, in addition to being downstream effectors of PPARγ activation (such as SCD), are also critically involved in PPARγ activity since they comprise the predominant enzymes in AA metabolism. For example, a well-characterized activating ligand of PPARγ is 15d-δ12,14-PGJ2, a prostaglandin metabolite derived from AA through reactions mediated by the heme-binding enzymes PTGS1 (COX-1) and PTGS2 (COX-2). When fibroblasts are supplemented with this metabolite, they will undergo adipocyte differentiation and upregulate lipogenic pathways (Kliewer et al., 1995).

While data are not available for cancer cells on iron-mediated increases in eicosanoid (AA metabolite) production, Mattera et al. (2001) showed that iron-overloaded cardiomyocytes produced increased eicosanoids when stimulated by cytokines in comparison to control, and the increase in prostaglandins was due to iron-mediated COX-2 induction, showing support for our hypothesis.

NOX-mediated reactive oxygen species production

NOX proteins are heme-binding proteins that generate reactive oxygen species (ROS) from NADPH and have been found to be upregulated in several cancers, including ovarian and breast (Block and Gorin, 2012). In cancer cells undergoing de novo FA synthesis, NADPH production is increased by the pentose phosphate pathway (PPP), which generates ribonucleotides and NADPH from glucose. The PPP is accelerated in cancer cells through increased expression and activity of several enzymes of the PPP, including G6PD and PGD, although other mechanisms for PPP acceleration exist (Jiang et al., 2014; Patra and Hay, 2014). Key PPP and NOX genes are upregulated in the datasets (Table 3).

We thus derive the hypothesis that an increase in NADPH production through the pentose phosphate pathway, partially for the purpose of de novo FA synthesis, also provides more substrate for heme-binding NOX proteins, which are increased in HGSOC (Fig. 3, relevant connections in red).

In support of this hypothesis, we consider that NOX activity vis-a-vis iron has been most thoroughly investigated in the original phagocytic NADPH oxidase, which includes NOX2 as part of its complex. Inhibition of heme synthesis in a cell model of differentiating myeloid cells results in downregulation of NOX2 and disruption of NADPH oxidase activity compared to control (Yu et al., 1997). At the organismal scale, the iron-chelator desferrioxamine (DFO), but not iron-loaded DFO, was shown to suppress inflammation-induced NADPH oxidase activity in mouse dorsal air pouches as well as suppress expression of another heme-containing member of the complex, p22phox (levels of NOX2 were not measured) (Li and Frei, 2006).

In the neutrophils of iron-deficient anemic patients, NADPH oxidase activity was reduced relative to controls, and increased after iron supplementations (Kurtoglu et al., 2003). Since all NOX homologs share conserved structural and functional properties, the studies presented here suggest that the impact of iron levels on NOX family activity may be significant in HGSOC.

In addition, both the PPP and NOX proteins act as positive regulators of ferroptosis (Dixon et al., 2012). Therefore, an increase in activity of these pathways/enzymes in OVC cells may influence their susceptibility to this form of cell death. Indeed, recent work has shown that FTtr cells are more susceptible to ferroptosis than FTstem cells (Basuli et al., 2017), and as shown in Table 3, the FTtr cells also show upregulation of the PPP enzymes G6PD and PGD.

FA desaturases, FADS2 and SCD

FADS2, FA desaturase 2, also known as delta-6-desaturase controls the rate-limiting step in the conversion of the omega-6 LA to AA and the omega-3 alpha LA to eicosapentaenoic acid. FADS2 has a conserved heme-binding domain at its N-terminus that has been shown to be essential for its function (Pereira et al., 2003; Sayanova et al., 1999).

FADS2 activity has been shown to be increased in breast cancer (Pender-Cudlip et al., 2013) by measurement of the ratio of LA metabolites to AA. FADS2 was also found to be increased in melanoma and lung tumors implanted in mice, and its inhibition led to decreased tumor growth in both models (He et al., 2012). Similarly, inhibition of FADS2 in a mouse model of intestinal cancer inhibited tumorigenesis, an effect that was rescued by the addition of exogenous AA (Hansen-Petrik et al., 2002).

FADS2 and the downstream gene for the enzyme ELOVL5 are increased in the Mes subtype over the Diff subtype in datasets [1] and [2]. In these two datasets, we also observe increases in the genes encoding FA import proteins CD36 and FABP4 in the Mes subtype (Table 3). If Mes tumor cells, potentially due to influence and/or availability of extra free FAs from the surrounding stroma, import more free FAs, then they will have more LA substrate to convert to AA by FADS2 and its downstream enzymes.

Indeed, co-culture with adipocytes results in an increased FA import and FABP4 expression in OVC cells, and FABP4 is increased in expression in metastases to the omentum (which is composed mainly of adipocytes) over primary ovarian tumors (Nieman et al., 2011). While the Mes subtype was not assayed for increased adipocyte infiltration, increased FA uptake has also been observed in cancer cells subjected to hypoxic conditions (Kamphorst et al., 2013), indicating that a change in the tumor microenvironment can influence FA uptake in cancer cells.

Although we do not observe an increased FADS2 expression in dataset [4], an increased metabolism of AA (for example by COX-2) will position FADS2 as one critical contributor (in addition to PLA2G4A-mediated phospholipid cleavage) to AA-mediated proinflammatory responses by providing the substrate for this signaling cascade. Note that the LA-AA-eicosanoid pathway involves iron-binding proteins at each step, indicating that iron chelation can abrogate or alter its activity.

SCD is another iron-binding desaturase and its expression is increased in all the cell lines (Table 3). It catalyzes desaturation of the saturated FAs stearic acid and palmitic acid to form oleic acid and palmitoleic acid, respectively, which are both monounsaturated FAs (MUFAs). The expression level of SCD has been positively correlated with several different cancers, due to the protumorigenic properties of MUFAs (Igal, 2010).

We hypothesize that the activity of the protumorigenic and inflammatory pathways stimulated by AA and monounsaturated FAs are partially mediated by iron-related genes, and hence modulation of iron levels will reduce the activity of these pathways (Fig. 3, relevant connections in purple). While there have not been studies performed on the effect of iron chelation on FADS2 or SCD activity in cancer, rats fed an iron-deficient diet were found to have a decreased AA to LA ratio in their hepatic phospholipids, a result that the authors attributed to decreased FADS2 activity (Stangl and Kirchgessner, 1998). In the same study, the authors also found lower palmitic and oleic acid hepatic phospholipid proportions and higher stearic acid proportions, indicating decreased SCD activity.

Activity of SCD with respect to iron has been more extensively studied than FADS2: expression and activity have been shown to be increased in liver and muscle during iron overload (Pigeon et al., 2001; Rodriguez et al., 2007); it was downregulated in livers of rats fed an iron-deficient diet (Kamei et al., 2010), but it was downregulated in heart after iron overload (Rodriguez et al., 2007). An experimental investigation into FADS2 and SCD activity in OVC, and how it is impacted by iron chelation, would be an important first test in examining the connection between iron and the FA desaturation pathways in OVC.

Discussion

In this work, we have synthesized heterogeneous data representing high-grade serous ovarian cancer (HGSOC) to develop a network focusing on novel connections between perturbed pathways in FA import and synthesis and iron, which can serve to guide experimentalists and modelers in developing and testing hypotheses related to the involvement of iron and iron-related genes in perturbation of FA metabolism in HGSOC (Fig. 3).

Importantly, the network was developed using both HGSOC-specific data and literature, as well as literature from sources outside of OVC. Thus, it is possible that a subset of these results and hypotheses will transfer to other cancers. Indeed, Hall et al. (2016) developed a transgenic mouse model that showed that activation of the Myc proto-oncogene in lung adenocarcinoma resulted in an increase in AA and its derivatives (eicosanoids) with concomitant activation of the AA metabolizing COX and LOX enzymes, which are iron-related genes that we have been identified as playing a role in increased activity of the AA pathway in our network (Fig. 3).

The authors showed that this effect was partly due to Myc-mediated increase in cytosolic phospholipase (PLA2G4A) activity, and expression of PTGS2 (COX-2) and ALOX5 (also known as 5-LOX). When the COX and 5-LOX pathways were blocked, tumor growth was reduced, as predicted by our network. Notably, of the four datasets in our study, PLA2G4A and PTGS2 were most significantly increased in the Myc-transformed HGSOC stem cell line from dataset [4] (Table 3), indicating that a similar process may occur in these cells.

AA-derived eicosanoids play a role in two of the hypotheses generated by our network: the existence of a positive feedback loop between PPARγ and AA metabolism, and the activation of the LA-AA-EA pathway. FADS2 and the COX/LOX enzymes rely on iron to perform their enzymatic activities, and thus targeting the iron dependence of these proteins by chelation may be a mechanism to inhibit these pathways at multiple sites. Importantly, the Hall et al. (2016) result shows that the protumorigenic relationship between iron-binding proteins and FA metabolism may be a phenomenon that extends to a broader subset of cancers.

Due to the inherent difficulties in analysis, it is not common that heterogeneous data are synthesized as in this study. Nevertheless, we felt it an important goal to work toward, since each dataset has both positive aspects and drawbacks. Datasets [1] and [2] have large sample sizes, but suffer from not including data from normal tissue or distinguishing between different cell types.

Indeed, while the Diff subtype expresses markers of ovarian tumors known to be expressed by the tumors (MUC1 and MUC16) as well as the secretory fallopian tube marker SLP1, indicating that the majority of the cells in the Diff subtype represent bonafide tumor cells, higher expression of markers for myofibroblasts (FAP), and microvascular pericytes (ANGPTL2 and ANGPTL1) in comparison to the other subtypes suggest that there is a stromal infiltrate in the Mes subtype (Cancer Genome Atlas Research Network, 2011).

The consistency among all the datasets in the changes in expression of iron-related genes and the FA synthesis network components indicates that the relationship between the two may extend to cells of the tumor microenvironment. Dataset [3] has a smaller sample size, but the laser microdissection performed on the cells ensures that only epithelial cells have been analyzed. While dataset [4] has the smallest sample size, as an in vitro model, it is directly amenable to experimental validation.

Since each dataset represents a different aspect of HGSOC, an overlap in results between the different datasets may be considered a measure of their robustness. We have attempted to minimize heterogeneity by using data from similar gene expression platforms (Affymetrix HGU-133A or -133plus2), which allowed us to develop an analysis pipeline that was identical for all data.

Mutual perturbation of the metabolism of xenobiotics by cytochrome P450 (hsa00980) in the datasets indicates that cytochrome p450 enzymes, which are heme proteins involved in oxidation reactions, are altered in HGSOC. Since these enzymes are involved in drug metabolism, their alterations in cancer have been proposed to promote drug resistance (Rodriguez-Antona and Ingelman-Sundberg, 2006). Thus, it will be of interest to study how their activity may be impacted both by increased intracellular iron found in cancer cells and by iron chelation therapy.

The final mutually perturbed pathway, lysine degradation (hsa00310), showed changes in expression of the PLOD genes (Table S2). PLOD1, PLOD2, and PLOD3 are lysyl hydroxylases that require iron as a cofactor and are critical for collagen crosslinking, and can contribute to extracellular matrix deposition during cancer development and progression (Xiong and Xu, 2016). Thus, a further investigation of the impact of iron on PLOD activity may lead to a connection between iron levels and extracellular matrix (ECM) development, which may be of special interest since the PLOD genes are upregulated by hypoxia-inducible factor 1 (Gilkes et al., 2014), which also plays a central role in intracellular iron homeostasis (Peyssonnaux et al., 2008).

There are a number of directions in which the network developed in Figure 3 can be further expanded. We provide some potential avenues here. FA ß-oxidation (FAO), which occurs in the mitochondria, has also been implicated in various cancers (Carracedo et al., 2013). The classical consideration that FAO and FA synthesis are mutually exclusive is giving way to a more nuanced view that the processes may contribute to the other's activity, and/or act independently (Beloribi-Djefaflia et al., 2016; Caro et al., 2012; Carracedo et al., 2013). Since FAO is regulated by PPAR-associated pathways (Varga et al., 2011), one may expand the network in Figure 3 to include FAO pathways.

Interestingly, carnitine, an ammonium compound that is essential for the transport of FA into the mitochondria for FAO (through the enzyme CPT1, carnitine palmitoyltransferase 1), can provide a potential link between iron and FAO. Carnitine can be absorbed from the diet, or synthesized in cells from lysine (one of the mutually perturbed pathways, Table S2, and one of the top 25 KEGG pathways in which iron-related genes are overrepresented, Table S9). Indeed, two enzymes critical for its synthesis, BBOX1 and TMLHE, use iron as a cofactor (Vaz and Wanders, 2002). Therefore, iron may directly impact the activity of these enzymes, and thus the availability of carnitine for FA transport. While we did not observe consistent changes in BBOX1 or TMLHE in the datasets, the impact of iron on their activity, and thus the activity of CPT1, may be significant.

More generally, the results established herein may be furthered to identify connections between iron-related genes and perturbed pathways not considered in this work. For example, one may look at potential relationships between the mutually perturbed pathways—since they share the similarity of being perturbed in multiple datasets, and each contain several differentially expressed iron-related genes.

For example, it has been shown that human vascular endothelial cells respond to PPAR ligand treatment by upregulating heme oxygenase 1 (HMOX1), an essential enzyme in heme catabolism that is part of the KEGG mineral absorption pathway. HMOX1 was found to have two PPAR-responsive elements in its promoter (Krönke et al., 2007). Appropriately, HMOX1 is significantly upregulated (data not shown) in all the dataset comparisons performed. This indicates that, while we have discussed the influence of iron and iron-related proteins on PPARγ activity, PPARγ and other PPARs may also act to alter iron metabolism.

Conclusion

We have identified several novel connections between iron, iron-related genes, FA synthesis and import, and HGSOC. With the development of a network incorporating all these factors (Fig. 3), as well as several hypotheses derived from it, we aim to stimulate inquiry into the interdependence between iron and FA metabolism HGSOC and other cancers. This will enhance not only understanding of the basic biology of HGSOC but also potentially identify vulnerabilities that can be exploited therapeutically.