Outcome-Related Differences in Gene Expression Profiles of High-Grade Serous Ovarian Cancers Following Neoadjuvant Chemotherapy

High-grade serous ovarian cancer (HGSC) is a poorly differentiated, fast growing subtype of ovarian cancer that accounts for 70% to 80% of ovarian cancer–related deaths with little improvement in long-term survival in the past 30 years (1). The use of neoadjuvant chemotherapy (NACT) followed by interval cytoreduction is increasing, because upfront complete cytoreductive surgery is often difficult to attain (2). An estimated 40% of women with HGSC undergo NACT as a primary treatment strategy, although numbers vary depending on the institution and country (2–4). It is increasingly accepted that some patients benefit from this treatment, but because most patients are cytoreduced to no visible disease in the NACT setting (5), other biological factors play a role in determining outcome for patients who underwent NACT.

Given that approximately 80% of HGSCs arise outside the context of a known hereditary predisposing factor, uncovering the somatic drivers of the disease is a priority (6). The Cancer Genome Atlas (TCGA) project revealed that somatic mutations in HGSC are generally confined to very specific genes and pathways, especially the homologous recombination DNA repair (HRD) and cell-cycle pathways (7, 8). Such studies have led to an increased understanding of response and resistance to treatment in the first line setting (5, 9–11). Other studies have similarly drawn a parallel between gene expression and clinical outcomes (12, 13). However, the majority of genomic studies of HGSC have excluded samples from NACT-treated, or chemotherapy-exposed patients, and only focused on patients undergoing primary debulking surgery (PDS) followed by adjuvant chemotherapy. Patients undergoing NACT are clinically selected due to the extent of their disease, and have been understudied compared with patients undergoing PDS. One group suggested that BRCA1-mutated cells are preferentially eradicated during NACT, leaving HRD-proficient cells unscathed (14). The same group reported changes in TP53 mutation status in the residual disease as well. Another group reported changing levels of tumor-infiltrating lymphocytes (TIL; ref. 15), whereas the most recent study on pre- and post-NACT samples showed significant changes in the expression of numerous genes (16). However, none have included an analysis integrating the molecular features of post-NACT samples with patient outcome. Understanding the molecular drivers of the residual tumor cells in HGSC after NACT exposure may allow us to better understand chemo-resistant HGSC to enable the development of tailored treatment options following NACT.

The aim of this pilot study was to characterize the transcriptional landscapes of residual tumor after NACT in patients with HGSC and correlate those features with patient outcome.

