The tumor material used in the present study was retrieved from three different cohorts and is referred to throughout as cohorts 1, 2, and 3. Cohort 1 comprises 72 ovarian carcinomas of mixed histopathologic subtypes, partly outlined in Jönsson et al. (29). The tumors in cohort 1 were collected among Swedish and Danish Lynch syndrome carriers (n = 28) and Swedish sporadic cases (n = 44). Formalin-fixed, paraffin-embedded (FFPE) tissue was available from all tumors. Cohort 2 consists of 11 OCCCs, 9 of which were obtained from the Skåne University Hospital (Sweden) ovarian tumor biobank, and 2 of which overlapped with cohort 1. FFPE tissue was available from all tumors. Cohort 3 consists of 43 OCCCs from Skåne University Hospital arranged in a tissue microarray (TMA). Clinicopathological data for all cohorts are summarized in Table 2. All tumor samples were collected at primary surgery, and the patients had not received prior chemotherapy. Histopathological subtype and grade were determined according to Silverberg and the WHO 2003 classification, which was used at the time of diagnosis (30, 31). For tumors in cohorts 2 and 3, the histopathological subtype was determined according to the more recent WHO 2014 classification (32). All tumors were staged according to the International Federation of Gynecology and Obstetrics criteria (33). Ethical approval for the study was granted from the ethics committee in Region Hovedstaden, Denmark, and from the Lund University ethics committee, Sweden, and all patients provided written, informed consent.

RNA was extracted from all FFPE samples in cohort 1 (n = 72 ovarian carcinomas). As outlined previously, non-necrotic tumor areas with >70% tumor cell content were selected from three to five 10 µm sections of Hematoxylin and Eosin stained tissue per tumor and RNA was extracted from these areas using the High Pure RNA Paraffin Kit (Roche, Castle Hill, Australia) (29). RNA concentrations were determined using a NanoDrop Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA); 68/72 tumors met the quality criteria (300 ng of RNA with a 260/280 ratio >1.8). Gene expression analyses were performed at the SCIBLU Genomics Centre, Lund University, Sweden. The cDNA-mediated Annealing, Selection, extension and Ligation (WG-DASL) assay (Illumina Inc., San Diego, CA, USA), containing 24,526 probes representing 18,626 unique genes, was used for whole-genome expression analysis. The samples were randomized on the chips and were profiled following the manufacturer’s instructions. BeadChips were then scanned on a BeadArray™ Reader using BeadScan software (v4.2), during which fluorescence intensities were read and images extracted. A raw average signal intensity >250 for each probe and >8,000 detected genes/sample were required for further analysis; 67/68 samples met these criteria. The raw data are freely available in NCBI’s Gene Expression Omnibus (34) through GEO series accession number GSE37394. The raw data were quantile normalized, log2 transformed and subjected to a presence filter of 80% across probes, with a detection P-value ≤0.01, leaving 12,747 probes. Biological duplicates were averaged and when multiple probes identified the same gene, the probe with the highest variance across tumors was chosen to represent the gene, resulting in 10,000 remaining probes. The data were then mean centered across samples using R version 3.2.2 (35). Thereafter, the data were imported into the MeV version 4.9.0 software (36) and a multiclass significance analysis of microarrays (SAM) analysis (37) using the histological subtypes as groups (serous, EM, MUC, and clear cell) was performed. A SAM analysis, using 1,000 random permutations, was used to identify the 5% most significant differentially expressed genes between subgroups, with a false discovery rate (FDR) <5%. OCCCs were compared to the other subtypes individually using the same settings. Pearson correlation was used for distance metric selection and complete linkage was used for cluster generation. Enrichments of gene ontology (GO) terms were explored using the online Panther gene list analysis tool with default settings for Homo sapiens, using a binomial statistical overrepresentation test and Bonferroni correction for multiple testing with the 10,000 genes as background (38). GO terms were compared using the online Gene Semantic Similarity Analysis and Measurement tool (G-SESAME) with default settings (39), generating a semantic similarity score between GO terms of 0 (low similarity) to 1 (high similarity). Gene set analysis using gene sets from the REACTOME database [1,764 pathways, ConsensusPathD (40)] was performed using the GS-Reg package version 1.8 in R with 1,000 permutations, a minimum gene count of 7 and a P-value cutoff of 0.001 for all genes (41). Genes from the SAM analysis were analyzed using the ConsensusPathDB, with a minimum gene count of 7 and a q-value cutoff of 0.05 for REACTOME pathways. Venn diagrams were generated using the online Venn Diagram tool from the Bioinformatics and Evolutionary Genomics web page.¹ Gene module score comparisons according to gene modules defined by Desmedt et al. were performed as previously described (42, 43). Statistical tests (Kruskal–Wallace and Dunn’s test) and boxplots were performed using R version 3.2.2.

