Transcriptomic Characterization of Endometrioid, Clear Cell, and High-Grade Serous Epithelial Ovarian Carcinoma

Epithelial ovarian cancer (EOC) makes up approximately 90% of all ovarian tumors and is the fifth leading cause of cancer-related death among women in the United States (1); five-year overall survival remains around 45% for all stages (27% for distant disease). One reason for the lack of success with the treatment of EOC can be attributed to disease heterogeneity from both the morphological and molecular perspectives (2, 3). The most common histological subtype (histotype) of EOC is high-grade serous carcinoma (HGSC), representing 70% to 74% of EOC cases followed by endometrioid carcinoma (EC; 7%–24%) and clear cell carcinoma (CC; 10%–26%; ref. 4). HGSC has been extensively studied in terms of its molecular profile and genomic landscape (3, 5–9), in part due to its higher frequency. Some features of HGSC tumors are the following: believed to originate from fallopian tube secretory epithelium and/or ovarian surface epithelium cells; majority have somatic TP53 mutations; and several germline risk alleles associate almost exclusively with serous EOC (10, 11). In contrast, features of ECs and CCs include: appear to arise from endometrial epithelium (12); often harbor ARID1A mutations (13); and appear to have some unique risk alleles (10, 14, 15). Because of the rareness of these two EOC histotypes, limited research has been conducted on EC and CCs tumors.

We sought to characterize tumor transcriptomes of EC and CC in relation to HGSC, and to evaluate clinical associations using RNA-Seq in a clinically annotated EOC patient population, followed by replication of findings in publicly available expression array datasets. The determination of the molecular features that differ between these histotypes and how these transcriptomic differences relate to clinical response might lead to new insights in histotype-specific treatments for EOC patients.

Eligible patients were women age 20 years or above ascertained between 2000 and 2014 at the Mayo Clinic within one year of diagnosis with pathologically confirmed primary EOC. Clinical diagnoses and histotypes were confirmed by immunohistochemistry-guided re-review by a gynecologic pathologist. The pathologist also verified tumor grade and reviewed each fresh-frozen tissue to ensure 70% tumor content before RNA extraction. All cases provided informed written consent to protocols approved by the Mayo Clinic Institutional Review Board. Supplementary Table S1 summarizes the characteristics of the 186 EOC patients included in this study. Transcriptomic sequencing of RNA was performed in four batches using Illumina TruSeq library preparation (Stranded Total RNA Library Prep Kit or RNA Library Preparation Kit v2) and sequenced on the Illumina HiSeq 2000 sequencer with paired end reads (50- or 100-nucleotide, Supplementary Table S2). No technical replicates were included in the study. Primary analysis and de-multiplexing was performed using Illumina's CASAVA software, followed by alignment to the GRCh37 genome assembly using TopHat2 (16) and abundance estimation at the gene level using RSEM (17).

The normalization across RNA-Seq batch effects consisted of two steps: estimation of surrogate variables and normalization. The first step involved surrogate variable analysis (SVA) conducted on the residual matrix generated from a linear model fit of the gene expression abundance against histotype. We adopted the methodology proposed by Buja and Eyuboglu to identify surrogate variables that represent known and unknown confounders other than the primary factor of histotype (18). The algorithm to carry out the analysis was provided in the Bioconductor package SVA via function svaseq, with two approaches to determine the number of significant confounders or surrogate variables (18, 19). For the second step involving across sample normalization, we employed the ComBat function in Bioconductor to adjust for the surrogate variables and/or confounders other than histotype, using an empirical Bayesian approach proposed by Johnson and colleagues (20). Finally, we adjusted for the known factor (i.e., batch) and visually inspected the data with principal component analysis (PCA) to determine whether the batch effects were adequately removed (see Supplementary Methods).

Following normalization of the data to account for technical artifacts, we determined whether any differences in tumor purity remained between the samples and associated with histotype. Estimates of tumor purity were calculated using the program ESTIMATE (21), where we observe the association of tumor purity estimates with histotype to be nonsignificant (P > 0.10). SVA and ComBat normalization effectively removed the differences in library sizes between samples and difference in tumor purity between histotypes, thus eliminating the need for additional normalization of these factors (see Supplementary Methods).

