Distinct molecular pathways in ovarian endometrioid adenocarcinoma with concurrent endometriosis

Since OEA is a rare histotype of epithelial ovarian cancer 11 , 12 and a majority of women with OEA do not have endometriosis at time of staging 13 – 16 , we used multiple tissue banks to obtain adequate number of samples for analysis. Supplemental table 2 details the metadata for each sample including clinical demographics, tissue source, and RNA quality control metrics. Table 1 summarizes clinical and pathological characteristics. Women without OEA were significantly ( P =0.0001) younger [median: 38 years, range (25–48)] than women with OEA [median 52 years, range (37–76)]. In women with OEA, women without concurrent endometriosis had a trend towards higher stage disease ( P =0.051) but did not have a significant difference in grade. Thus, the samples used for molecular studies are clinically similar. We performed genome-wide transcriptome profiling on RNA isolated from specimens ( n =20) collected from endometriomas ( n =4), endometriomas with epithelial atypia ( n =3), proliferative endometrium ( n =4), OEA without endometriosis ( n =4), and OEA with concurrent endometriosis ( n =5). Considering the mix of pathologies ( e.g ., benign and malignant), tissue types ( e.g ., endometrium, cancer, endometrioma cyst wall), number of experimental batches ( e.g ., only six samples on a single microarray), multiple sources of specimens ( e.g ., tissue bank), we could not adequately conduct batch effect removal. Small total RNA abundance may cause a bias in microarray data, in particular Illumina microarrays 27 , 28 . Our analysis of the whole dataset suggests the samples with too small amount of total RNA are more like outliers compared to the other samples of the sample disease class. Hence, the samples with total RNA amount less than 30 µg at time of RNA isolation from tissue were excluded from the final microarray analysis ( Supplemental Figure 1 ). To validate this exclusion analysis, we compared our datasets to published datasets. We compared our benign datasets ( e.g ., endometrioma and proliferative endometrium) with four publicly available and previously published data sets 17 , 23 – 25 . We found that the differentially expressed genes in endometriomas versus endometrium ( Supplemental Tables 3 and 4 ) were highly enriched among differentially expressed genes identified by other published datasets (pooled P -value <1E-30) ( Supplemental Table 5 ). Only one previously published dataset exists for OEA with concurrent endometriosis, but the data are not publically available 29 . Therefore, our datasets for OEA with and without concurrent endometriosis represent unique publically available research datasets. Next, gene expression profiles ( n =15) were evaluated using principal component analysis ( Figure 1 ). Graphical representation of principal components (PC) 1 and 2 shows that gene expression profiles from malignant samples clustered separately from benign samples ( Figure 1A ). Graphical representation of PC1 and PC3 shows that OEA samples with endometriosis cluster separately from those without endometriosis ( Figure 1B ). Thus, at a global level, OEA with endometriosis is molecularly distinct from OEA without endometriosis. To identify the biological characteristics of each principal component, we conducted pathway enrichment analysis on the 200 genes that are most positively/negatively associated with the loading of each principal component ( Supplemental Table 6 ). Glycolysis ( P =0.04), pyruvate ( P =0.019), cysteine ( P =0.013), and alanine/glutamate ( P =6E-5) metabolism; P53 ( P =0.04), WNT ( P =0.04), and Ephrin receptor (axon guidance, P =0.004) signaling; and RNA polymerase ( P <0.01) and DNA repair ( P =0.002) pathways are negatively correlated with PC1. Extracellular region and cell adhesion related pathways including cell adhesion ( P =2E-12), cytokine-cytokine receptor ( P =0.01), extracellular matrix (ECM) receptor interaction ( P =5E-10), vascular endothelial growth factor (VEGF, P =0.001), hypoxia-inducible factor 1-alpha (HIF1A, P =0.04), and glycoprotein metabolic ( P =8E-13) pathways are positively correlated with PC1. Cell cycle ( P =4e-5), DNA replication ( P =3E-5), and glycolysis ( P =0.04) pathways are negatively correlated with PC2. Apoptosis ( P =0.001), cytokine-cytokine receptor ( P =0.01), coagulation ( P =7E-08), extrinsic inflammation ( P =0.001), cell-matrix adhesion ( P =7E-05), and chemokine signaling ( P =0.01) pathways are positively associated with PC2. Transforming growth factor beta (TGFβ, P =0.01) and cell cycle ( P =0.01) pathways are negatively associated with PC3. Immune response ( P =2E-4), cytokine signaling in immune system ( P =0.002), T cell signaling ( P =0.02), and interleukin (IL1, IL6) signaling pathways ( P <0.01) are positively associated with the PC3 ( Supplemental Table 6 ). Such pathway associations suggest that OEA with and without endometriosis are most similar in terms of glycolysis and P53 and WNT signaling (PC1) and cell cycle (PC2), while proliferative endometrium and endometrioma are most similar in terms of pathways involved in dysregulated cell adhesion, extracellular matrix (PC1), and an inflammatory response (PC2, Figure 2A ). Figure 2B shows a distinct difference between OEA with and without endometriosis (PC1 and PC3). This distinction is mediated by PC3 with significant contributions from TGFβ, cell cycle, and immune response ( Supplemental Table 6 ). Our pathway analysis of the PC-associated genes suggests that OEA with endometriosis has increased TGFβ signaling compared to OEA without endometriosis. OEA without endometriosis, endometrioma, and proliferative endometrium have similar level of dysregulated immune and inflammatory response. A close distance between atypical endometrioma and OEA without endometriosis relies heavily on this similarity of dysregulated immune and inflammatory response. Detailed lists of the genes and their enriched pathways are given in Supplementary Table 6 . We directly compared gene expression profiles of OEA samples without endometriosis ( n =4) to OEA with concurrent endometriosis ( n =3) to determine molecular distinctions between these histologically similar tumors (log-fold change > 0.5 or < −0.5, P <0.01; Supplemental Tables 7 and 8 ). We discovered 1022 probes differentially expressed, with 239 probes downregulated and 783 probes upregulated. Of these, 526 mapped IDs corresponding to 497 unique genes were upregulated and 206 mapped IDs corresponding to 184 unique genes were downregulated in OEA with concurrent endometriosis compared to OEA without endometriosis. Figure 2 shows the volcano plot of differentially gene expression analysis and heat map of differentially expressed genes. To explore potentially impactful pathways and processes related to the differentially expressed genes, we conducted pathway enrichment analysis by using IPA, DAVID 20 , 21 , and gene set enrichment analysis 22 . From IPA, the most statistically significant canonical pathways with a positive Z-score were leukocyte extravasation signaling ( P =5.13E-05) and acute phase response signaling ( P =4.07E-04), two pathways involved in inflammation. The most statistically significant canonical pathways with negative Z-scores were gonadotropin releasing hormone (GNRH) signaling ( P =1.20E-03) and bone morphogenetic protein (BMP) signaling ( P =6.17E-03). For these four statistically significant canonical signaling pathways, Table 2 lists the differentially expressed genes and fold-change. Figure 3 depicts a waterfall plot of the statistically significant canonical pathways ( P <0.05) dysregulated in OEA with concurrent endometriosis. Out of those 35 canonical pathways, 25 have contributions from nuclear factor kappa-light-chain-enhancer of activated B cells (NFkB) signaling ( Supplemental Table 9 ). Overlaying this gene list on the molecular mechanisms of cancer canonical pathway shows activation of RAS signaling ( Supplemental Figure 2 ). DAVID functional clustering analysis identified 153 and 92 functional clusters enriched by the up- and downregulated genes, respectively. From DAVID in support of RAS signaling, the top upregulated functional clusters included RAS association, protein kinase, response to stimulus, immune response, and regulation of cell development. The top downregulated functional clusters included endoplasmic reticulum, response to hormone stimulus, anti-apoptosis, and cell migration and adhesion. A complete list of the identified functional clusters is given in Supplementary Tables 10 and 11 . Pathway enrichment test against MsigDB canonical pathways identified 117 pathways enriched by the upregulated genes, including interleukin signaling ( P =2E-04), RAS homolog family member A (RHOA) regulation ( P =0.001), C-X-C motif chemokine receptor 4 (CXCR4) signaling ( P =0.003), and kinase activity ( P =2E-04), pathways that are related to immune and inflammatory responses. We also identified 325 pathways enriched by the downregulated genes, including integrin ( P =1.21E-06), cell-cell adhesion ( P =1.6E-05), cell proliferation ( P =1.4E-04), NOTCH ( P =2.3E-04), WNT ( P =3.8E-04), and hypoxia inducible factor 1 alpha subunit (HIF1A) ( P =4E-04) signaling pathways. A complete list of the MsigDB enriched pathways is provided in Supplementary Tables 12 and 13 . We validated the expression of selected individual genes using real-time QPCR with a focus on genes in unique molecular pathways in OEA or genes not studied in ovarian cancer ( Figure 4 ). Most of the genes followed a similar trend in expression via QPCR compared to microarray analysis. Plasminogen activator, urokinase ( PLAU ) was 3.8-fold upregulated on microarray analysis and followed similar direction of regulation in QPCR validation studies (1.6-fold change, P =0.008). There was no change in gene expression of the PLAU receptor ( PLAUR ), although studies suggest that PLAUR serum protein levels may play diagnostic roles for women with adnexal masses 30 . Gene members involved in BMP signaling ( Table 2 ), including CREB binding protein ( CREBBP ), protein kinase cAMP-activate catalytic subunit beta ( PRKACB ), and mitogen activated protein kinase 13 ( MAPK13 ), showed a statistically significant difference in gene expression compared to samples without endometriosis. A similar trend in appropriate direction was observed for SMAD family member 6 ( SMAD6 ) and nuclear factor kappa B subunit 2 ( NFKB2 ). Bone morphogenetic protein 8B ( BMP8B ) showed a 4-fold downregulation that was statistically significant ( P =0.002) but in opposite direction of microarray (2.3-fold upregulation). In silico inhibition of BMP leads to no predicted downstream signaling differences ( Supplemental Figure 3 ). Although catenin beta 1 ( CTNNB1 ) was not a statistically significant change, it was downregulated in QPCR studies similar in direction to microarray results. P21 (RAC1) activated kinase 2 ( PAK2 ) has shown to decrease migration of ovarian cancer cells in vitro 31 . PAK2 showed a 5.3-fold downregulation in microarray analysis and 2.2-fold change downregulation by QPCR ( P =0.0009). The gene expression changes were validated by independent QPCR.

