Spheroid growth in ovarian cancer alters transcriptome responses for stress pathways and epigenetic responses

Ovarian cancer is the most lethal gynecological cancer, with over 200,000 women diagnosed each year and over half of those cases leading to death. These poor statistics are related to a lack of early symptoms and inadequate screening techniques. This results in the cancer going undetected until later stages when the tumor has metastasized through a process that requires the epithelial to mesenchymal transition (EMT). In lieu of traditional monolayer cell culture, EMT and cancer progression in general is best characterized through the use of 3D spheroid models. In this study, we examine gene expression changes through microarray analysis in spheroid versus monolayer ovarian cancer cells treated with TGFβ to induce EMT. Transcripts that included Coiled-Coil Domain Containing 80 (CCDC80), Solute Carrier Family 6 (Neutral Amino Acid Transporter), Member 15 (SLC6A15), Semaphorin 3E (SEMA3E) and PIF1 5'-To-3' DNA Helicase (PIF1) were downregulated more than 10-fold in the 3D cells while Inhibitor Of DNA Binding 2, HLH Protein (ID2), Regulator Of Cell Cycle (RGCC), Protease, Serine 35 (PRSS35), and Aldo-Keto Reductase Family 1, Member C1 (AKR1C1) were increased more than 50-fold. Interestingly, EMT factors, stress responses and epigenetic processes were significantly affected by 3D growth. The heat shock response and the oxidative stress response were also identified as transcriptome responses that showed significant changes upon 3D growth. Subnetwork enrichment analysis revealed that DNA integrity (e.g. DNA damage, genetic instability, nucleotide excision repair, and the DNA damage checkpoint pathway) were altered in the 3D spheroid model. In addition, two epigenetic processes, DNA methylation and histone acetylation, were increased with 3D growth. These findings support the hypothesis that three dimensional ovarian cell culturing is physiologically different from its monolayer counterpart.


Introduction
Despite recent improvements in surgery and chemotherapy, ovarian cancer is still the leading cause of death from gynecological malignancy [1]. Due to poor detection methods and a lack of symptoms, most patients are diagnosed at advanced stages, when the tumor has metastasized and spread [2]. Studies suggest that in order for metastasis to occur, the cancer cells must undergo phenotypic changes modulated by the epithelial-mesenchymal transition (EMT) [3].
EMT is a distinct process whereby epithelial cells undergo changes in favor of mesenchymal properties [4]. This process is most commonly observed during developmental stages when epithelial cells must migrate and dedifferentiate, such as in the formation of the mesoderm during gastrulation [5]. In order to properly study this phenomenon, scientists have discovered multiple factors and signals which can induce EMT. Of these, the most popular inducer is transforming growth factor β (TGFβ) [6,7]. The addition of TGFβ to epithelial cells induces transient EMT within hours of treatment through activation of the Smad pathway [8].
Although two dimensional (monolayer) tissue culture models are largely used to study the EMT process, evidence suggests that three dimensional (spheroid) culturing may be more physiologically relevant as it better emulates oxygen levels, pH conditions, glucose levels, extracellular matrix strength, and overall morphology of in vivo solid tumors [9][10][11][12]. This is especially the case when focusing on metastasis, tissue invasion, angiogenesis, and drug sensitivity [13][14][15].
At least a third of ovarian cancer patients present with ascites [16]. Ascites is the accumulation of fluid in the peritoneal cavity which may contain ovarian cancer cells, lymphocytes, and mesothelial cells in the form of single cells and aggregates [17]. Further studies revealed that ascites spheroids may cause secondary tumors due to their ability to adhere to extracellular matrix proteins via interaction between multiple integrins and their ligands [18,19]. Here, we conducted a comprehensive gene expression analysis for the process of culturing HEY epithelial ovarian cancer cells in 3D vs. 2D cultures during the TGFβ-induced EMT process. Using subnetwork enrichment analysis, we identified stress pathways, DNA integrity pathways, and epigenetic processes as those most affected by 3D vs. 2D growth.