Statistical analyses used the normalized count data, treating it as a continuous measurement, where both positive and negative values are possible. To assess similarity of global tumor transcriptomes, PCA was completed, along with unsupervised clustering using a parametric model-based approach implemented in mclust (22) using the top 75% most variable genes. Determination of number of optimal clusters was completed using the Bayesian Information Criteria (BIC). After assessing model fit for a random sample of 10 genes using the normalized gene expression data and three modeling methods [i.e., negative binomial model using edgeR (23), ANOVA, and nonparametric Kruskal-Wallis (KW) tests], we choose to use the conservative nonparametric KW tests to assess per-gene differential expression between all histotypes (2 degree of freedom tests; see Supplementary Methods). In addition to testing whether mean gene expression abundances differed between the three histotypes, we also completed three sets of pairwise differential expression analyses (1 degree of freedom tests) using non-parametric two-sample t tests (Wilcoxon tests).

Association of gene expression with outcome was completed with progression-free survival (PFS) defined as time from the date of diagnosis to the date that second-line therapy was initiated for a clinically actionable tumor recurrence or death, if progression were unknown. Cox proportional hazards models were used to estimate hazard ratios (HRs) and 95% confidence intervals (CIs) for association of each gene expression with PFS. Analyses adjusted for covariates associated with PFS (P < 0.05 based on score test), including stage (coded as I/II and III/IV) and surgical debulking status (coded as optimal, sub-optimal). Censoring at 10 years was performed to minimize competing causes of mortality. We excluded four patients with missing or unknown debulking status from the PFS analysis. We applied a less stringent statistical significance threshold for PFS analysis compared with the differential gene expression analysis due to the smaller effect sizes one expects to see for gene expression association with clinical outcome. To determine if a different relationship exists between gene expression and PFS by histotype, a gene expression by histotype interaction term was also included in the Cox proportional hazards model. Q values were computed to control for the false-discovery rate (24). Because of limited sample size, molecular subgroup analyses were not completed.