After expedited IRB approval, de-identified flash-frozen specimens were obtained from the Human Tissue Acquisition and Pathology Core of the Dan L. Duncan Comprehensive Cancer Center, the Gynecologic Tissue Biorepository for the Department of Obstetrics and Gynecology at Baylor College of Medicine, and the Gynecologic Oncology Group (GOG). Histopathologic and demographic data were abstracted from de-identified surgical pathology reports for each study subject. Benign samples were obtained as previously described 17 . Total RNA was extracted using the mir Vana kit (Applied Biosystems, Foster City, CA). Only RNA samples that passed strict quality control using an Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA) underwent whole genome expression profiling using Illumina’s Human WG-6 version 3.0 BeadChip, which contains 48,804 probes covering 27,455 genes. This was performed as a fee for service at the Texas Children’s Cancer Genomics and Proteomics Core Lab. Raw CEL microarray data was processed and normalized by Bioconductor R package “beadarray” 18 under default parameter setting. Considering the sample size of the experiment, differential gene expression analysis was conducted by using empirical Bayes-based tests in the Bioconductor R package “siggenes” 19 . Differential gene expression was determined by a significant cutoff with P -value <0.001 and log-fold change 0.5. To fully identify the biological pathways and processes related to the differentially expressed genes, we conducted pathway enrichment analysis by using (1) Ingenuity Pathway Analysis (IPA), (2) functional clusters in DAVID 20 , 21 , and (3) a hypergeometric test-based approach against gene sets from MsigDB (molecular signatures database) 22 , with a significance cutoff P <0.05. Multidimensional scaling (MDS) plot was made by the top two principle components derived from expression profile of all probes. Normalized gene expression profiles of data sets GSE7305 , GSE5108 , GSE11691 , and GE23339 were retrieved from GEO database 17 , 23 – 25 . Due to the relatively large number of samples in each database and the fact that each of these four data sets were generated using different microarray platforms, we used the non-parametric Mann Whitney test for differential gene expression analysis, with P <0.01 as the cutoff for significance. Significance of overlap of the differentially expressed genes identified in different data sets was tested by Fisher’s exact test. RNA was treated with Turbo DNase (Life Technologies, Foster City, CA). DNase-treated RNA (1000 ng) was reverse-transcribed in a 50-µL-reaction volume with Superscript III (Life Technologies) and random primers (Life Technologies). Samples were diluted to 100 µL and two µL was used for each quantitative real-time PCR (QPCR) reaction. QPCR was performed using a QuantStudio 3 real-time PCR system, inventoried TaqMan gene expression assays ( Supplemental Table 1 ), and TaqMan universal PCR master mix II (Life Technologies) in 10-µL-reaction volume with reaction conditions as published 17 . Each sample was analyzed in duplicate and a no template control sample was included on each plate for each primer-probe set. Expression of 18S RNA was used as an endogenous control. The relative quantity (RQ) of individual transcripts was calculated using the 2 –ΔΔCT method 26 , plotted as mean ± SEM, with a Student’s t test used to generate P values for statistical significance.

