Discrete Molecular Classes of Ovarian Cancer Suggestive of Unique Mechanisms of Transformation and Metastases

Rapid metastases and disease progression are a major cause of mortality in high-grade serous epithelial ovarian cancer (HGSC) and with which, a mechanistic association of epithelial to mesenchymal transition (EMT) is implied (1). The process of EMT involves cooperation of niche factors through autocrine–paracrine signaling with a transcriptional program driven by E12/47 (TCF3), TCF4, Twist, ZEB1, ZEB2, Snail (SNAI1), Slug (SNAI2), FOXC1, or FOXC2 (refs. 2, 3; collectively termed as EMT-TFs). Besides EMT, passive shedding, cooperative migration of multicellular aggregates from the surface of growing tumors, and peritoneal seeding (4), and associations of lysophosphatidic acid (LPA) and/or sphingosine 1-phosphate (S1P) signaling with HGSC metastases (5) are reported. The latter are increasingly being associated with cancer cell survival, tumor progression, angiogenesis, metastasis, and chemoresistance (6).

The ovarian surface epithelium (OSE), an innocuous, primitive mesothelium belies the fact that its transformation results in aggressive HGSC. Incessant ovulation and release of the ovum into the fallopian tubes leads to localized OSE injury and consequent healing through EMT and mesenchymal to epithelial transition (MET; 7, 8). Recently, involvement of FT distal fimbrial epithelial (FTE) transformation in HGSC was suggested based on similarities with Müllerian tract tumors (9, 10). Evidences and contradictions for both cells-of-origin exist; (i) a pathologic association of OSE-derived inclusion cysts and inflammation in association with cyclic ovulation-wound healing with HGSC is established (8, 11). However, higher level of differentiation in tumor epithelia than OSE contradicts the dogma of tumor-associated maturation arrest. (ii) Oxidative stress, cell hyperactivation, and DNA damage in FTE lead to preneoplastic tubal intraepithelial lesions in p53/BRCA1/BRCA2 mutation harboring patients in whom prophylactic oophorectomy does not reduce risk of recurrence (12, 13); however, HGSC incidence in patients without these mutations remain unaccounted for. This conundrum suggests subsistence of at least 2 subtypes in HGSC that at present is diagnosed as a single disease based on clinical characteristics.

In the present study, we attempted tumor classification through two approaches; (i) supervised metastases-gene signature based clustering, and (ii) unsupervised, Weighted Gene Correlation Network Analysis (WGCNA) analyses. Both methods defined similar tumor groups in which molecular systems networks-based resolution of class-specific tumor-associated biologic functions revealed unexpected pathway associations. Class1 tumors survive, metastasize in an EMT-independent manner, and are associated with aberrant differentiation, genetic instability, and immune evasion. Their molecular profiles and cellular behavior suggest possible involvement of FTE and ovarian stroma. Class2 tumors express extracellular matrix (ECM)-EMT–driven invasive gene networks and IFN signaling that put forth a plausible cell-of-origin being the OSE. Preliminary validation of these suggestive etiologies was carried out in experimental models using distinctive markers purported to be associated with either the FTE (oxidative stress, DNA damage, genetic instability, and apoptosis) or OSE (EMT, ECM).