Cells culture and treatments
The HEY human ovarian cancer cell line (kindly provided by Dr. Meera Nanjundan, University of South Florida, Tampa, FL) was authenticated by short tandem repeat (STR) DNA profiling (Genetica, Inc.) and was compared to ATCC and previously published profiles [20]. Cells were cultured in a humidified incubator containing 5% CO 2 at 37˚C in RPMI with 10% fetal bovine serum and 1% Pen-Strep-Glutamine. Spheroid formation was accomplished through the hanging drop method. Briefly, trypsinized HEY cells were resuspended at 1 x 10 6 cells/mL in supplemented RPMI. Multiple 25 μl droplets of the cell solution were then placed onto plate lids, inverted, and incubated for 72 hours to allow cells to aggregate into spheroids. To induce EMT, a final concentration of 10 ng/μl TGFβ was added to cells and incubated for 72 hours as monolayer cell culture or hanging drops during the creation of spheroids. Following incubation, monolayer cells and spheroids were photographed using an Evos1 FL Cell Imaging microscope (ThermoFisher Scientific) and then collected in 1X PBS.
(Applied Biosystems) per the manufacturer's instructions. Subsequent samples were then diluted to 50 ng/ μL and used as a template for quantitative PCR (qPCR). qPCR was accomplished with a Step One Plus Real-time PCR system (Applied Biosystems) and SYBR1 Green Supermix with ROX (BioRad) according to the manufacturer's protocol. Relative mRNA levels were quantified for ahnak2, akr1c1, ccdc80, hspa1a, hsph1, prss35, rgs2, and rrad using genespecific primers (S1 Table).

Affymetrix GeneAtlas Platform and 3'IVT compatible U219 probe arrays
The oligonucleotide probe arrays used were the Affymetrix HG-U219 human array strips. These arrays consist of more than 530,000 probes detecting over 36,000 transcripts and variants-representing more than 20,000 genes.

Sample processing and validation
HEY Human ovarian carcinoma cells treated with TGFβ were cultured as 2-dimensional monolayers or as 3-dimensional spheroids using the hanging drop method. Replicates of 4 separate experiments were performed on 4 different culturing days. Cells were then harvested and total RNA was isolated from all 8 samples using Qiagen's RNeasy Mini Kit according to the manufacturer's instructions. RNA integrity was verified with 1 μl of total RNA on RNA 6000 nanochips (Agilent) using The Agilent 2100 Bioanalyzer. All results for all 8 samples reported RNA Integrity Numbers (RIN) > 9. 100 ng of polyadenylated RNA was converted to cDNA and then amplified and labeled with biotin using Affymetrix 3'IVT Expression System according to the manufacturer's instructions. Hybridization with the biotin-labeled RNA, staining, and scanning of the chips followed the procedure outlined in the Affymetrix technical manual. You can put this at the end of the methods on the micro array hybridization. All microarray data have been deposited into Gene Expression Omnibus (GSE80373).

Data analysis
Scanned output files were visually inspected for hybridization artifacts and then annotated and normalized using Affymetrix Expression Console v1.3.1. Additionally, all QC metrics reported no outlier samples. Signal intensity was scaled to an average intensity of 500 during comparison analysis. Annotated expression data were assigned an ANOVA P value for the likelihood that any perceived difference was due to chance. The P values for all probe sets were exported to a text file and all pairwise comparisons of Bi-weight average signals were then aligned in MS Excel. For the comprehensive analysis, P < 0.05 was identified as having a linear change (increased or decreased) for the comparison. This analysis reported~3,329 genes that showed differential regulation greater or less than 2-fold. To increase stringency, only genes whose fold-change was greater than or less than four (p < 0.05) were considered to be differentially expressed, reducing the pool to 493 candidate genes.
Data were analyzed using different visualization techniques that included hierarchical clustering, principle component analysis (PCA), and volcano plots. All analyses were conducted in JMP Genomics v7.0. Transcripts that showed p<0.01 were used in the two-way clustering using the Fast ward algorithm after each row was centered to a mean of zero (0) and variance scaled to one. Results from hierarchical cluster analyses were visualized using heat map dendrograms, and biological replicates were partitioned into groups based on similarity (i.e., individuals clustered together are most similar). PCA used spatially weighted averages and were conducted on normalized expression values, the same dataset used in the cluster analysis. Volcano plots were generated in JMP Genomics 7.0 following an ANOVA to identify differentially expressed transcripts.