Ovarian cancer is the seventh most common cancer in women, and the eighth leading cause of cancer-related death worldwide, claiming over 151,000 lives in 2012 32 . Currently, standard-of-care therapy treats nearly all women with ovarian cancer similarly, even though studies have shown that individual histotypes of ovarian cancer ( i.e ., high-grade serous, endometrioid, clear-cell, and mucinous) are molecularly distinct and likely arise from a unique cell of origin 1 . Thus, there is a clinical need to understand the key molecular drivers of cancer subsets such as those arising from endometriosis, which may drive treatment recommendations. Currently, whole genome sequencing analysis of ovarian cancer samples allows stratification into histotypes, including sub-classification of OEA into tissues that are microsatellite stable and those with microsatellite instability 33 . While large amounts of genomic data have been generated based on histotype, there is limited information regarding the molecular contribution of the endometriotic tumor microenvironment on OEA. While this subtle pathologic diagnosis is often overlooked in research studies, the improvement of outcomes of women with endometriosis at time of ovarian cancer staging is not 7 – 9 . Therefore, the molecular contributions of the endometriotic tumor microenvironment may affect the overall tumor biology. We have performed a comprehensive assessment of gene expression signatures in OEA from women with concurrent endometriosis. Our data suggests that OEA from women with endometriosis has a distinct molecular signature when compared to OEA from women without endometriosis. We have identified signaling pathways that appear to uniquely contribute to the pathogenesis of this subset of ovarian cancers with inflammatory, NFkB, RAS, and TGFβ signaling pathways playing a significant role. In support of our data, KRAS oncogenic mutations have been discovered in 29% of ovarian cancers with concurrent endometriosis 34 . Transgenic mice expressing oncogenic Kras develop both endometriosis-like lesions and OEA when a conditional Pten deletion is introduced in the ovarian bursa 35 . However, activated signaling cascades through RAS have not been specifically described until this report. The TGFβ superfamily has been implicated in ovarian cancers, particularly sex cord-stromal tumors 36 – 40 , but not specifically in epithelial ovarian cancers. Thus, these results support further study into the molecular contributions of the endometriotic tumor microenvironment. To date, only one study has used whole transcriptome microarray analysis to evaluate specimens of normal ovary, endometriomas, and OEA with and without concurrent endometriosis 29 . Results of this work revealed only a small group of cytokines dysregulated in OEA with concurrent endometriosis and endometriomas, consistent with the known inflammatory milieu of endometriosis. While those results suggest that the endometriotic tumor microenvironment may contribute to OEA, we were unable to compare our results to those datasets directly, as they are not publically available. Nonetheless, our current observations build on this earlier work by expanding and identifying additional gene products and signaling pathways that may contribute to the distinct clinical behaviors associated with subsets of patients diagnosed with OEA with concurrent endometriosis. A strength of the current study is that the specimens of OEA used in our study are histologically similar, containing >80% tumor cells, without mixed histology, and all OEA with concurrent endometriosis are confirmed by histopathology reports. Similar to other studies 7 , OEA without endometriosis represented a trend towards higher stage disease ( P =0.051). While there was no significant difference in high-grade disease, our studies were not powered for this analysis ( Table 1 ). However, due to our strict inclusion criteria, only a small number of tissue samples were available, reflecting the rarity of OEA and OEA associated with endometriosis 11 – 16 . Although OEA with and without concurrent endometriosis may arise from endometriosis 1 , recent studies show that OEA samples may be stratified into distinct biological features within the OEA histotype based on the somatic genome 33 . Our study is limited by the use of relative insensitive microarray technology compared to other technologies now available for whole exome profiling. For example, a key player in OEA, AT-rich interaction domain 1A ( ARID1A ) 14 , 41 is not represented on this microarray platform. The validity of our observations is supported by the fact that we were able to validate them by QPCR in a second, independent cohort of well-characterized tissue samples. Nonetheless, future studies with larger sample sizes on next-generation sequencing platforms may highlight additional signaling pathways and help to provide an even more complete picture of OEA and the molecular contributions of endometriosis. The molecular mechanism of transformation of endometriosis or atypical endometriomas into malignant ovarian cancers is being actively studied. While mutations in the oncogene, KRAS , and ARID1A have been frequently documented in endometriosis-associated ovarian cancers 42 , these mutations are frequently discovered in deeply infiltrating implants of endometriosis that do not typically progress to ovarian cancer 43 . Thus, the correlation between genotype and clinical phenotype still needs to be determined 42 . Clinically, it is thought that women with long-term untreated endometriosis are at highest risk of developing ovarian cancer 42 . For example, the protection of combined oral contraceptive therapy on ovarian cancer risk is more robust for women with endometriosis [odds ratio 0.21 (0.08–0.58, P =0.003) compared to non-endometriosis population 0.47 (0.37–0.61, P <0.001)] 44 . Both endometriosis and ovarian cancer are hormone responsive diseases 45 , 46 and thus while the overall steroid hormone suppression offered by combined oral contraceptives makes logistical sense, the actual molecular mechanism has not been elucidated and likely involves molecular, genetic, and hormonal factors. An interesting target within the context of the endometriotic tumor microenvironment is the contribution of NFkB to OEA. A majority of our canonical signaling pathways had significant contributions of NFkB ( Supplemental Table 9 ). Studies have shown that endometriosis progression relies on constitutive active NFkB, and NFkB signaling pathways are potential targets for non-hormonal therapies for endometriosis 47 , 48 . NOTCH, WNT, and HIF1A signaling pathways were also significantly enriched. Whether these signaling cascades are downstream of NFkB signaling cascades or if they arise independent of NFkB signaling from the cell’s response to the endometriotic tumor microenvironment are unknown. Deeper understanding on the role of NFkB signaling cascades, inflammation, and even microenvironmental stress in OEA with concurrent endometriosis may drive novel treatment recommendations with future studies. As an additional area under study, the origin of ovarian cancer is still relatively controversial. While some high-grade serous ovarian cancers are thought to arise from the malignant transformation of fimbria of the fallopian tube with subsequent early metastasis to the ovary 1 , 49 , others may arise from the ovary itself 1 , 50 . OEA is thought to arise from secretory epithelial cells that are frequently found in eutopic endometrium and ectopic endometriomas 2 . Further, the cell of origin may predict aggressive disease. For example, high-grade OEA without endometriosis may develop from the secretory epithelial cells of the eutopic endometrium 2 and already have a metastatic phenotype by the time it is discovered 42 . On the other hand, OEA with concurrent endometriosis may develop from secretory epithelial of endometriomas and represent a non-metastatic phenotype 42 . While our results examine the gene expression differences between OEA with and without endometriosis, these differences may be from the unique endometriotic tumor microenvironment or from the unique cell of origin. Due to low frequency of these clinical samples, additional model systems may be necessary to answer these questions. Our studies delineate key molecular pathways in OEA with concurrent endometriosis. Future studies should be undertaken to detail the role of inflammation including the contribution of constitutive active NFkB in endometriosis and RAS and BMP signaling in the clinical features of OEA. We believe that future studies targeting key signaling pathways may have implications for novel treatment. Thus, understanding the molecular features of OEA with concurrent endometriosis may have an impact on a significant number of women’s lives.