DNA was extracted from all FFPE samples in cohort 2 (n = 11 OCCCs) using the AllPrep DNA/RNA FFPE Kit (Qiagen, Venlo, Netherlands). Two 10 µm tissue sections/tumor were used for DNA extraction and >1 μg DNA was obtained from all samples. Three micrometer hematoxylin stained sections, retrieved immediately before the sections used for DNA extraction, were reviewed by a gynecological pathologist (SWF) to establish tumor cellularity, which ranged from 65 to 99%, with a median of 90% (data not shown). DNA concentrations were measured using the Qubit Fluorometric Quantitation^® (Life Technologies, Thermo Fisher Scientific, Waltham, MA, USA), and a qPCR-based DNA quality control was performed using the Trusight Tumor Sample Preparation Kit (Illumina Inc.), with ΔCT values <6.5 considered sufficient. One sample was excluded due to poor DNA quality. For the remaining 10 samples, 600 ng genomic DNA/sample (50 ng/µl) was used for targeted DNA sequencing. Sample preparation and sequencing was performed at Oxford Gene Technology™ (OGT, Oxford, Great Britain). The Next-Generation Sequencing SureSeq™ solid tumor panel, a hybridization-based enrichment tool targeting all codons of the exons in the included genes, was used. The gene panel consists of 60 key cancer genes and has been validated for research use on FFPE samples² (Table S1 in Supplementary Material).

Alignment was performed using the human build GRCh37 reference genome and Burrows–Wheeler Aligner (version 0.7.10). Local realignment was carried out using the Genome Analysis Tool Kit (GATK, version 1.6), and duplicate reads were removed using Picard version 1.107. Variant calling was performed using the BCbio-nextgen tumor only prioritization pipeline³ with Mutect2 from GATK as variant caller and default settings (44). Variants were annotated using the Variant Effect Predictor (45). Variants were classified using the GEnome MINIng (GEMINI) framework to filter out possible germline variants (46), with a variant allele frequency filter cutoff of >5%.

Immunohistochemical and CISH evaluation of ERBB2 were performed using standard procedures for HER2 testing in breast cancer. The anti-HER2/neu rabbit monoclonal antibody 4B5 (Ventana, Tucson, AZ, USA) was used for IHC and the 2,4-dinitrophenyl (DNP)-labeled HER2 probe and digoxigenin (DIG)-labeled Chr17 were used for CISH. The HER2 gene and Chr17 signals were detected using the ultraView SISH DNP and ultraView Red ISH DIG detection kits (Ventana). Tumors with an IHC score of 3+ and/or a gene/centromere ratio ≥2 were considered positive (Figure S1 in Supplementary Material). Anti-AKT-pS473 (Dako, Denmark) staining was performed using Dako EnVision™ System/HRP. A 1:10 dilution was used with Dako Retrieval Solution (pH 9). Any staining was considered positive.

We wanted to investigate the gene expression profiles of OCCC compared to the other ovarian cancer subtypes. In a previous study of Lynch syndrome-associated ovarian cancer, we reported a distinct gene expression pattern of OCCC compared to other subtypes using unsupervised hierarchical clustering. This was found to be irrespective of Lynch status (29). Hence, we performed a supervised multiclass SAM analysis without a variance filter using all 10,000 genes and histological subtypes as groups, thereby identifying 505 differentially expressed genes (FDR < 5%). OCCCs were found to display a distinct transcriptional pattern compared to all the other histological subtypes (Figure 1A). The pathway analysis revealed four altered pathways (Table 3), with the extracellular matrix (ECM) organization, axon guidance, and developmental biology pathways exhibiting many genes in common. We next applied a variance filter to select the 2,500 genes with the highest variance across tumors. This was done to select as many targets as possible throughout the signaling network hierarchy, as proposed by Komurov and Ram (47). This SAM analysis identified 504 significant genes (FDR < 10%), and the subsequent pathway enrichment analysis revealed genes belonging to both the NCAM signaling and TFAP2A transcription pathways, as well as pathways involved in MAPK signaling. In addition, pathways connected to ECM organization, axon guidance, and developmental biology were also differentially expressed (Table 3). Although three pathways were found in common between the two SAM analyses, only 48% (271/504) of the genes overlapped (data not shown).