Pathway Studio analysis
Pathway Studio 9.0 (Elsevier) and ResNet 10.0 were used for sub-network enrichment analysis (SNEA) of cell processes [21]. The option of "best p value, highest magnitude fold change" in Pathway Studio was used for duplicated probes. Transcripts were successfully mapped using GeneBank ID. SNEA was performed to identify gene networks that were significantly different. A Kolmogorov-Smirnov test with 1000 permutations was conducted to determine whether certain networks were preferentially regulated compared to the background reference probability distribution. Networks were constructed based on common regulators of expression and regulators of specific cell processes. The enrichment P-value for a gene seed was set at P < 0.05. Additional details on the use of SNEA can be found in Langlois and Martyniuk [22].

Results and discussion
HEY cells treated with TGFβ have distinct gene expression profiles when grown as 3D spheroids vs. 2D monolayers To examine whether growth of HEY ovarian cancer cells as 3D spheroids vs. 2D monolayers could influence gene expression, we performed microarray analysis on four biological replicates of cells grown under each condition. To create spheroids, the hanging drop method [23] was used as outlined (Fig 1). We find that HEY cells treated with TGFβ have distinct gene expression profiles, depending on whether they are cultured under 2D vs. 3D conditions ( Fig  2). Hierarchical clustering reveals that each of the two groups form into distinct clades based on gene expression (Fig 2A). Additionally, principle component analysis shows that biological quadruplicates for HEY 2D and HEY 3D samples are more similar to themselves than to each   other ( Fig 2B). Volcano plot analysis shows that there are numerous genes between the two groups that have a high fold change and that are also statistically significant (S1 Fig).
RT-qPCR was performed on eight different genes which were found to be up or downregulated by varying degrees upon 3D culturing to validate our microarray analysis (Fig 3). The Variability in transcriptome response separates strongly along the first PCA1. Red color is HEY 3D biological replicates and green color is HEY 2D biological replicates. The four biological replicates for HEY 2D are more closely related in expression as compared to HEY 3D.
https://doi.org/10.1371/journal.pone.0182930.g002 Altered gene expression upon spheroid growth in ovarian cancer eight genes we chose to validate include two highly upregulated genes, the serine protease PRSS35 and the aldo-keto reductase AKR1C1, as well as the highly downregulated E-cadherin regulator CCDC80. Additionally, we chose some genes that were more moderately regulated, including the chaperones HSP110 and HSP70, the G-protein signaling inhibitor RGS2, the fibroblast growth factor secretion regulator AHNAK2 and the apoptotic inducer RRAD. Our results show that similar trends are observed via RT-qPCR for these genes as compared to our microarray results. Overall, our data show that HEY ovarian cancer cells, treated with TGFβ to induce EMT, show dramatic differences in gene expression profiles when grown as 3D spheroids vs. standard 2D culture.
Transcripts affected by 3D culturing may enhance the tumorigenicity of HEY cells Growth as 3D spheroids leads to both the upregulation and downregulation of gene expression (Table 1. For a complete list of gene expression data see S2 Table). The most highly upregulated gene under 3D growth conditions in HEY cells is AKR1C1 (Aldo-Keto Reductase Family 1, Member C1), part of a family of cytosolic NADP(H)-dependent oxidoreductases that is involved in detoxification of xenobiotics, steroids, and polycyclic aromatic hydrocarbons [24]. Increased expression of AKR1C1 is associated with the development of cisplatin resistance in human ovarian carcinoma cells [25]. Previous work shows that ovarian cancer cells, grown in 3D, are more resistant to chemotherapy treatment [26]. It is therefore plausible that the increased expression of AKR1C1 upon ovarian cancer cell 3D growth may promote chemotherapy resistance. PRSS35 (Protease, Serine 35), another highly upregulated gene upon 3D growth, belongs to the trypsin class of serine proteases. This protease is highly expressed in the mouse ovary [27]. Proteases play an important role in proteolysis that is essential for tissue remodeling that occurs during EMT by breaking down the extracellular matrix. While metalloproteases play the largest role in this process, serine proteases also contribute [28]. We postulate that enhanced expression of PRSS35 under 3D growth conditions may promote ovarian EMT through the proteolytic digestion of cell attachments. A transcript that was highly downregulated upon 3D growth is CCDC80 (Coiled-Coil Domain Containing 80). Interestingly, CCDC80 null mice develop thyroid adenomas and ovarian carcinomas, and CCDC80 gene expression is decreased in human ovarian cancer samples as compared to normal ovarian samples [29]. The mechanism for the tumor suppressor properties of CCDC80 may be through its ability to induce E-cadherin expression [29]. E-cadherin provides crucial cell-cell adhesion to hold epithelial cells tightly together, and enhanced E-cadherin expression promotes the epithelial phenotype over the cancer-promoting mesenchymal phenotype [30]. Collectively, we propose that regulation of these and other genes upon 3D growth may enhance the tumorigenic properties of ovarian cancer cells.
Gene-set enrichment analysis reveals novel pathways that are enriched upon HEY cell growth as 3D spheroids Gene set enrichment analysis suggested that pathways downregulated by more than 15% include branched chain amino acid metabolism, folate biosynthesis, fatty acid oxidation, and the mevalonate pathway, while pathways upregulated in the 3D cells include those related to tumor necrosis factor receptor (TNFR) superfamily member 1A and 6 (for a complete list of gene set enrichment pathways identified see S3 Table). These pathways were increased more than 20% with 3D growth in HEY cells, and they have all been related to cancer growth. Branched chain amino acids, such as leucine, positively regulate the mammalian-target-ofrapamycin (mTOR) pathway [31]. The mTOR pathway is involved in regulating cellular functions involved in growth and proliferation, and upregulation of this pathway is commonly observed in human cancers [32]. Folate plays a role in nucleotide synthesis and methylation, making it essential to rapidly growing cancer cells [33]. Fatty acid oxidation is required for Transcripts are organized as those down-regulated and up-regulated between the two groups. The gene, as well as the role of the protein in the cell is provided in addition to the relative fold change and p-value. For a complete list of gene expression changes see S2 Table. https://doi.org/10.1371/journal.pone.0182930.t001 functional angiogenesis and is utilized by cancer cells to overcome metabolic stress to proliferate [34]. The mevalonate pathway is responsible for converting acetyl-coenzyme A into isoprenoids, which are required for cholesterol and steroid synthesis [35]. LDL-cholesterol accumulates in cancer tissues and is associated with migration, proliferation, and loss of adhesion from the primary tumor which are vital steps in the epithelial to mesenchymal transition [35]. Upregulation of TNFR family members allows signaling via the NF-κB transcription factor, which can promote tumorigenesis through the activation of the expression of genes involved in processes including cell proliferation, migration and anti-apoptosis [36]. Therefore, we conclude that multiple cellular pathways that are altered upon 3D growth could increase tumorigenicity.