Most ovarian cancers arise from cells that are not normally found in the ovary 1 . For example, ovarian endometrioid adenocarcinoma (OEA) is thought to arise from secretory epithelial cells 2 of the eutopic endometrium or endometriosis. Endometriosis is a hormonally responsive, pathologic growth of endometrium outside the uterus ( i.e ., ectopic location), frequently discovered as a benign cyst on the ovary called an endometrioma 3 . The presence of endometriosis, itself a benign disease, increases the risk of ovarian endometrioid and clear-cell adenocarcinoma up to 8.9 fold, depending on genetic admixture and environmental exposures 4 – 6 . Clinically, studies suggest that co-occurrence of endometriosis with ovarian cancer is associated with an improved prognosis 7 – 9 . Although Dr. John Sampson first hypothesized the malignant potential of endometriosis in 1925 10 , the relationship between endometriosis and ovarian cancer is still being deciphered. Previous reports have identified a limited number of oncogenes or tumor suppressors that may be involved in the pathogenesis of OEA 3 . Mutations alone cannot explain the clinical and phenotypic differences, and thus, the endometriotic tumor microenvironment may play a significant role in the pathogenesis of OEA 2 . However, specific molecular contributions of the endometriotic tumor microenvironment to OEA remain poorly understood. At least in part, insight into the molecular contributions of endometriosis to OEA remains limited due to the rarity of well-characterized tissue samples. OEA accounts for fewer than 10% of all epithelial ovarian cancers 11 , 12 , and less than 43% of women with OEA have endometriosis at time of staging 13 – 16 . Consequently, few existing studies have focused on molecular features unique to OEA with concurrent endometriosis. To study the unique molecular contributions of endometriosis in OEA, we have directly compared the transcriptome of OEA with and without concurrent endometriosis with the goal of identifying previously unappreciated aspects of this disease. Our results provide insight into the critical gene networks important in OEA in the context of the endometriotic tumor microenvironment.