To further explore this finding, we performed a gene set analysis of the whole dataset (10,000 genes), revealing genes of the axon guidance pathway to be differentially expressed between OCCC and the other subtypes combined (P < 0.001, 156/251 genes in pathway present in dataset). However, when OCCCs were compared to the other subtypes individually, the genes of the axon guidance pathway were only found to be differentially expressed between the OCCC and the EM subtypes (P < 0.001; File S1 in Supplementary Material).

Comparing the individual differentially expressed pathways between OCCC and the other subtypes, we found the highest number of differentially expressed pathways when comparing OCCC to the EM subtype. A total of 324 pathways, encompassing 4,599 genes, included genes that were differentially expressed between the OCCC and EM subtypes (Figure 2). However, a significant overlap of genes between 295 of these pathways (91%) was observed. Furthermore, several pathways were found to be general and were probably a result of these pathways containing other smaller pathways, resulting in an elevated number of significant pathways. However, 30% (96/324) of these pathways had genes which were differentially expressed between OCCC and the serous and EM subtypes combined. Among these 96 pathways were the “TFAP2 (AP-2) family regulates transcription of growth factors and their receptors” pathway, the “Chromatin organization” pathway, and several PI3-kinase and AKT associated pathways (File S1 and File S2 in Supplementary Material). A comparison of the genes in the differentially expressed pathways revealed three genes of interest: the Rho GTPase CDC42, the growth factor receptor ERBB2, and the transcription factor TFAP2A. These genes were differentially expressed between subtypes, with overexpression in OCCC (Figure 1B), and were also present in many of the altered pathways in both the SAM and whole gene set analyses. Genes in the Rho GTPase pathways were also differentially expressed between OCCC and both the EM and serous subtypes; however, no direct overlap in pathways between subtypes was observed.

The overexpression of ERBB2 in the OCCC subtype prompted us to investigate the protein expression and copy number status of ERBB2 using a TMA. Only one of 43 OCCC tumors (2%) showed overexpression and amplification of ERBB2 (Figure S1 in Supplementary Material), and one tumor displayed amplification (ratio ≥ 2) with weak protein expression (1+). Three tumors (7%) displayed weak protein expression (1+) without amplification. This suggests that while ERBB2 may be overexpressed in OCCC on the mRNA level relative to other histological subtypes of ovarian cancer, this does not result in overexpression on the protein level, nor does it appear to be caused by amplification of the ERBB2 gene. Further support of this was found when applying a module score for the ERBB2 pathway from breast cancer to our data (P = 0.2896, two sided t-test, 17/28 genes present) (43), which did not provide evidence for ERBB2 mediated activity in OCCC (Figure 1C). However, there was a significant difference between the OCCC and serous subtypes (P = 0.0147, Dunn’s test, 17/28 genes present), with a lower module score in the serous subtype, in line with a previous study from our group (42).

In a previous study, we evaluated the protein expression of mTOR, PTEN, and EGFR in OCCCs in cohort 1 (n = 12). EGFR was expressed in two tumors (17%), PTEN in four tumors (33%), and mTOR in 7 tumors (58%) (29). pAKT was evaluated using the previously mentioned TMA of 43 OCCCs. pAKT protein expression was observed in 23 tumors (58%).

Next, we performed a GO term similarity score analysis using the genes from the SAM analysis without variance filter and compared those to a previous study of OCCC (including the clear cell, serous and EM subtypes) (16) using G-SESAME. The similarity score was 0.683, confirming our findings, and GO terms above this threshold were related to anatomical structure development, cell differentiation, and signaling (Table S2 in Supplementary Material).

Taken together, these findings indicate that different mechanisms involved in structural modifications, such as Rho GTPases and chromatin organization genes, may play an important role in tumor development and progression in OCCC. In addition, upregulation of the transcription factor TFAP2A may be involved in the distinct transcriptional differences observed between OCCC and other subtypes. Furthermore, differential activation of the PI3K and mTOR pathways, implicated in the development of type I ovarian carcinomas, between OCCC and both EM and serous tumors, may indicate deregulation of the PI3K/mTOR/AKT pathway in OCCC.