EMT genes are altered upon spheroid culturing in ovarian cancer
Unique sub-networks underlie the transition from 2D to 3D growth in HEY cells (Table 2, for a complete list of enriched cell networks see S3 Table). We note that many of the genes involved in EMT are differentially expressed upon formation of spheroids as compared to growth in monolayer (Fig 4). EMT transcription factors downregulate the expression of cellcell adhesion proteins and upregulate metastatic proteins, leading to cancer progression [37]. Importantly, a number of EMT transcription factors are enhanced upon growth in 3D, including Snail 1 (SNAI1), Snail 2 (SNAI2), ZEB2, TCF3, and SIX1 (S4 Table). Supporting this finding, in recently published work, we show via qPCR that EMT transcription factors, while slightly upregulated upon TGFβ treatment under 2D growth conditions, are upregulated more significantly when grown in 3D [38]. Thus, 3D growth, through enhanced expression of EMT transcription factors, is likely to allow for a stronger induction of EMT by TGFβ.

Transcriptional networks involving stress pathways are altered in HEY cell growth as 3D spheroids
Interestingly, a number of the identified sub-networks are related to the response to stress ( Fig  5, for a complete list of gene networks related to stress see S4 Table). The finding that stress responses are activated upon 3D growth is not surprising, given that growth under these conditions is likely to cause nutrient limitation and hypoxia, among other cellular stresses. Transcripts in the network that are related to the oxidative stress response include peroxiredoxin 2 (PRDX2) (Fold change = +2. Also changed in 3D cells are transcripts that are associated with the heat shock response. Some examples include heat shock protein 90kDa alpha (cytosolic), class A member 1 (HSP90AA1) (+1.16), heat shock 27kDa protein 1(HSPB1) (+1.20), and heat shock transcription factor 1 (HSF1) (+1.28). Recent studies have suggested that many cancer types utilize activation of heat shock factor 1 (HSF1), the master regulator of the heat shock response, to stabilize oncogenes within the cell [44]. The heat shock response is a highly conserved response to specific environmental stressors such as heat shock, heavy metals, and oxidative stress [45]. HSF1 promotes the transcription of heat shock protein genes, many of which encode molecular chaperones. Increased chaperone expression provides a cytoprotective response in advanced cancers with acidotic, hypoxic, and nutrient-deprived microenvironments [46]. Elevated levels of one or more chaperones are commonly found in both solid cancer tumors and hematological malignancies [47][48][49][50][51]. In fact, overexpression of HSP27, HSP70, or HSP90 Table 2 correlates with poor prognosis and may contribute to drug resistance in some cancer types [52][53][54][55]. Activation of the heat shock response in spheroids may promote cell survival, in the face of multiple cellular stresses, through enhanced expression expression of HSF1 target genes, including molecular chaperones. Further analysis of this spheroid model could offer valuable insight into how HSF1 and HSPs could be utilized as targets for chemotherapy in ovarian cancer patients. Other examples of stress-related genes that show a high magnitude of response in terms of downregulation include the PIF1 5'-to-3' DNA helicase homolog (S. cerevisiae) (-13.0), epidermal growth factor receptor (-4.6), far upstream element (FUSE) binding protein (-4.4), and cysteine-rich, angiogenic inducer (-4.4). Increasing transcripts in the network included nuclear factor of kappa light polypeptide gene enhancer in B-cells inhibitor, alpha (+10.3), FBJ murine osteosarcoma viral oncogene homolog (+14.7), prostaglandin-endoperoxide synthase 2 (prostaglandin G/H synthase and cyclooxygenase) (+29.5), and interleukin 8 (+49.8). Many of these aforementioned transcripts have been shown to play critical roles in the progression of cancer. For example, PIF1 is a DNA-dependent adenosine triphosphate (ATP)-metabolizing enzyme that is required for proper replication and repair during cell division. Studies show that this protein inhibits S-phase progression and reduces proliferation rates of RAS  Table. https://doi.org/10.1371/journal.pone.0182930.g004 oncogene-transformed fibroblasts [56], and may therefore be a novel drug target for cancer therapy. Here, the downregulation of PIF1 mRNA may lead to increased proliferative activity during the transition from a non-cancerous cell to a cancerous cell. Our network analysis also identifies DNA duplex unwinding (down regulated -1.49) and S-M checkpoint (down-regulated -1.59) as processes significantly down-regulated in 3D spheroids, processes which may be related to changes in PIF1 expression. In terms of ovarian cancer, there are many studies implicating epidermal growth factor receptor (EGFR) as a key regulator of cell differentiation Specifically, EGF-induced EMT increases phosphorylation of Akt, ERK1/2, and S6 ribosomal protein, which alters ovarian cancer cell proliferation and differentiation [57]. Additionally, elevated EGFR expression in tumor stroma is linked to aggressive epithelial ovarian cancer in patients and relates to Ki-67 expression in tumor cells [58]. In terms of up-regulated transcripts, notable transcripts included FBJ murine osteosarcoma viral oncogene homolog, a member of the Fos family of transcription factors. In cancer, this gene is a proto-oncogene implicated in cell proliferation and transformation [59]. Lastly of interest in the network was  Table. https://doi.org/10.1371/journal.pone.0182930.g005

Theme Gene Set Seed Total # of Neighbors # of Measured Neighbors Median change p-value
Altered gene expression upon spheroid growth in ovarian cancer IL8, a gene that showed a dramatic increase in expression of~50-fold in 3D cells. Polymorphisms in this gene have been shown to be associated with a significantly higher risk of ovarian cancer [60,61]. Studies also show that increases in IL8 is associated with increased tumor growth and metastases [62].
Various responses related to DNA integrity were altered, including the response to DNA damage, genetic instability, nucleotide excision repair, and the DNA damage checkpoint pathway. The hypoxic environment created within spheroids may lead to increased DNA damage. Spheroid culturing has shown to alter chromatin packaging which in turn improves DNA repair through what is known as the cell "contact" effect [63]. DNA damage via metabolic products and by-products, such as ROS, may decrease replication fidelity, resulting in increased mutagenesis. Mutations and chromosomal abnormalities can increase the risk of cancer through the activation of oncogenes or the inactivation of tumor suppressor genes. Complex DNA-damage response mechanisms have evolved in order to isolate and repair these mutations. One such mechanism is the nucleotide excision repair (NER) pathway. NER is able to eradicate a variety of DNA lesions due to its ability to circumvent recognition of the lesion itself and focus on a set of commonalities shared by many different lesions [40]. Our analysis showed a significant enrichment of the NER pathway when cells are cultured as 3D spheroids.
Sub-network enrichment analysis identifies that genes involved in epigenetic processes are enriched upon HEY cell growth as 3D spheroids Epigenetic alterations can play a key role in promoting transformation and tumor growth [64][65][66], although the underlying mechanisms are still being worked out. We find two epigenetic  Table. https://doi.org/10.1371/journal.pone.0182930.g006 Altered gene expression upon spheroid growth in ovarian cancer processes, DNA methylation and histone acetylation, to be enriched upon 3D growth. A number of genes are both up-and downregulated related to DNA methylation, and this process is affected by about 4-5% (Fig 6, for a complete list of all genes in the DNA methylation pathway see S4 Table). DNA methylation changes are seen in several cancer types and have been linked to changes in gene expression in highly metastatic tumors [67].
The process of histone acetylation is also significantly affected in both positive and negative fashions (Fig 7, for a complete list of all genes in the histone acetylation pathway see S4 Table). Histone acetylation has important roles in diverse processes including gene regulation, DNA damage repair, and DNA replication. Significant correlations between cancer patient survival and histone acetylation have been shown in several studies, although it has been associated with both better and worse prognoses depending on the specific cancer type [68,69]. Thus, we conclude that growth as 3D spheroids may affect the epigenetic status of cancer cells, which could alter the tumorigenic properties of the cells. Pathway for histone acetylation. A number of histone modifying enzymes were increased in expression in the 3D group, and this may be reflective of chromatin remodeling during cancer progression. Red indicates the gene or process is up-regulated and blue indicates the gene or process is downregulated. All genes for histone acetylation pathways are listed in S4 Table. https://doi.org/10.1371/journal.pone.0182930.g007

Conclusions
Spheroid formation more closely mimics that of an in vivo tumor due to the cells ability to form an extracellular matrix and cell adhesions. Here, using the HEY ovarian cancer cell line, we show that 3D growth affects a number of cellular processes, including the EMT process, multiple cellular stress pathways, DNA integrity pathways, and epigenetic pathways. As these pathways all could affect tumorigenesis and the response to chemotherapies, our studies suggest that using 3D culture instead of 2D monolayers may be more informative in studying the properties of ovarian cancer cell lines.