Supplemental Figure 1: Flow diagram of samples in microarray analysis. Strikethrough indicates samples that were removed from final analysis, based on total RNA<30µg. Supplemental Figure 3: Predicted signaling effects of BMP8B on BMP signaling pathway. A, BMP signaling pathway with high expression of BMP8B as indicated by microarray analysis. B, in silico predicted signaling effects with inhibition of BMP8B . Note little difference in downstream signaling. Orange shading, predicted activation; blue shading, predicted inhibition. Supplemental Table 1: TaqMan assay IDs for assays used in QPCR Supplemental Table 2: Details of clinical samples used in study Supplemental Table 3: Endometrioma down genes Supplemental Table 4: Endometrioma up genes Supplemental Table 5: Comparison of endometrioma differentially expressed genes to published studies Supplemental Table 6: Multiple tabs containing contributions of each principle component Supplemental Table 7: OEA with concurrent endometriosis down genes Supplemental Table 8: OEA with concurrent endometriosis up genes Supplemental Table 9: Listing of IPA canonical pathways unique to OEA with concurrent endometriosis, highlighting NFkB pathways. Supplemental Tables 10 and 11: List of DAVID functional clusters enriched by the differentially expressed genes in OEA with concurrent endometriosis versus. Supplemental Tables 12 and 13: List of MsigDB canonical pathways enriched by the differentially expressed genes in OEA with concurrent endometriosis. Supplemental figure 2: Molecular mechanism of cancer canonical pathway dysregulated in OEA with endometriosis. Overlay shows predicted signaling activation (orange) and inhibition (blue). Note predicted activation of RAS (left side) through inhibition of RAS-GTPase accelerating protein (GAP), leading to predicted activation of MAPK signaling including cyclin D1 (CCND1) activation.