The transcriptional differences observed across the ovarian cancer subtypes, specifically the appearance of a distinct OCCC gene expression profile, prompted further analyses of the potential underlying mutational drivers of OCCC. We therefore performed targeted deep sequencing of 60 cancer-related genes in the 10 OCCCs from cohort 2; >99% of the target bases had at least 30X coverage across all 10 samples, with a mean target coverage of 534X (range 314X-698X). A total of 8,595 variants were called, spanning 3,849 unique variants across the 10 tumors. Of these, 1,244 were tagged by GEMINI as potential germline variants, while Mutect2 filtered out 7,319 variants, leaving 32 variants (0.3%) using a minor allele frequency cutoff of 5%.

The 32 mutations were manually compared to the COSMIC, ExAC, and ENSEMBL (dbSNP148) databases to identify previously reported variants and filter out variants with a global minor allele frequency >0.0001 (i.e., 1/10,000, suspected to be germline SNPs). Twenty probable somatic variants were identified in 11 genes, with a variant range from 0 to 5 (median 2), supporting the notion of a genomically stable subtype. Two variants were inframe indels, 10 were missense variants and eight were truncating variants. Forty-five percent (9/20) of the variants have not been previously reported (Figure 3A and Table 4). Using the gene family nomenclature provided by the online Molecular Signatures Database (MSigDB V5.1; Broad Institute, MA), 2 of these genes are tumor suppressors (CDKN2A, TP53), 6 are oncogenes (ERBB2, PIK3CA, KRAS, ALK, NOTCH1, ROS1), 2 are transcription factors (TP53, SPOP), and 3 are protein kinases (ERBB2, ALK, ROS1). Next, the PantherDB online tool was used to explore potential affected pathways. Several pathways, including the TP53, WNT, EGF, and PIK3CA/KRAS pathways were found to be affected (Figure 3B).

ARID1A and PIK3CA variants, previously identified in OCCC, were present in 4/10 (40%) tumors, and co-occurred in three tumors (30%). Samples 2 and 9 had >1 ARID1A variant each. The six ARID1A variants included three small indels and three missense variants; five of these variants resulted in protein truncation downstream of the variation site (Table 4). Only one variant, the truncating variant Gln2176 in Sample 9, has previously been reported, although not in ovarian cancer. Of the three PIK3CA variants, the missense variant His1047Arg occurred in two samples and has been associated with ovarian cancer in 91 tumors in the COSMIC database. Two variants had previously been reported in COSMIC, while the EDLLNPI453del variant in Sample 1 has not been reported.

The ERBB2 variant (p.989I > M) has not previously been reported in any of the databases used in this study. However, it is located outside the functional domains of the ERBB2 gene, and the significance of this finding is therefore unknown (data not shown).

Two TP53 single nucleotide deletion variants were found, both of which are reported in COSMIC (Table 4). The missense variant Arg181Cys has been identified in breast and endometrial cancer according to COSMIC, but not in ovarian cancer.

The KRAS variant found in one tumor (c.35G > T) has previously been reported as a somatic variant in ovarian cancer (COSM520).

The CDKN2A variant has previously been identified in COSMIC, however not in ovarian cancer. Finally, variants in the genes ROS1, NOTCH1, ALK, KMT2D, and SPOP have to our knowledge not previously been reported in ovarian cancer.

In summary, genes previously reported to be mutated in OCCC, including ARID1A, PIK3CA, and ERBB2, were identified in the present cohort. Variants in samples 3 and 7 were validated using available exome data with paired tumor and blood, supporting the use of FFPE tissue for confident variant calling using strict criteria as implemented in the BCbio pipeline (Table 4). Overall, genes affected by potential somatic variants in OCCCs in the present study included KMT2D, KRAS, ARID1A, PIK3CA, TP53, ERBB2, NOTCH1, and SPOP. These genes, spanning variants over 8/10 tumors, are associated with a comprehensive list of functions that relate to multiple cancer driving cellular processes. However, in common are chromatin remodeling and transcriptional regulation relying on chromatin, providing a possible explanation to a common underlying cause in OCCC pathogenesis.