Patients considered for the study were diagnosed and treated for HGSC at the Jewish General Hospital between 2003 and 2015 (Supplementary Table S1). Tissue and blood samples were collected at the time of surgery and stored in the gynecologic oncology tumor bank (protocol #03-041). This study was approved by the Jewish General Hospital Research Ethics Board (protocol #15-070). All patients participating in biobanking and research activities gave informed written consent. Sixty-five samples were selected for this study (NACT n = 57, paired pre-NACT biopsies n = 8). All samples were selected after histologic confirmation and high tumor content (>70% tumor content) by a gynecologic pathologist (A. Ferenczy and M. Pelmus). For each patient, information such as age, body mass index (BMI), histologic type, tumor grade, FIGO stage, extent of cytoreduction, and chemotherapy treatment history, was collected in a prospectively maintained clinical database. During the post-surgical surveillance period, follow-up examinations were performed at four-month intervals during the first 2 years from diagnosis, at 6-month intervals during the fourth and up to the fifth year, and yearly thereafter.

Progression-free survival (PFS) was defined as the time from diagnosis to evidence of recurrence by imaging, or death. Overall survival (OS) was defined as the time from diagnosis to death, or the latest follow-up date. Platinum resistance was defined as progression within 6 months after the last cycle of platinum-based chemotherapy. Responses were assessed using the definitions provided by the Gynecological Cancer Intergroup (17).

Before DNA/RNA extraction and inclusion in the study, 12-μm sections were cut from each sample and subsequently stained with hematoxylin and eosin (H&E). Slides were reviewed by a gynecological pathologist (A. Ferenczy and M. Pelmus) and samples with low cancer cell content were excluded. RNA and DNA were serially extracted from the selected samples using the Bio Basic “All-in-One” DNA/RNA/Protein Mini-Prep Kit (Bio Basic Inc.). Concentrations for each sample were measured using the NanoDrop ND-100 spectrophotometer (NanoDrop Technologies) and Qubit (Thermo Fischer Scientific). Biopsy samples represent available samples from patients with at least 2 years of follow-up from the time of diagnosis that passed DNA/RNA quality controls.

Gene expression was measured using the NanoString PanCancer Pathway Panel (NanoString Technologies). The panel contains probes for 770 genes implicated in carcinogenic pathways, curated from data from TCGA. Raw data were processed and normalized using standard protocols suggested by NanoString using the nSolver Analysis Software (NanoString Technologies) provided. Detailed protocols and normalization procedures can be found in the Supplementary Methods.

For DNA analysis, sequencing was performed on the Illumina MiSeq (Illumina Inc.). The Roche Nimblegene TruSeqLT (Roche) preparation kit was used to create a library of 400 targeted regions in 168 genes of interest (List in Supplementary Table S2). Results were aligned to the hg38 UCSC Genome Browser assembly (NCBI Assembly ID: GRCh38, GCA_000001405.15). The alternative allele frequency was estimated using the somatic variant caller provided by Illumina in the “MiSeq Reporter Software.” If multiple variants were identified at a single locus, the variant with the maximum frequency value was used. The average sequencing depth was 687. All sequencing reads within two standard deviation from the mean total were removed. Reads with an estimated Quality Score <30, which is a prediction of a false-positive call by the Reporter Software of 1 in 1,000, were filtered out. We performed principal component analysis (PCA) for all of the alternative allele frequencies for the variants remaining at this point. We used the eigenvector of the first principal component to create a tumor content value and then compared the value of each sample with one another, as described in a previously published method (18). Stromal tissue samples should have significantly different somatic mutation allele frequencies than tumor tissue samples, and the differences should be exposed using PCA. The tumor content variable was also used to control the effect of non-tumor cell content on the gene expression results.

To determine the most likely pathogenic mutations that may affect disease progression, the Ensembl Variant Effect Predictor (19) was used to annotate the remaining mutations. Alleles with a population allele frequency below 1.5% in the gnomAD (20) database were kept for further downstream analysis based on the assumption that rare alleles are more likely to be pathogenic. For selected mutations the raw BAM files were visualized using the Integrated Genome Viewer (21). Synonymous and intronic mutations were removed, unless the later occurred within three base-pairs of a coding exon that could result in a splice-region mutations. The potential pathogenicity of missense mutations was assessed using the following in silico computational tools: PolyPhen-2 (22), Sift (23), M-CAP (24), MutationAssessor (25), and REVEL (26). Mutations of interest were those predicted to be pathogenic by at least three of the five predictors, as previously described (27). Frameshift mutations were assumed pathogenic, unless proven otherwise in the literature. Because of limitations in the sequencing panel design, MiSeq was performed on TP53 only at exons 5 to 9. All other exons were evaluated by established PCR-based assays followed by Sanger Sequencing using methods and primers described in the IARC TP53 Database (28).

To discriminate molecularly relevant survival groups, we used a previously described semi-supervised 5-fold cross-validation classification technique based on the gene expression data (29–31). In summary, the full dataset F was randomly separated into five approximately equal folds (F1,…, F5). Using data from a training set T1 = F – F1, variables were selected by fitting a univariate Cox proportional hazard for all genes to OS. The expression values of genes with a P value of lower than 0.003 in the univariate analysis were then fitted to T1 using a multivariate Cox proportional hazard model, and the resulting coefficients were used to predict the survival risk in patients within F1. The P value threshold choice captures approximately 1 to 10 genes in each of the five fitted prediction models. These genes' P values are generally smaller than the point of inflection in each QQ-plot (Supplementary Fig. S1) supporting their selection. The process was reiterated for each of the five folds, where T2,…,5 = F – F2,…,5. Each case was thus assigned a risk score using a predictive model using gene expression from which it was absent. All cases were then grouped and assigned to a predicted worse (n = 34) or a predicted better outcome group (n = 23) based on their risk score in a way that optimized the separation of the computed Kaplan–Meier curves between the two risk groups. Clinical data were integrated using Cox Proportionate Hazard multivariate models. Differential gene expression between the two survival groups was performed using the negative binomial generalized linear model and controlled for tumor content using the tumor content variable. P values were adjusted by the Benjamin–Hochberg FDR method. Relevant molecular pathways were uncovered using Gene Set Enrichment Analysis and the Reactome Pathway database (32, 33). The difference in gene expression between biopsy and post-NACT samples were compared using paired t tests. All data manipulations were performed using the R statistical software (34).

Patients were discriminated into either the predicted better or predicted worse outcome group as previously described in the methods. The clinical characteristics between the two predicted outcome groups were not significantly different, except for second line platinum resistance, where more patients from the worse outcome group were platinum resistant following second line treatment with platinum-based chemotherapy (Table 1). The outcome groups did not have significantly different PFS (Fig. 1A; log-rank test P = 0.26, median PFS = 14.0 months vs. 15.8 months), but did have significantly different OS (Fig. 1B; log-rank test P = 0.0019, median OS = 31.8 months vs. 48.6 months), This observation is not surprising given the similarities in response rates to platinum-based therapy in the first line setting between the two groups, followed by an increase in platinum-resistance in the worse outcome group after second line therapy. When controlling for clinical variables in a multivariate analysis, the outcome group division remained significantly associated with OS [Table 2; HR, 3.16; 95% confidence interval (CI), 1.56–6.40; Wald P < 0.001]. Surgical outcome (HR, 2.57; 95% CI, 1.25–5.28; Wald P < 0.01) and age (HR, 1.04; 95% CI, 1.01–1.06; Wald P < 0.007) at diagnosis were also significantly associated with OS.

A differential expression analysis between the two outcome groups was performed. Out of the 770 genes profiled for expression using the NanoString panel, 42 differentially expressed genes were found between the two groups at an FDR below 10% (Fig. 2A–B; Supplementary Tables S3–S4), where 11 genes were overexpressed in the better outcome group and 31 were overexpressed in the worse outcome group (Fig. 2B–C). The majority of differentially expressed genes exhibited a >2-fold change in gene expression (Fig. 2A). Figure 2C illustrates the relative differences and overlap between the expression of the differentially expressed genes in the two survival groups. Gene Set Enrichment Analysis showed significant expression of cell cycle and DNA repair pathway genes in the worse outcome group (Fig. 2D). The core enriched genes uncovered in the analysis tend to belong to multiple of the top five enriched pathways, which is consistent with the significant interrelations between the cell cycle and DNA repair pathways (Fig. 2E).

We had the opportunity to investigate eight fresh-frozen biopsy samples obtained before NACT, together with the clinical paths of each patient for whom the biopsy sample was available. An analysis of the differences in gene expression between the chemo-naïve pre-NACT biopsy samples and the chemo-treated post-NACT sample revealed that 12 of the 42 genes previously found to be differentially expressed between our two outcome groups showed a statistically significant change in direction of gene expression (Fig. 3). BRIP1, CDC25A, CHEK1, EPOR, CCNE1, and SKP2 expression levels were lower in the post-NACT sample in the better outcome group only. CCND2, ITGA7, and LTBP1 expression levels increased in the better outcome group only, whereas H2AFX and E2F1 levels lowered in the worse outcome group only. CCNA2 expression levels were significantly lower in both outcome groups. All changes in gene expression were consistent with the difference shown in the differential expression analysis, which may suggest that, for those genes, the difference in gene expression and their relation to outcome may be due to chemotherapy.

Out of the 168 genes assayed, the three most recurrently mutated genes in our patient cohort were TP53, BRCA1, and BRCA2, consistent with previous studies on the mutational landscape of HGSC (Fig. 4A; Supplementary Table S5). However, the frequency of TP53 mutations was lower than the 95% observed in previous reports (7). In the NACT worse outcome and NACT better outcome groups, the mutation rates were 70% and 59%, respectively, though this difference in frequency was not significant (χ² P = 0.41). Studying post-chemotherapy samples might raise the risk of sampling residual stromal tissue from the tumor rather than cancer cells as the latter is affected by treatment. Therefore, we created a tumor content variable by performing PCA on the alternative allele frequencies for all variants and ascribing the value of the first eigenvector to the new variable. Tumor samples should share similar somatic mutation loads, and samples with excessive stromal contamination and fewer somatic mutations as a result should stand out through this analysis When we compared the tumor content variables of our two groups, the difference was not statistically significant (Fig. 4B; Wilcoxon P = 0.063), suggesting that tumor purity was similar between the two groups. Similarly, there was no difference in tumor content value between TP53-mutated and TP53-WT samples, in either outcome groups (Fig. 4C), or TP53 gene expression values among the same groups (Fig. 4D). Note the similarity in TP53 expression patterns between both groups based on mutation effect: Missense variants exhibiting at higher TP53 expression value in contrast with those with a frameshift, splice site and nonsense variants (Fig. 4D). Our mutational profiling suggest that differences in gene expression profiles between the two outcome groups are unlikely due to excessive stroma tissue contamination.

To the best of our knowledge, our study is the first to report on the molecular features of post-NACT HGSC in the context of survival. Using high quality, post-NACT tumor samples from a single-institution biobank and its associated prospectively maintained clinical database, we uncovered two molecularly different subgroups associated with OS. Unsupervised clustering methods of gene expression values have previously been used with some success in identifying patient groups with distinct survival outcomes (12), but their clinical application remains elusive. Alternatively, supervised molecular analyses that use survival cutoff values to identify patient groups may obscure subtle molecular differences among tumors. Because patient survival is most likely due to a mix of clinical and molecular characteristics of the tumor, we applied a cross-validation approach that incorporated both gene expression and survival. This method was previously shown to provide a good balance between bias and variability (29, 30). Because of the small number of patients included in this pilot study, overfitting is a possibility and validation in independent datasets would strengthen our conclusions. However, no other studies have publicly available data on post-NACT patients and direct comparison of patients that received NACT with patients that were primarily cytoreduced would not be appropriate. Indeed, chemotherapy treatment is known to alter tumor cell populations in patients and thereby affect the expression profile of cancer cells, rendering direct comparisons with tumor samples from unexposed patients difficult (15, 16, 35). In addition, patients are clinically selected for NACT because of diffuse carcinomatosis, increased tumor burden, worse symptoms, and presumed poorer outcomes (36) and thus represent a different patient population compared with patients undergoing primary debulking who were not yet exposed to chemotherapy (36).

Differential expression analysis between the two survival groups revealed an increase in the expression of genes from the cell cycle and homology-directed DNA repair pathways. Aberrant cell cycle function is necessary for escaping the cell's checkpoints, evading apoptosis, and thus driving tumorigenesis. FOXM1, CHEK1 and CDC25A are important genes overexpressed in the worse outcome group and are interesting due to their roles affecting both the cell cycle and HRD pathways. CHEK1 relays the stress responses due to double strand breaks from ATR, arresting the cell cycle and allowing HRD-mediated repairs in the G₂ phase. Interestingly, many of the genes found in our analysis are under FOXM1 transcriptional control (Supplementary Fig. S2; ref. 37). FOXM1, similarly to CCNE1, drives the cell cycle by transcribing genes necessary for G₂–M phase transition. It has been linked to invasion and chemoresistance in HGSC (37, 38, 39). Together, FOXM1 and CHEK1 seem to possess antagonistic functions, but cell cycle regulation is complex and the interplay between the different proteins remains poorly understood. Interestingly, targeting FOXM1 and CHEK1 in HRD-proficient HGSC is increasingly being studied. In HGSC cell lines, FOXM1 downregulation was shown to confer increased susceptibility to cytotoxic agents (40). Moreover, due to its ubiquitous overexpression in HGSC, it may be an interesting biomarker if it does serve as an early indicator of tumor response through the use of NACT. In a recent Phase II study, CHEK1 and CHEK2 inhibition in HRD-proficient, heavily pre-treated HGSC cases showed promising clinical activity (41). On the basis of our results, such agents may represent early alternatives to traditional cytotoxic agents in patients with high expression of those genes after NACT and surgery, such as those identified in our worse outcome group that are more likely to be platinum resistant in the second line setting.

Arend and colleagues (16) recently reported on the differential expression of 20 paired pre- and post-neoadjuant HGSC samples. The genes most downregulated after NACT were uncovered in our analysis of pre- and post- NACT samples as well: E2F1, BRIP1, CHEK1, and CCNA2. In both studies, the downregulation was not necessarily seen among all patients. Our analysis suggests that it may be featured most prominently in patients with better OS. This is consistent with previous observations that high expression of genes in the cell-cycle pathway and homologous recombination DNA repair pathway, to which these genes belong, are detrimental to patient outcomes (7). It may also suggest the possibility that differences in gene expression post- and pre-NACT may differ from patient to patient in a way that reflects outcome. In other words, the change in expression during treatment, in addition to baseline levels, may also contribute to outcomes. However, the sample size is too small in both studies to derive definitive conclusions. It reflects the difficulty in acquiring high quality pre- and post-NACT tumor samples from patients. Serial analysis of additional biopsies are invasive and post-NACT intraoperative samples often contain nodules of low cancer content due to chemotherapy-related tissue necrosis.

Contamination of RNA by stromal cells in the studied samples is a major issue in transcriptional studies. To further investigate this possibility, we assayed cancer cell content by multiple methods. First, gynecologic pathologists estimated the cancer cell percentage for each sample before processing, choosing only samples with a proportion of cancer cells above 70%. Second, we created a tumor cell content covariate using mutational data that allowed us to compare cancer cell content between our two outcome groups to insure that our results were due to changes in gene expression, not sample purity. Third, because most HGSC are TP53 mutated or deficient in some way, we compared the levels of TP53 expression. The similarity in TP53 expression patterns between the two groups are in keeping with known mutation effect on gene expression: Missense variants usually exhibit higher TP53 expression value due to stability of mutant mRNA in contrast to frameshift, splice site and nonsense variants, which are often degraded due to nonsense-mediated decay (42, 43). Alternatively, the similar tumor cell content covariate and TP53 expression between both groups could also suggests that, in the event of contamination, it occurred equally in both groups and should not affect differential expression results.

Both outcome groups had a lower TP53 mutation rate than expected on the basis of other large scale studies of tumor samples primarily derived from chemotherapy-naïve patients (7). This may be due to our targeted approach to sequencing that specifically targeted exons 5–9 of TP53 by MiSeq, and of the other exons by PCR amplifications. The latter is less reliable method to detect low allele frequency mutations (42, 43). The majority of TP53 mutations in HGSC can be found in exons 5–9 (43), but many may have been missed in the exons not covered by MiSeq analysis. Other cellular processes, such as methylation (44) and EME2 overexpression (45) may also be responsible for a TP53-null phenotype in the absence of detectable mutations. However, because the vast majority of mutations are point mutations, reversions are possible and the lower frequency of mutations seen in the post-NACT group remains an interesting observation. Reversion of TP53 mutations has previously been observed in a study looking at matched pre- and post-NACT breast cancer tissue, and was associated with a better survival (46). Similarly, in one of the only studies of pre- and post-NACT HGSC, loss of TP53 mutations was observed in 5 out of 11 initially TP53-mutated samples suggesting a restoration or selection of tumor cells with wild-type alleles post-NACT (14). Interestingly, gain of TP53 mutations was also observed in 3 out of 11 TP53-wild type samples. This may be explained by the clonal heterogeneity of HGSC, and the Sanger sequencing methods used by the authors. Nonetheless, TP53 mutation status, which is a defining feature of primary, pre-chemotherapy derived HGSC to many, may play a role in the chemotherapy induced selection during NACT and should be further investigated.

Although the complete absence of disease on pathological evaluation after NACT is a rare occurrence, it is common to observe no macroscopically visible disease at the time of surgery in patients who respond remarkably to chemotherapy. One possible explanation may be that BRCA1/2-mutated cells may be selectively killed during NACT, especially if the mutation is accompanied with somatic loss of the WT allele as suggested by previous reports (14). Therefore, our cohort may be enriched for patients who do not respond well to chemotherapy in the first place. Patients with post-NACT, debulked residual disease exhibiting lower levels of cell cycle and HR genes, especially when compared with the levels before chemotherapy, may have a lowered risk of early death compared with those patients in which few molecular changes had occurred. In the latter population, the differences in expression of cell cycle and HR genes may be mediated by key genes, such as FOXM1 and CHEK1, for which known antagonists exist (40, 41). Future studies should strive to include pre-NACT biopsy samples for all patients. On the basis of our current understanding, the cancer cells remaining after NACT are the most likely candidates to be responsible for recurrences in the future, and a better understanding of this residual cancer cell population might help in the selection of the most appropriate therapies.

We recognize the limited scope of our analysis, but we hope that the results from this pilot study will encourage the longitudinal assay of the progression of HGSC and promote the use of survival data when analyzing the temporal molecular landscape of the disease.

W.H. Gotlieb is an advisory board member, reports receiving a commercial research grant, and has received speakers bureau honoraria from AstraZeneca. No potential conflicts of interest were disclosed by the other authors.

Conception and design: D. Octeau, R. Kessous, L. Kogan, C.M.T. Greenwood, S. Salvador, S. Lau, P.N. Tonin, A. Yasmeen, W.H. Gotlieb

Development of methodology: D. Octeau, R. Kessous, K. Klein, C.M.T. Greenwood, S. Salvador, P.N. Tonin, A. Yasmeen, W.H. Gotlieb

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D. Octeau, R. Kessous, L. Kogan, L.C. Van Kempen, S. Salvador, S. Lau, W.H. Gotlieb

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Octeau, R. Kessous, K. Klein, L. Kogan, C.M.T. Greenwood, L.C. Van Kempen, S. Salvador, S. Lau, P.N. Tonin, W.H. Gotlieb

Writing, review, and/or revision of the manuscript: D. Octeau, R. Kessous, K. Klein, L. Kogan, A. Ferenczy, C.M.T. Greenwood, L.C. Van Kempen, S. Salvador, S. Lau, P.N. Tonin, A. Yasmeen, W.H. Gotlieb

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D. Octeau, P.N. Tonin, W.H. Gotlieb

Study supervision: S. Salvador, P.N. Tonin, A. Yasmeen, W.H. Gotlieb

Other (selection of slides and review of pathological features of ovarian cancers submitted for this study.): M. Pelmus

This work was supported by grants from AstraZeneca Canada Inc., the Israel Cancer Research Fund, Gloria's Girls of the Jewish General Hospital Foundation, the Anne-Marie and Mitch Garber Fund and the Susan and Jonathan Wener Fund.

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.