Data from eight EOC studies where more than one histotype group were included from the curated OvarianData (9, 25) were used to assess replication of key findings. The TCGA ovarian cancer study, which included only serous EOC tumors, was not included as one of the replication studies. We did not restrict serous tumors to high-grade due to the limited sample size and lack of information related to grade in some of the replication studies. A summary of these studies are presented in Supplementary Table S3. These eight studies consisted of the following studies in GEO (Gene Expression Omnibus): GSE51088 (26), GSE44104 (27), GSE26193 (28), GSE14764 (29), GSE2109 (IGC's Expression Project for Oncology), GSE30161 (30), GSE6008 (31–34), GSE9891 (5). Five of these eight studies also had clinical outcome(s): GSE51088 (overall survival, OS), GSE26193 (OS and PFS), GSE14764 (OS), GSE30161 (OS and PFS), GSE9891 (OS and PFS). The gene expression data in the GEO studies were all generated from microarray technologies, whereby all the datasets have been previously normalized for analysis, as described in Waldron and colleagues (9) and Ganzfried and colleagues (25).

Statistical analyses were completed in a similar manner as completed for the Mayo Clinic transcriptome study. We used a rigorous threshold for determining replication of histotype-specific differentially expressed genes, whereby a gene would be initially considered as “replicated” if more than one replication study had a KW P < 0.001. As the KW test does not look at direction of the effect, we next looked to verify similarity of gene expression pattern for the three histotypes between the discovery set (Mayo Clinic) and the replication set (GEO studies). Replication for PFS was considered if at least one of the 5 GEO studies with outcome data had an association P < 0.05.

We completed robust normalization of RNA-Seq gene expression data measured on 186 Mayo Clinic EOC patient tumors, resulting in 22,234 Ensembl genes suitable for differential expression analyses. SVA and normalization accounted for both known (i.e., tumor purity, library size, library prep) and unknown differences between samples. To determine the extent to which EC, CC, and HGSC histotypes clustered separately based on gene expression, unsupervised model-based clustering and PCA were performed, with results presented in Fig. 1. The three histotypes clearly separated from one-another (Fig. 1A), and the model-based cluster assignments agreed well with the histotypes. The HGSC tumors primarily grouped into two clusters (Cluster 2 and Cluster 5), whereas Cluster 1 primarily comprised EC tumors and Cluster 3 and Cluster 4 represented CC tumors (Fig. 1B; Supplementary Table S4). In comparison of the model based cluster groups with previously reported subtypes related to outcome (i.e., CLOVAR signature by Verhaak and colleagues; ref. 6), we observed a moderate association between the two different cluster assignments (all EOC tumors P value = 0.01; HGSC tumors only P = 0.35; Supplementary Table S4). In particular, Cluster 4 is enriched for “differentiated” (DIF) tumors (7 out of 9 or 78% of EOC tumors in Cluster 4 were denoted as DIF tumors), with the other cluster groups having moderate representation of EOC tumors across the four CLOVAR subtypes. However, it should be emphasized that the motivating questions and the approaches for generating the two cluster assignments were different. Measures of global expression, defined as the median gene expression level, differed between the three histotypes (P = 2.7 × 10⁻⁹; Fig. 2A and B), but were not related to PFS. However, median global gene expression was found to be associated with debulking status (P = 0.004) with higher expression observed in tumors that were optimally surgically debulked (Fig. 2C), noting that debulking status was not associated with histotype (P = 0.1, Supplementary Table S1).

We identified 1,254 individual genes that showed differential expression across the EOC histotypes (KW P < 1 × 10⁻¹⁰; Supplementary Table S5); 606 showed differential expression between HGSC and EC (474 upregulated in EC) and 244 genes showed differential expression between CC and HGSC (238 upregulated in CC) based on pairwise histotype comparisons (Wilcoxon P < 1 × 10⁻¹⁰). Only four genes were upregulated in HGSC compared with both EC and CC tumors: three protein-coding genes (GFRA3, ASCL5, and NPTX1) and a lncRNA (ERVH45-1). It should be noted that these three pairwise histotype comparisons have varied power to detect true differences due to the differences in the sample size in each histotype. Analysis involving the eight GEO publically available EOC studies, 73 of the 1,254 differentially expressed genes showed replication, defined as KW P < 0.001 in more than one study (Supplementary Table S6). Plots of the mean expression levels for the histotypes for these 73 genes can be found in Supplementary Fig. S1, with the top plot displaying the data from the Mayo Clinic study (discovery study) and the bottom plot showing the data from the various GEO study (replication studies). On the basis of this visual assessment, we determined that 32 protein-coding genes showed consistent direction of effect in the discovery and replication studies. The gene expression data from the Mayo Clinic study for these 32 genes are presented in Table 1, Fig. 3 and Supplementary Fig. S2, which include the following biologically relevant genes: CDH1, DRAM1, ENTPD5, FERMT2, GLRX, GPX3, KIAA1324, LAMB1, MAP2K6, PRKAB1, RHOB, SLC22A5, and TSPAN1. In addition, 6 of the genes denoted in Table 1 (ANXA4, FLOT1, GLRX, LAMB1, MITF, and SGPP1) were also found to be differentially expressed between histotypes by Zorn and colleagues (35).

Analysis of PFS with gene expression levels in the Mayo Clinic discovery study revealed that expression of BNIP3P17 was prognostic in all three histotypes and in the same direction (i.e., additive effect) and indicated that 22 additional genes were associated with PFS in a different manner by histotypes (i.e., interaction effect between histotype and gene expression; P < 0.0001; Supplementary Table S7). Therefore, we assessed the association of clinical outcome with these 23 genes in the five publically available GEO studies in which clinical outcome were reported; results are shown in Supplementary Table S8. It should be noted that many of the replication studies had small sample sizes and the model was only adjusted for stage (when available), as none of the studies had debulking status recorded. Association between PFS and BNIP3P17 did not replicate in any of the studies. However, the following nine genes showed evidence of a histotype-by-gene expression interaction effect (P < 0.05) replication in at least one replication study (Table 2; Fig. 3): CCNB2, CORO2A, CSNK1G1, FRMD8, LIN54, LINC00664, PDK1, PEX6, and USP31. Figure 4 presents the predicted PFS curves for these nine genes for the Mayo Clinic study for each histotype by high or low gene expression level. In particular, we observed that for 8 of these 9 genes (i.e., all but USP31), lower expression levels were associated with better outcome (longer PFS times) for EC and CC EOC patients, whereas the converse was observed for HGSC EOC patients (i.e., higher expression levels associated with longer PFS time), after adjusting for age of diagnosis, stage and debulking status. However, one of these eight genes, PEX6, was also found to have the expression level associated with disease stage (P = 0.001, Spearman correlation = −0.24).

Little progress has been made in the treatment of EOC, due in part to incomplete knowledge about disease heterogeneity with regards to histotype. Using the largest known next-generation sequencing study of the three most common histotypes, we sought to determine transcriptomic differences between EC, CC, and HGSC. We found 1,254 genes with P < 1 × 10⁻¹⁰ (FDR q < 3 × 10⁻¹⁰), of which 32 replicated with consistent direction of effect in publicly available data from GEO (Table 1). In addition, we replicated the association of nine genes with differing relationships with PFS between the three histologies (Table 2).

Overall, we found that the EOC tumors clustered together based on histotype. In particular, we observed most EC tumors (78%, 43 out of 55 EOC tumors) cluster together in one cluster (C1), whereas CC tumors cluster evenly into two clusters (C3 and C4) with 95% of CC tumors falling into C3 and C4 cluster groups. Finally, we observed that clusters 2 and 5 (C2, C5) are predominately HGSC clusters, with 84% (44 out of 52) and 93% (66 out of 71) of tumors in C2 and C5, respectively, being HGSC (Supplementary Table S4, Fig. 1). We also observed globally, as measured by median expression level, that CC tumors had higher global expression levels, followed by EC and HGSC (P = 2.7 × 10⁻⁹); optimally debulked patients had higher transcriptomic expression levels compared with tumors from non-optimally debulked patients (P = 0.0044; Fig. 2). This observed difference in global gene expression could be due to copy number and genomic instability, known to be a hallmark of serous EOC and triple-negative breast cancers (36). That is, the observed lower levels of median global gene expression in HGSC tumors could be do the higher genomic instability compared with the other tumor histotypes of EOC (i.e., higher genomic instability might lead to fewer “functional” genes and thus a lower level of global gene expression levels). Much less is known about how this compares in the rarer histotypes of ECs and CCs, with some evidence that ECs and CCs have lower levels of genomic stability (12). We did not observe any association of median global gene expression level with PFS (P = 0.32), after adjusting for debulking status, stage and histotype, and no observed histotype-specific effects for association of median global expression level with PFS (interaction effect, median global gene expression by histotype P = 0.73).

As the primary focus of this research, we identified 32 genes that were differentially expressed across histotypes, along with another nine genes associated with PFS differently between the three histotypes (i.e., interaction effect between gene expression level and histotype with respect to relationship with PFS). Out of the 32 replicated differentially expressed genes, 29 genes were significantly upregulated (Wilcoxon P < 0.0001) in CC compared with HGSC and EC (Table 1, Supplementary Fig. S2). KIAA1324, important in cellular response to stress, and MAP2K6, involved in cell-cycle arrest and apoptosis, were found to be upregulated in EC compared with both HGSC and CC tumors, whereas FLOT1, SLC3A1, and C12orf75 were found to be downregulated in EC tumors (Wilcoxon P < 0.0001). Mitogen-activated protein kinases (MAPK) have been found to be metastasis suppressor proteins (37); therefore, high expression of MAP2K6 observed in EC tumors might be related to reduced metastases and subsequent longer survival. Recently, miR-625-3p has been identified as a biomarker for oxaliplatin resistance in colon cancer, where MAP2K6 is the direct target for this microRNA (miRNA). The subsequent down-regulation of MAP2K6 impairs p38-MAPK stress signaling and thus leads to oxaliplatin resistance by preventing apoptosis (38). As most EOC patients are treated with combination platinum–taxane chemotherapies, a similar mechanism might underlie differences in clinical outcome, whereby EC EOC patients generally have better clinical outcome as compared with HGSC patients. Finally, TALDO1 and UPK3B were downregulated in HGSC compared with EC and CC tumors; little is known about how these two genes relate to carcinogenesis and EOC.

One limitation of replication study was the inclusion of both low- and high-grade serous tumors due to the limited sample size or lack of grade information in many of the studies. As the majority of serous EOC tumors are high-grade, with only 5% to 10% of serous tumors being low grade (39), we believe the impact of including possible low-grade serous EOC tumors in the replication study to be small. A second limitation of the replication study is that the majority of studies involved patients enrolled before the start of the enrollment for the Mayo Clinic study (Supplementary Table S3). A final limitation of the study is the focus on only transcriptomic differences between the three histologies, without orthogonal validation of downstream proteins known to be discriminant of the different histologies. None of our 32 validate genes included genes known to be related to histotype based on immunohistochemistry (IHC) studies (VIM, PGR, ARID1A, NAPSA, TP53, CDKN2A, TFF3, WT1; ref. 40); however, all but TP53 had a KW P < 0.05 for association with histotype (Supplementary Table S9).

Of particular interest among the nine PFS-related genes, two genes (LIN54 and CCNB2) are related to cell-cycle and two genes (CSNK1G1 and PDPK1) are kinase-related genes. For all nine genes, with the exception of USP31, we observed an increased risk of progression in EC patients with higher gene expression levels (Table 2; Fig. 4). To a lesser extent, we saw a similar association between high gene expression and increased risk of progression in CC for these genes. In contrast with EC and CC, high gene expression for these genes, with the exception of USP31, was associated with better clinical outcome in HGSC. For example, EC and CC patients with higher levels of CCNB2 had worse outcome than patients with lower expression. This finding is not unique to EOC, as it has also been observed in breast cancer (41). In contrast, HGSC patients with high expression of CCNB2 had a reduce risk of progression (41). Although elevated PDK1 has been found to promote tumor growth and metastasis in breast cancer (42), HGSC patients with high expression of PDK1 have been found to have improved survival compared with patients negative for protein expression (43), supporting our findings that HGSC patients with higher PDK1 expression have better outcome compared with patients with low PDK1 expression (Fig. 4). Interesting, we observed the opposite effect of PDK1 for patients with CC and EC tumors, with low expression levels resulting in better clinical outcome. Finally, CSNK1G1 (i.e., casein kinase I, gamma 1 or CK1γ1) is involved in many cancer related pathways, including hedgehog signaling, mTOR signaling and p53 interaction pathway, and as such, many CK1 inhibitors are currently under investigation (44, 45).

In conclusion, we present replicated findings on transcriptomic differences and how they appear to relate to PFS for two of the rarer EOC histotypes of EC and CC, along with comparison to the common histotype of HGSC. One of the strengths of this report is the incorporation of eight additional public expression datasets for replication. Future research is needed to assess these results in a larger cohort of EC and CC patients, along with the incorporation of epigenetic, copy number and protein expression information to determine molecular differences between the histotypes of EOC and their impact on clinical outcome. These findings suggest important biological insights into the differences between the three most common histotypes of EOC and may aid in the development of targeted treatment options, particularly those related to the treatment for patients with EC and CC EOC.

S.J. Weroha reports receiving commercial research funding from Genentech, Novartis, and Tesaro and is a consultant/advisory board member of KIYATEC. No potential conflicts of interest were disclosed by the other authors.

Conception and design: B.L. Fridley, S.A. Gayther, E.L. Goode

Development of methodology: B.L. Fridley

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): R. Raghavan, X. Hou, S.J. Weroha, K.R. Kalli, J.M. Cunningham, K. Lawrenson, E.L. Goode

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): B.L. Fridley, J. Dai, R. Raghavan, Q. Li, S.J. Weroha, C. Wang, J.M. Cunningham

Writing, review, and/or revision of the manuscript: B.L. Fridley, J. Dai, Q. Li, S.J. Winham, S.J. Weroha, C. Wang, J.M. Cunningham, S.A. Gayther, E.L. Goode

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): B.L. Fridley, R. Raghavan, X. Hou

Study supervision: B.L. Fridley

This research was supported in part by the University of Kansas Cancer Center (P30 CA168524) through a summer student research project award (to J. Dai), The Mayo Foundation (to B.L. Fridley; sequencing), R01 CA122443 (to E.L. Goode), P50 CA136393 (to K.R. Kalli, E.L. Goode, and B.L. Fridley), P20 GM103418 (to B.L. Fridley and Raghavan), R21 CA182715 (to B.L. Fridley), and R00CA184415 (to K. Lawrenson).

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.