Level 3 high-grade ovarian serous cystadenocarcinoma gene expression data was extracted from TCGA (14; https://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp) from two platforms, viz. Affymetrix (512 tumor and 8 normal samples) and Agilent (359 tumor and 8 normal samples) as available on October 2010 (Supplementary Table S1). Data was Lowess normalized and log2 transformed. Further correlation analysis of EMT-TFs SNAI1, SNAI2, TCF3, TCF4, TWIST1, ZEB, and ZEB2 with epithelial and mesenchymal markers involved applying Pearson correlation coefficient (−0.5 < r < 0.5; P < 0.05) followed hierarchical clustering of both datasets. ChIP-on-chip was carried out as described earlier (15). Briefly, 10⁸ cells (serous ovarian adenocarcinoma cell line - A4 established earlier in our lab; refs. 15, 16) were chemically crosslinked, harvested, lysed, sonicated to solubilize, and shear crosslinked DNA that was immunoprecipitated with Snail, Slug, or Twist antibodies, purified, labeled, and hybridized with Agilent human whole-genome promoter microarrays. Identification of global targets of Twist, Snail, and Slug was carried out using Agilent Genomic Workbench Lite version6.5.0.18; using a significance cutoff (P < 0.05). Expression levels of these genes were extracted from GSE18054 (gene expression profiles of A4 pretransformed and transformed cells). Genes with fold change (ratio of log intensity value between transformed and A4 pretransformed) expression value between 0 and 1 were considered low expressed and those with negative values were considered repressed. Furthermore, Pearson correlation coefficient (−0.4 < r < 0.4 with at least 2 gene partners; P < 0.05) was carried out in the A4 gene expression dataset. Finally, the 99 gene set that represented common targets of all three transcription factors were considered as EMT target genes. Using a similar approach, Pearson correlation analyses of known LPA and S1P associated receptors as well as components of LPA-S1P receptor-mediated signaling (identified from REACTOME, KEGG and Pathway Interaction Databases was carried out (http://www.reactome.org/ReactomeGWT/entrypoint.html, http://www.genome.jp/kegg/pathway.html, http://pid.nci.nih.gov/PID/index.shtml). The final list of metastases genes were subjected to median absolute deviation (MAD) analysis for highly differential genes [criterion (−1 < fold change >1)] in at least 10% of samples) followed by Pearson correlation coefficient (−0.4 < r > 0.4 with at least 2 gene partners; P < 0.05) that identified 39 genes referred to as the 39 metastases-associated gene signature.

We applied this 39 metastases-associated gene signature toward tumor classification using three methods, (i) Silhouette plots that identify silhouette-width–based clusters and associated outliers, (ii) Consensus average linkage clustering, and (iii) Hierarchical clustering. The 39 metastases-associated gene signature was applied for K-means classification of the cohort (n = 359) from K = 2 to K = 6 and silhouette plots were generated in which the Silhouette width (indicating clustering strength) is determined for each sample (17). Agglomerative average linkage hierarchical consensus clustering (R: ConsensusClusterPlus) was performed on positive silhouette width samples (n = 329) using a distance metric (1- Pearson correlation coefficient); the procedure was run over 1,000 iterations and a subsampling ratio of 0.8. Consensus clustering involves subsampling from a set of total tumor samples through calculation of pairwise consensus values based on the probability that two samples occupy the same cluster out of the number of times they occurred in the same subsample (18). A consensus matrix was the summarized form of all pairwise consensus values that then cluster over number of iterations. Unsupervised Principle Component Analyses (PCA) was carried out for class resolution, whereas Between Group Analyses (BGA) was performed to identify class-specific gene–sample association. Common tumor samples from the two methods were further subjected to hierarchical clustering.

For details, see Supplementary Data S1 and Refs. 19, 20.

For details, see Supplementary Data S2.

Gene set enrichment analysis (GSEA; ref. 21) was performed to identify enriched gene sets at P < 0.1 and q < 0.25. Class comparison networks were visualized using Enrichment map in Cytoscape_v2.8.2. Application of GSEA identified some biologic processes and signaling pathways (P < 0.05; FDR < 0.25).

Cells were cultured as per the American Type Culture Collection guidelines. Reverse transcriptase-PCR (RT-PCR) for highly co-relating genes of the top 25 list of each module was done as described earlier (1) using GAPDH as an internal control. Fold change > 1.0 was considered significant. Relative gene expression was plotted as a heatmap using MeV. Scratch assay was carried out as described (1); in brief, cells were grown till confluency and a wound was inflicted with a pipette tip. Subsequently two washes with 1XPBS were given to remove floating cells and media with additives was added. Additives to basal medium include FBS (5%), EGF (10 ng/mL), TGF-β (10 ng/mL), S1P (5 μmol/L), LPA (5 μmol/L), FBS (5%) + LPA (5 μmol/L). Images captured on Olympus IX71 microscope and analysed by TScratch Software (22). Significance of relative wound area covered for respective time points was determined by one way ANOVA using Bonferroni test. For three-dimensional growth analyses, cells were cultured in suspension (in ultra low adherent dishes for suspended culture in medium devoid of serum) or in Matrigel (23); clumps of 4 to 10 cells were obtained by dislodging cells through scraping and gentle pipetting. Cell clumps were observed under microscope for size and cell number counted. Larger clumps were fragmented by pipetting. At day 4 and 9, spheroids were fixed and immunofluorescence stained. Images were captured, using CellD image capture software.

Cell line models were immunostained for cell-specific expression markers MET (TCF21), EMT (Slug), ovarian stroma (PEG3) and secretory FTE (PAX8). Immunofluorescence staining and imaging was carried out as described earlier (15). In brief, cells grown on coverslips were fixed, blocked with 5% bovine serum albumin (BSA), and incubated with primary antibody for 1 hour. After two washes, cells were probed with labeled secondary antibody. Animal experimentation was done in accordance with the rules and regulations of the National Centre for Cell Science (NCCS) Institutional Animal Ethics Committee (IAEC). Xenografts were raised as described earlier (15). In brief, 2.5 × 10⁶ cells were injected in nonobese diabetic/severe combined immunodeficient mice subcutaneously in the thigh. Animals were maintained under pathogen-free conditions and tumor growth assessed every 2 days until the tumor diameter was approximately 1 cm, whereupon animals were sacrificed to harvest tumors.

Expression of TCF21, Slug, PEG3, PAX8, and ECM molecules THBS2, SPARC, FN1 was assessed in three xenograft models per cell line. In brief, sections were deparafinized, hydrated, then blocked with 5% BSA and incubated with primary antibody overnight. Next day, sections were incubated with secondary antibody, HRP-conjugate (Invitrogen), and color developed by DAB (Pierce). After counterstaining with hematoxylin (Sigma), sections were dehydrated, paraffinized, and mounted. Dilutions of primary and secondary antibodies used were as recommended.

Genomic instability was assessed by Immunofluorescence staining of (α-tubulin, H2AX pericentrin, Hoechst) in cell lines and xenograft models. Ten to 15 nuclei from 10 fields were randomly chosen and scanned for identifying different types of mitotic and cell division abnormalities and determining percent frequency of each aberration. Profiling of cellular Glutathione [GSH; determined by flow cytometry of monochlorobimane (Invitrogen) stained cells], H2AX, and Annexin V-FITC was carried as described earlier (24). In brief, 10⁶ cells were incubated with 25 μmol/L monochlorobimane for 30 minutes at 37°C, washed twice with 1× PBS, and intensities acquired using Amcyan filter on Canto II flow cytometer (BD). Since GSH acts as a primary antioxidant, its intracellular levels are often assayed as a measure of reactive species of oxygen (ROS) and nitrogen (RNS) in cells that lead to oxidative stress. Thus, we assayed GSH levels through competitive binding of Mono-ChloroBimane (mBCl) to glutathione-S-transferase (GST) that is involved in the oxidation of GSH by ROS and RNS. Lower levels of mBCl-GST indicate more ROS (as seen in OVCAR3 than in A4 cells), which in turn can induce DNA damage and increased H2AX foci formation. Profiling for DNA damage and apoptosis achieved through profiling H2AX and Annexin V-FITC as described. In brief, cells were fixed with 4% paraformaldehyde, permeabilized using 2% Triton X, blocked with 10% goat serum for 20 minutes, incubated with primary antibody for 30 minutes and secondary antibody for 20 minutes. H2AX was added to a final concentration of 1 μg per 1 million cells. A minimum of 10,000 cells were analyzed for each dataset and fluorescence intensities were measured on a logarithmic scale of fluorescence of four decades of log. Data were collected, stored, and analyzed with the FACStation software (BD Biosciences). BRCA1 and BRCA2 specific oligos were designed for amplification of genomic DNA; amplification products were processed for Sanger sequencing as described earlier (25). Primer sequences are available on request.

Tumor stratification was achieved using two approaches for gene expression pattern clustering, derivation of a core set of metastases-associated genes, and resolution of independent weighted correlation networks.

Expression of five EMT-TFs TCF4, TWIST1, SNAI2, ZEB1, and ZEB2 correlated positively with Vimentin (mesenchymal marker) and each other; and negatively with CDH1, KRT8, KRT18 (epithelial markers) in the serous ovarian adenocarcinoma data across two platforms (Affymetrix and Agilent) in The Cancer Genome Atlas Research Network (TCGA; Supplementary Table S1; Fig. 1A and Supplementary Fig. S1A; ref. 24). Twist, Snail, and Slug global targets were identified through genome-wide ChIP-on-chip. Significant enriched target probes and genes (P < 0.05) were identified for each transcription factor and correlated with their expression levels (extracted from GSE18054, the A4 gene expression dataset generated and deposited by us with GEO; refs. 26, 27) to identify repressed genes (strong targets) and genes with low level expression (that may possibly result from a dose dependent regulation and hence are considered as weak targets). Thus, 58 repressed and 135 low expressed Slug target genes; 91 repressed and 206 low expressed Snail target genes; and 107 repressed and 200 low expressed Twist target genes were identified (Supplementary Fig. S1D). Of these, only 33 repressed genes and 66 low expressed genes harbored an E-box (CANNTG, the consensus sequence recognized by the EMT-TFs), either within the probe or 500 bp up- and downstream of the probe. These 99 common target genes constitute the consolidated EMT targets in this study (Supplementary Table S2). Toward exploring the effects of LPA-S1P receptor–mediated signaling, an initial listing of the components was derived by manual literature annotation and from the REACTOME, KEGG, and Pathway Interaction Databases (http://www.reactome.org/ReactomeGWT/entrypoint.html, http://www.genome.jp/kegg/pathway.html, http://pid.nci.nih.gov/PID/index.shtml). This was followed by Pearson correlation analysis (0.4 < r > −0.4) across TCGA Affymetrix and Agilent platform generated data, to establish their relevance to serous ovarian adenocarcinoma. Thereby, 61 LPA-S1P associated and strong correlating genes were revealed (Fig. 1B and Supplementary Fig. S1B; Supplementary Table S2). Further hierarchial clustering, median absolute deviation (MAD) analyses and Pearson correlation between these 160 (99 EMT + 61 LPA-S1P) genes led to the derivation of a core set of 39 EMT-LPA-S1P signaling genes (Table 1; Supplementary Fig. S3).

Clustering identified class-specific gene-sample association and assigned a projection to each sample within individual classes (detailed in Supplementary Fig. S2). Class comparison between the silhouette- and consensus-based clustering suggested a high degree of overlap (Supplementary Fig. S1C) and further reaffirmed in hierarchical clustering. Unsupervised principle component snalyses and between group analyses lent support through similar class resolution and class-specific gene–sample associations (Fig. 1C and Supplementary Fig. S1C). The class-specific distribution of 286 samples classified was - Class1 (97), Class2 (119), Class3 (70; Fig. 1D).

Concurrently, we evaluated a metastases-independent classification scheme using Weighted Gene Correlation Network Analysis (WGCNA) that determines Pearson correlation coefficients across all possible pairs of genes by generating an adjacency matrix and organizing expression intensities to identify unique functional modules of highly correlating genes (28). Applying such a matrix at power 6 (weighted correlation = correlation⁶) to a core list of 6,840 significant genes generated 13 exclusive Blockwise Modules on the dendrogram using Dynamic Tree Cut algorithm (Fig. 2A). The detailed WGCNA methodology is provided as Supplementary Data S1; modules and genes listed in Supplementary Table S4; M0 represents genes unassigned to any module. Each module comprised of genes with high inter-connectedness and an expression pattern distinct from other module genes as revealed in the Topological Overlap Matrix plot (Supplementary Data S1 and Supplementary Fig. S1A).

K-means clustering in conjunction with Top 25 highest correlating genes from all modules led to the resolution of 3 tumor classes (Class1-127, Class2-137, Class3-95; Supplementary Data S1 and Supplementary Fig. S1B). A striking feature that was observed related to the similar distribution of specific tumor samples in the metastases- and WGCNA-based classifiers (Fig. 2B). This encouraged us to extract the core samples through an overlap of two approaches (Class1-77; Class2-99; Class3-51; Supplementary Data S1 And Supplementary Fig. S2) that were considered in the further analyses. A flow-chart of the derivation of these three classes and enriched samples is summarized in Fig. 2C.

Further derivation of a comparative median expression profile graph and heatmap of the common samples resolved class-specific functional module associations. Thus, enrichment of modules M1a, M6, M11 in Class1, that of M1b, M3, M9 in Class2, and M2, M4 in Class3 respectively was identified (Supplementary Data S2 and Supplementary Fig. S2; Fig. 2D). On profiling the distribution of 39 metastases-associated genes across the module genes, it was seen that LPA-S1P genes were represented in M1b and M3 (Class2), whereas EMT-TF targets were represented in M4 (Class3; Supplementary Data S1 and Supplementary Table S2). Toward further independent evaluation of the predicted class-module associations across platforms and data in different submissions, we assessed the same in 7 other publicly available datasets (Supplementary Data S2). Briefly, the Top 50 genes from enriched modules were used for clustering and class identification; these classes and module associations were further statistically affirmed through rigorous Module Preservation Analysis.

Class1 and Class2 module associations validated well and network based preservation statistics (Z summary and medianRank) of these across all datasets indicated a strong preservation of modules M1a, M1b, M3, M6, M9, M11. However, an inconsistent correlation of M2, M4 genes with Class3 tumors and poor module preservation statistics in the validation datasets led us to discontinue Class3 from further analyses. We also affirmed the gene exclusivity within enriched modules in Class1 and Class2 tumors through profiling of their interactor networks in ARACNe (29). Within a tumor class, each module displayed strongly coregulated nodes and their hubs of influence besides exhibiting intermodular connectivity through a few hub components (Fig. 3A and B), that could possibly identify cross-talk and regulatory features.

We further speculated the anticorrelative expression patterns of modules M1a-M6-M11 (Class1) and M1b-M3-M9 (Class2) tumors (Fig. 3C) to possibly instruct class-defining features. An initial exploration using GSEA (21) indeed revealed class-specific biologic functions that unfortunately were not too comprehensive (Supplementary Data S3). We thus resorted to detailed literature-based functional annotation of Top 50 genes in enriched modules to identify functional class associations (Supplementary Table S5).

In Class1 tumors, M1a genes are involved in altered chromatin remodeling; Wnt and HH signaling; metastases; regulation of quiescence, cell-cycle and differentiation; cell lineage commitment toward adipogenesis, chondrogenesis, myogenesis, neurogenesis, and spermatogenesis. M6 genes represent 3 functional categories, (i) Cancer Germline/Testes Antigens (CGA/CTA), whose association with meiosis and unexpected expression in tumors suggests abnormal chromosome segregation (aneuploidy), (ii) Cilia and Dynein related genes essential for ciliary architecture and mitotic spindle stability, (iii) Tumor suppressor genes. M11 genes are associated with integral ovarian functions, including steroid hormone biosynthesis, metabolism, stem cells, and epithelial differentiation.

Within Class2 modules, M1b genes define either hematopoietic cell surface molecules linked to inflammation, innate/adaptive immune mechanisms, or escape of tumor cells from immunosurveillance. M3 genes include extracellular matrix (ECM), EMT-associated, and dedifferentiation signature components known to be associated with TGF-β, PDGFRB, and Wnt signaling. Some M9 genes are components of the immunoproteasome and HLA Class I/Class II molecules involved in antigen processing and presentation while others define interferon networks (either being regulated by interferons or trigger their production) to mediate diverse biologic effects, including apoptosis, endothelial progenitor regulation, early immune stress response, posttranslational modifications (ubiquitination, DNA damage-dependent PAR), immune tolerance, and survival.

The above networks thus impart unique identities to the tumor groups in support of their being defined as molecular classes. Further validation of predicted biologic functions necessitates development of appropriate experimental models vis-à-vis class-specific cell lines that could serve to be convenient validation tools. On generating gene expression data-based silhouette plots of eight human serous ovarian adenocarcinoma cell lines (Supplementary Table S6) along with the TCGA tumors, it was observed that OAW42, OVCAR3, HEY, OV90, 2008 patterns clustered with Class1 while those of A4, CP70, SKOV3 clustered with Class2 tumors, albeit SKOV3 being represented as an outlier (Fig. 4A). The availability and an earlier familiarity of handling OAW42, OVCAR3 (Class1) and A4, CP70 (Class2) cells in our lab led to their selection in validation. The predicted expression patterns concurred fairly within the 4 cell lines tested (Fig. 4B).

Module M3 (enriched in Class2) functions suggest niche-epithelial interactions that link an altered microenvironment milieu with transcriptional regulation in tumor cells primarily mediated by the EMT-TF Slug. This presents a major differential modality of metastases between the two tumor classes, that was confirmed through the following assays:

• In vitro scratch assays—OVCAR3 and OAW42 (Class1) and A4 and CP70 (Class2) cells demonstrated differential invasive behavior vis-à-vis time required for wound healing, growth characteristics at the migratory edge, and growth factor influences. Both Class1 cells required around 120 to 144 hours in the absence of media supplements to heal at least 80% of the wound; A4 and CP70 cells were more efficient and achieved complete wound healing within 48 and 96 hours, respectively under similar conditions (Supplementary Figs. S4 and S5). The response of Class1 cells was enhanced by FBS + LPA, whereas that of Class2 cells was maximally triggered by EGF, TGF-β, FBS+LPA (alone or in combination). An important qualitative difference between the two classes was that wound healing in Class1 cells was effected exclusively through rapid cell proliferation that achieved extension of the scratch edge into the wound without dissolution of cell junctions; whereas A4 and CP70 cells at the scratch edge undergo EMT, actively migrate into and proliferate to swiftly fill up the wound. In support of these observations, lack of vimentin and strong expression of E-cadherin was seen in OVCAR3 cells, whereas strong vimentin with lowered E-cadherin expression was observed in A4 cells (Fig. 4C).

• Evaluating influence of ECM components on tumor behavior in non-adherent cell systems - Monitoring the growth of representative cells in suspension culture (no ECM), or ECM-enriched Matrigel-based sandwich cultures revealed novel effects of niche-tumor cell cross-talk. Single, suspended Class1 cells in either culture system failed to divide and underwent apoptosis within 10 to 12 days. However, when seeded as groups of 4 to 10 cells with intact cell junctions, they proliferate to generate compact, organized noninvasive spheroids in both systems, in which, a continual expression of E-cadherin but not vimentin is associated with these spheroids (Fig. 4D). In direct contrast, single Class2 cells survive, proliferate in both culture systems, and generate invasive mesenchymal cells as well as multicellular spheroids. The former settle/invade and form colonies at the bottom of the culture dish (Supplementary Fig. S6). Suspended Class2 spheroids in liquid culture are compactly organized, express E-cadherin but not vimentin—representative of MET (30). Spheroids in Matrigel range from small-compact, irregular-stellate, or large-organized shape and have invasive, vimentin-expressing cells at their borders (Fig. 4D, ii). All effects were more pronounced in A4 cells. The differential behavior between the two classes lend strong support to the purported correlation of Class1 being EMT-deficient and Class2 being EMT-enabled.

We proceeded to probe the differential effects of M1a genes (chromatin remodeling, cell cycle) and M6 genes (Cancer Testes Antigens, dyneins, tumor suppressors) that are preferentially enriched in Class1 over Class2, on altered cell cycle/division. This was carried out in OVCAR3 (Class1) and A4 (Class2) cells that were preferred due to a shorter latency periods of tumor generation in mouse models over OAW42 and CP70 xenografts. Immunofluorescence staining (α-tubulin, pericentrin, Hoechst) revealed more frequent anomalies of chromosomal and nuclear material rearrangements, including aberrant metaphase chromosome alignments, multipolarity, lagging anaphases, ring-shaped chromosomes, micronuclei, and multinucleation following cytokinesis defects in OVCAR3 than A4 cells in culture and in tumors (Fig. 5A, i and ii). While a higher frequency of metaphases in OVCAR3 tumors suggest them to be proliferative, increased nuclear fragmentation and apoptosis leads to delayed tumor latency and growth in these tumors over Class2 xenografts [Fig. 5A (i) and B]. Together, the data suggest more frequent cell-cycle defects and genetic instability in Class1 tumors.

A p53 signature, BRCA1/BRCA2 mutations, oxidative stress–driven DNA damage, genetic instability, and cell hyperactivation define the recently identified etiology of preneoplastic fallopian tube intraepithelial lesions leading to HGSC (6, 7) that challenged the existing dogma of the OSE being its cell-of-origin. Considering this conundrum in view of the differential biologic behavior of the two tumor classes, we speculated the etiology of Class1 and Class2 tumors to be the FTE and OSE respectively, and delved deeper into class-specific features in searching for proof-of-concept support as follows:

• p53-signature—Mutated p53 is widely associated with HGSC (24); OVCAR3 established from a sporadic case of HGSC has been earlier associated with a loss-of-function p53 mutation (31), whereas A4 harbors wild-type p53 (26). Accumulation of nuclear p53 is also evident in these tumors along with DNA damage (γ-H2AX) and a high proliferative index (Ki-67; Fig. 6A)

• BRCA1-BRCA2 mutational status—Novel point mutations in BRCA2 leading to nucleotide alterations 82800g>A and 82801a>G in exon 27 were identified in OVCAR3 cells that effect the non-synonymous amino acid change E3256R (Fig. 5C). This site lies in the RAD51 binding domain of the protein, which is critical for its functionality in DNA repair (32).

• Oxidative stress-DNA damage—We probed this purported etiologic linkage in both cell systems by profiling their intrinsic ROS and γ-H2AX expression (early DNA double-stranded break, DSB marker). Since glutathione (GSH) acts as a primary antioxidant, its intracellular levels are often assayed as a measure of ROS and RNS in cells that lead to oxidative stress. Thus, we assayed GSH levels through competitive binding of Mono-ChloroBimane (mBCl) to glutathione-S-transferase (GST) that is involved in the oxidation of GSH by ROS and RNS. Lower levels of mBCl-GST indicate more ROS (as seen in OVCAR3 than in A4 cells), that in turn can induce DNA damage and increased H2AX foci formation. As seen from the assays, both effects viz. higher ROS levels (lower levels of GSH) and DNA damage (γH2AX foci formation) were more pronounced in OVCAR3 than A4 cells (Fig. 5D, i–iii). Complementation of DSB with the earlier observed cell cycle and division defects also were more frequently observed in OVCAR3 than A4 cells, and suggest an intrinsic DSB propensity possibly triggered through an upstream release of ROS that together with the p53 and BRCA2 mutations ultimately culminates in genetic instability in these cells.

• Cell-specific expression patterns—We further profiled cell-specific expression markers including those of MET - EMT [TCF21 and Slug, respectively; ref. 33), ovarian stroma (PEG3; ref. 34), and secretory FTE (PAX8; ref. 35) in the representative models. A clear association of expression of these markers in OVCAR3 and A4 cells and experimental tumors with the putative etiologies of the two molecular classes was observed (Fig. 6B). We also probed the expression of ECM molecules purported to drive EMT (defined by the Class2-enriched modules M1-M3 association). As expected, Thrombospondin2 (THBS2), Osteonectin (SPARC), and Fibronectin (FN1) were expressed only in Class2 experimental tumors (Fig. 6B). Together, these data provide preliminary support to our hypothesis of differential cell-of-origin between the two classes, with Class1 involving FTE-ovarian stroma cross-talk and Class2 exhibiting features of OSE. However, further validation through lineage tracking and class-affirmation through correlations in patient tumors between all the risk factors from an epidemiologic viewpoint are essential toward a more complete validation of this predicted etiology.

Molecular classification of HGSC has been reported earlier; our findings are comparable with and in fact complemented by two recent studies.

Comparing actual tumor samples stratified in our classifier and those in the four HGSC TCGA subtypes (24) revealed Class1 tumors as being either proliferative or differentiated (76.62% and 12.99% respectively), Class2 tumors as mesenchymal or immunoreactive (53.53% & 31.31%, respectively), and Class3 tumors as differentiated or immunoreactive subtypes (54.9% and 25.49% respectively; Supplementary Table S7).

While our manuscript was under compilation, two HGSC classes were reported (36) based on their phenotype as iE and iM classes (epithelial and mesenchymal respectively) resolved through miRNA regulatory network analysis of the TCGA data. Some extent of phenotypic and sample correlation of this classifier with our classes was revealed (Supplementary Table S8). However, although Class1 exhibits epithelial features, our identification of an association of cancer testes antigens, stem cells-, differentiation-, development-, cell cycle-, chromatin modulation- and ovarian function-associated genes suggest unique complementation of cellular programs in this group of tumors that, in the context of metastases, are incapable of undergoing EMT-mediated invasion. Similarly, although our Class2 tumors in direct contrast are EMT-driven, they cannot be assigned a simple “mesenchymal” phenotype since our study also predicts significant cross-talk of the Class2-specific EMT program with extracellular matrix components and other gene networks associated with inflammation, innate/adaptive immune mechanisms, escape of tumor cells from immunosurveillance, antigen processing, and IFN responses. Together, all these features contribute to and are crucial determinants of their identity.

The relevance of EMT-based signatures as functional classifiers in identification of aggressive disease, and that of genetic instability and CTA expression in the choice of therapy are already established (37, 38). HGSC classification is currently being investigated intensely in view of developing personalized approaches (39). Our approach in the field led to the resolution of three molecular classes, two of which exhibit characteristics that could be important considerations in clinical management. Similarities exist between our classes with earlier relevant classifiers. Comparison with the TCGA subtypes suggests that our classification defines tumors in which the TCGA defined phenotypes may coexist in the same sample, yet assigning such tumors to a single phenotype may be an incomplete assessment of its biological functions. The iE-iM classification on the other hand is a limited description of the biological functions resolved in our stratification as indicated above. Since EMT is a tightly regulated program that is reversed during disease progression in the development of micrometastases (and involves re-establishment of an epithelial phenotype), Class2 tumors can never entirely be considered as being mesenchymal. Furthermore, the fact that a large majority of HGSC tumors present with mixed epithelial–mesenchymal features is overlooked in the iE-iM stratification but are likely to belong to our Class3 and “non-core” samples. Further resolution of such heterogeneity in unclassified samples may be achieved through integration with other functional signatures or profiling differential mechanisms of gene regulation (mutational profiles, methylation, and proteomics) that are increasingly applied in addressing tumor heterogeneity.

On another front, a commonality between some class-specific features and predisposing events of FTE or OSE transformation in HGSC etiology was realized. FTE transformation involves—(i) post-ovulatory oxidative stress, (ii) DNA damage followed by defective repair due to BRCA1/2 mutations in patients that lead to neoplastic changes, genetic instability, and malignant transformation, (iii) precursor lesions characterized by a p53 signature (p53 mutations, high nuclear p53 protein expression and proliferative index), and (iv) expression of the differentiated Müllerian epithelial marker (PAX8) in inclusion cysts within the ovarian stroma. In exploring such features in our experimental models, we did identify high oxidative stress, DNA damage, BRCA2 mutations in its RAD51 binding domain, p53 signature and genetic instability along with PAX8 expression in OVCAR3 cells (Class1). Further association of TCF21 and PEG3 and ovarian functions (steroid hormone biosynthesis—FOXL2, STAR, AMRH2) and germ cell differentiation (CTA) suggests cross-talk between FTE and ovarian stroma, which is supported by earlier reports of both fallopian tube or ovarian stroma-derived mesenchymal stem cells being capable of MET (40, 41) and niche interactions between preneoplastic lesion–derived FTE and the ovarian cortex being crucial in HGSC (42). Collectively considered, these suggest that Class1 tumors may originate from the FTE (Fig. 6C).

Events associated with OSE transformation include (i) OSE rupture during ovulation accompanied by EGF and TGFβ production triggers a signaling cascade activating Slug/Snail/ZEB1 to initiate EMT-mediated wound healing (42). Incessant ovulation leading to aberrant EMT is thus implicated in OSE transformation (7), (ii) OSE secretes several ECM components, including collagens, integrins, laminins, fibronectin, and vitronectin (43); some of which are implicated in epithelial cell transformation, development of a tumor niche and EMT. These include SPARC (enhanced Slug expression; 44); Decorin (downregulates E-cadherin; 45), Fibronectin (enhanced Slug expression; 46), Periostin (activates Snail via β-catenin; 47), and activation of TGF-β (48). The altered ECM milieu thereby triggers an EMT program that generates mesenchymal cells with a capacity of invasion and metastases; such tumor-derived mesenchymal cells contribute to a “mesenchymal” signature. Module M3 enriched in Class2 tumors of our study reflects on such ECM-EMT cross-talk and represents a different modality of metastases from that in the Class1 tumors that are EMT-deficient. (iii) OSE also triggers cytokine production (IL-1, IL-6, G-CSF, M-CSF, GM-CSF) toward follicular growth, differentiation, and ovulation; altered cytokine profiles are associated with transformation (49, 50). Effectively, aberrant derivations of OSE-functions appear to be recapitulated in Class2 tumors as a Slug-ECM-stromal signature, inflammation, IFN signaling and altered host immune responses (Modules M1b, M3, M9; Supplementary Table S5). While their stringent regulation maintain normal organ functioning and homeostasis, these network perturbations in Class2 tumors reflect on intrinsic OSE properties gone awry culminating in transformation (Fig. 6C).

In conclusion, our resolution of HGSC tumor classes generates a much wider and in depth class-specific understanding of the transformation-progression process than that achieved earlier, while their associated networks reveal a novel possibility regarding cell of origin for each class. Further affirmation of the latter through lineage-tracking approaches in elegant disease models are now necessitated to provide affirmation and final validation. Irrespective of this finding, our HGSC stratification has important ramifications toward improving HGSC outcomes through achieving a long-term goal of patient stratification for targeted therapy.

No potential conflicts of interest were disclosed.

Publicly available TCGA and GEO datasets were used in this study and are referenced.

Development of methodology: N.L. Gardi, T.U. Deshpande, S.A. Bapat

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.C. Kamble, S.R. Budhe

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): N.L. Gardi, T.U. Deshpande, S.C. Kamble, S.R. Budhe, S.A. Bapat

Writing, review, and/or revision of the manuscript: N.L. Gardi, T.U. Deshpande, S.C. Kamble, S.A. Bapat

The authors also acknowledge excellent technical assistance from Mr. Avinash Mali.

This work was supported by the Department of Biotechnology, Government of India, New Delhi (BT/PR11465/Md/30/145/2008, BT/IN/FINNISH/25/SB/2009; to S.A. Bapat). S.C. Kamble receives a research fellowship from the Council of Scientific and Industrial Research, New Delhi, India.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.