Multiomics Characterization of Low-Grade Serous Ovarian Carcinoma Identifies Potential Biomarkers of MEK Inhibitor Sensitivity and Therapeutic Vulnerability

Low-grade serous ovarian carcinoma (LGSOC) is a rare subtype of epithelial ovarian cancer. Women with LGSOC are usually diagnosed with advanced-stage disease and only 10% to 20% of such patients survive more than 10 years after diagnosis (1, 2). Research on LGSOC is challenging due to its low prevalence, uncertain etiology, and the limited availability of research models. However, in the last 5 to 10 years, investigators have elucidated many key genomic aberrations, leading to major advancements in the molecular characterization and classification of LGSOC (3, 4). In contrast to other ovarian cancer subtypes, LGSOC contains few point mutations (5), yet are characterized by a high frequency of oncogenic mutations in KRAS, NRAS, and BRAF (20%–40%, 7%–26%, and 5%–33%, respectively; refs. 6–9), a very low prevalence of TP53 mutations (<8%; refs. 10, 11), loss of CDKN2A (15%–53%; refs. 12, 13), and high expression of estrogen receptor (ER) and progesterone receptor (>90% and >50%, respectively; refs. 14–16). More recently, less frequent mutations in USP9X (15%) and EIF1AX (8%) have been described in a small proportion of LGSOC tumors (12, 17).

Historically, the same therapeutic agents have been used for the treatment of the many different histotypes of epithelial ovarian cancer. We now recognize that these histotypes are characterized by unique molecular features that define each of them as different disease entities. Novel treatments that successfully exploit the molecular vulnerabilities in cancers are leading to the development of different drug therapies specific to each molecular subtype (18). Patients with LGSOC are often treated with cytoreductive surgery, followed by platinum/paclitaxel chemotherapy. However, unlike the most common subtype of epithelial ovarian cancer (high-grade serous ovarian carcinoma or HGSOC), LGSOC is relatively resistant to standard chemotherapy, likely due to its low proliferative activity relative to HGSOC (2, 19). Antiestrogen therapy is also a common treatment, although response rates in patients with recurrent disease are low (9%–14%; ref. 20). In this context, the identification of oncogenic mutations affecting MAPK genes (RAS/RAF/MEK) in LGSOC led to the evaluation of targeted agents, such as MEK inhibitors (MEKi), in clinical trials (21–24). In a recently reported randomized trial (NRG-GOG281), patients with recurrent LGSOC were found to have a statistically significant improvement in progression-free survival (PFS) when treated with the MEKi, trametinib, compared with treatment by physician's choice standard of care (24). Trametinib monotherapy resulted in an improvement in PFS of almost 6 months (13 vs. 7.2 months) with objective response rates of 26% (compared with 6% in the standard treatment arm; ref. 24). While promising, these results highlight a clear need for more effective therapies.

Advancements in the treatment of rare cancers, such as LGSOC, rest on a better understanding of their unique molecular features and how they can be exploited for therapeutic gain. Unfortunately, results from genomic precision oncology clinical trials have shown that only 5% to 7% of patients with cancer derive benefit (25). However, in this small group of patients, response rates were more than 50% and durable response durations of 30 months were observed. To accelerate progress in precision oncology, future research must go beyond the genomic aberrations found in tumors and explore other complex aspects of tumor biology. Thus, multiomics tumor profiling with data integration analysis (26) provides a new window of opportunity for the identification of molecular disease drivers leading to more successful drug therapies. This approach can be extremely useful in the study of rare cancers, such as LGSOC, where there is an urgent need to develop more effective therapies.

In the last few years, we have established a large collection of cancer cell line models from patients with advanced/recurrent LGSOC (27). These models have since been used for preclinical drug testing and biomarker discovery in several studies (27–29). With the introduction of MEKi therapy for LGSOC, there is a pressing need to identify biomarkers that predict tumor response. These predictive biomarkers will help us identify those patients who are more likely to benefit from MEKi therapy, and through further research, identify novel treatment targets for the MEKi-resistant patient population (29).

In this study, we performed an integrative analysis of the genomes, transcriptomes, and proteomes of 14 clinically characterized LGSOC cell line models with the goal to identify novel tumor vulnerabilities for the development of future drug treatments. Here, using multiomics approach, we aimed to: (i) characterize the molecular complexity of our LGSOC models, (ii) elucidate novel potential drivers associated with LGSOC disease, and (iii) identify molecular biomarkers that predict response to MEKi therapy.

Advanced or recurrent LGSOC samples (tumor and ascites) were obtained from the gynecologic tumor bank [Vancouver General Hospital (VGH), Vancouver, British Columbia, Canada and BC Cancer (BCC)], and London Health Sciences Centre (LHSC; London, Ontario, Canada). Tumor bank protocols, cell line derivation, and all research relating to this study were approved by Institutional Human Ethics Review Boards at BCC (H14-02859), the University of British Columbia (Vancouver, British Columbia, Canada; R05-0119), and the Western University (London, Ontario, Canada; HSREB 12668E). Clinical information was extracted retrospectively from patient records. At the time of the original diagnosis, all cases were reported by a gynecologic pathologist to confirm LGSOC tumor histology. All cases were also recently rereviewed to confirm the original pathologic diagnosis. Cases were pathologically evaluated according to the 2014 World Health Organization criteria, which rely on histologic criteria (nuclear features and mitotic count) for diagnosis. Images of the study cases are shown in Supplementary Fig. S1A–S1N.

Low-grade serous ovarian cancer patient–derived cell lines were established through continuous in vitro culture of patient material obtained through VGH (Vancouver, British Columbia, Canada) or LHSC (London, Ontario, Canada; cell line CL-01) tumor banks. Cell line development and the confirmation of their identity using short tandem repeat testing are described previously (27, 29). Details of the cell lines, including the passage number used for experimentation, are shown in Supplementary Table S1. As described in our previous publication (29), we used a binary system to classify LGSOC cell lines in response to in vitro MEKi testing. Cell lines were classified as sensitive when they showed complete cell death, at the end of the experiment (120 hours), regardless of the MEKi drug tested. MEKi-resistant lines were defined as resistant when cells survived drug treatment and continued to proliferate. While this binary classification is not conventional for characterizing cell line drug sensitivity, the utility of this classification has been clinically validated. Using this classification, we identified KRAS mutations as a predictive biomarker. The finding was subsequently validated in the MILO/ENGOT-ov11 clinical trial (NCT01849874; ref. 23).

DNA was extracted from LGSOC cell lines using an AllPrep DNA/RNA Mini Kit (Qiagen), and from matched patient buffy coat samples using a DNeasy Blood and Tissue Kit (Qiagen), according to the manufacturer's protocol. DNA concentration was quantified using a Qubit dsDNA HS Assay (Thermo Fisher Scientific). About 0.5 μg of genomic DNA was fragmented by Hydrodynamic Shearing (Covaris, Inc.) to generate 150- to 280-bp fragments. After end repair, fragments were ligated to Illumina barcoded adapters, cleanup, amplified by PCR to enrich for adapter-ligated fragments, and controlled for quality. The library was then enriched by liquid-phase hybridization using Agilent SureSelect XT Human All Exon v6 (Agilent Technologies), following the manufacturer's recommendations. The captured library was amplified by PCR using indexing primers, cleanup, and quality assessment was done with the TapeStation 4200 (Agilent Technologies). Libraries are submitted to PE100 sequencing on Illumina HiSeq4000.

Sequence alignment and mutation calling were performed in Partek Flow Environment (Partek Inc). Sequence reads were aligned to the human genome hg38 build using bwa 0.7.2 (30). Variants were called using Strelka 1.0.15 (31) for all cell lines, and the potential germline variants present in buffy coat samples from matched patients were filtered out, except for the cell line CL-02, which was missing a matched buffy coat sample. CL-02 variant calling was performed using LoFreq 2.1.3.a (32). The called variants were annotated using the Annovar software (33). Annotated calls were then filtered to show only protein-changing single-nucleotide variations (SNV) that were present in cell line DNA at allele frequencies greater than 0.1 and coverage higher than 16 ×. For CL-02, all calls not present in dbSNP (version 138) were retained, while of the calls that were present in dbSNP, calls with (average heterozygosity + aveHet SE) less than 0.1 were retained. These were additionally filtered using the same criteria as for the Strelka calls.

Data analysis was performed using Nexus Copy Number Discovery Edition Version 9.0 (BioDiscovery, Inc.). Samples were processed using the Nexus next-generation sequencing functionality with the FASST2 segmentation. The log ratio thresholds for single-copy gain and single-copy loss were set at +0.18 and −0.18, respectively. The log ratio thresholds for gain of two or more copies and for a homozygous loss were set at +0.6 and −1.0, respectively. Tumor sample binary alignment map (BAM) files were processed with corresponding normal tissue BAM files. Probes were normalized to the median.

Total RNA was isolated using the mirVana Isolation Kit (Thermo Fisher Scientific). Strand-specific RNA sequencing (RNA-seq) was performed on quality controlled high RNA integrity number value (>7) RNA samples (Bioanalyzer Agilent Technologies). In brief, 200 ng of total RNA was first treated to remove the rRNA, and then purified using the Agencourt RNA Clean XP Kit (Beckman Coulter) prior to analysis with the Agilent RNA 6000 Pico Chip to confirm rRNA removal. Next, the rRNA-depleted RNA was fragmented and converted to cDNA. Subsequent steps included end repair, addition of an “A” overhang at the 3′ end, ligation of the indexing-specific adaptor, followed by purification with Agencourt Ampure XP beads. The strand-specific RNA library was prepared using TruSeq Stranded mRNA Kit (Illumina). Libraries were quality checked and sized with a TapeStation 4200 (Agilent Technologies), and then quantified using the KAPA Library Quantification Kit for Illumina Platforms (Kapa Biosystems). Libraries were submitted to PE75 sequencing on Illumina NextSeq500. About 20 to 50 × 10⁶ sequencing reads per library were analyzed for each sample.

RNA-seq alignment and transcript/gene quantification were performed in Partek Flow Environment (Partek Inc). The RNA-seq reads were first trimmed to remove low quality bases from the 3′ end. The trimmed reads were then mapped to the human reference genome hg38 and transcriptome (RefSeq release 83) using splice-aware aligner STAR (2.5.3a; ref. 34). On the basis of the reads that can only be mapped to a single-genomic location, the transcript/gene expression quantification was performed using HTSeq-count (35). Cross-sample normalization of expression values was done by DESeq (36).

Seven LGSOC cell lines were selected for global proteome analysis (Supplementary Table S1) in triplicates using three separately cultured samples for each. Cells were collected by trypsinization and washed with 1 × Dulbecco's phosphate buffered saline. A minimum of 1 × 10⁶ cells per cell line were used in these analyses. Cell pellets were lysed in 200 μL of a buffer containing guanidine hydrochloride (4 mol/L), HEPES (pH 8.5, 50 mmol/L), 2-chloroacetamide (40 mmol/L), tris(2-carboxyethyl) phosphine (10 mmol/L), ethylenediaminetetraacetic acid (5 mmol/L), 1 × phosphatase inhibitor, and 1 × protease inhibitor. Bead beating of the cells was performed in lysing matrix D tubes on a FastPrep24 [6 meters per second (m/s) for 45 seconds, 60 seconds rest, and 6 m/s for 45 seconds]. Cell lysate was then heated at 95°C for 15 minutes with 1,200 rpm mixing. A total of 50 μg of each sample was diluted 10× in 50 mmol/L HEPES, pH 8.5 (concentration was determined by using BCA assay) and digested using a 50:1 ratio of protein:trypsin/Lys-C mixture at 37°C overnight with shaking at 1,200 rpm. A small aliquot of each LGSOC sample was mixed and analyzed as LGSOC-mixed control. The five HGSOC samples were mixed and used as a HGSOC-mixed control sample. An aliquot of all the LGSOC and HGSOC samples was mixed and used as an all-sample pooled internal standard. Tandem mass tag (TMT10) labeling was performed by adding 100 μg TMT label in acetonitrile and incubating at room temperature for 30 minutes twice. Excess label was quenched with 10 μL of 1 mol/L glycine. Individual TMT channels were combined, and the volume of the combined TMT-labeled sample was reduced to 10% to 20% of the original volume. Three 10-plex batches were prepared, where each batch contained one replicate of each of the LGSOC samples, the LGSOC-mixed control, the HGSOC-mixed control, and the all-sample pooled internal standard. The combined samples were fractionated by high pH reverse-phase high-performance liquid chromatography in a gradient of acetonitrile and 20 mmol/L ammonium bicarbonate in water into 48 individual fractions. These 48 fractions were concatenated into 12 fractions by combining every 12th fraction, each of which was analyzed by low pH nanoLC mass spectrometry (MS) using a MS/MS/MS method on a Thermo Fisher Scientific Easy nLC coupled to a Thermo Fisher Orbitrap Fusion MS. The UniProt Human Proteome was used as a reference proteome. Only peptides that were identified in all three batches were retained. ComBat, a function from sva R package (37), was used to adjust for batch effects using an empirical Bayes adjustment. Detection of differentially expressed proteins from the collected peptide data was determined by the probe-level expression change averaging in R (38). No enrichment for specific posttranslational modifications was performed.

We used deconstructSigs (39), a multiple regression approach to statistically quantify the contribution of mutational signature for each tumor. The 30 mutational signatures were obtained from the Catalogue of Somatic Mutations in Cancer (COSMIC) mutational signature database (40). In brief, deconstructSigs attempts to recreate the mutational pattern using the trinucleotide mutation context from the input sample that closely resembles each of the 30 mutational signatures from COSMIC mutational signature database (v2, March 2015; ref. 40). In this process, each mutational signature is assigned a weight normalized between 0 and 1 indicating its contribution. Only those mutational signatures with a weight more than 0.06 were considered for analysis.

HIT'nDRIVE (41) measures the potential impact of genomic aberrations on changes in the global expression of other genes/proteins that are in close proximity in a gene/protein interaction network. It then prioritizes those aberrations with the highest impact as cancer driver genes. Both nonsilent somatic mutation calls and copy-number aberration (CNA) gain or loss were independently collapsed in a gene–patient alteration matrix with binary labels. mRNA and protein expression values were used to derive expression-outlier gene–patient outlier matrix using GESD test. STRING ver10 protein interaction network was used to compute pairwise influence value between the nodes in the interaction network. We integrated these genome and transcriptome, as well as genome and proteome data using HIT'nDRIVE. First, we integrated genome and transcriptome data. Here, we ran HIT'nDRIVE separately using SNV-mRNA expression data and CNA-mRNA expression data. For this, we used the parameters: α = 0.9, β = 0.7, and γ = 0.7. Next, we integrated genome and proteome data. Here, we ran HIT'nDRIVE separately using SNV-protein expression data and CNA-protein expression data. For this, we used the parameters: α = 0.9, β = 0.7, and γ = 0.8. We used IBM-CPLEX as the integer linear programming solver. The results were later combined to downstream analysis.

Differential expression analysis of MEKi response phenotypes was performed by applying linear empirical Bayes model using sva R package (37). Gene (mRNA) and protein expression data of the MEKi response phenotypes were used for this purpose. We used the following threshold values to select genes and proteins for MEKi response analysis and downstream pathway analysis: (i) for mRNA expression data, we included variables in the analysis if statistical testing resulted in both P ≤ 0.05 and an expression fold change difference of greater than 1.03 or (ii) for protein expression data, we included variables in the analysis if testing resulted in both P ≤ 0.5 and a fold change of greater than 1.5. Furthermore, to identify differentially expressed protein complexes, we used Wilcoxon rank-sum test on the average protein expression profiles of protein complex members. Protein complexes were selected for protein complex analysis if testing resulted in P ≤ 0.05 with a fold change of greater than 1.

The differentially expressed set of genes or proteins was tested for enrichment against gene sets of Kyoto Encyclopedia of Genes and Genomes pathways present in Molecular Signature Database v6.0. A hypergeometric test–based overrepresentation analysis was used for this purpose (https://github.com/raunakms/GSEAFisher). A cutoff threshold of FDR ≤ 0.01 was used to obtain the significantly enriched pathways. Only pathways that were enriched with at least four differentially expressed genes were considered for further analysis. To calculate the pathway activity score, the expression dataset was transformed into standard normal distribution using the “inverse normal transformation” method. This step was necessary for a fair comparison between the expression values of different genes. For each sample, the pathway activity score was the mean expression level of the differentially expressed genes linked to the enriched pathway.

For a given protein complex, its coexpression score was computed as the average Pearson correlation coefficient of all pairwise protein–protein interactions (as represented in STRING ver 10 database) in the complex. The coexpression scores were separately calculated using mRNA expression data (⁠|{R_{mRNA}}$|⁠) and protein expression data (⁠|{R_{protein}}$|⁠). Protein complexes with (⁠|| {{R_{mrna}} - {R_{protein}}} | \le 0.05$|⁠) were defined as correlated protein complexes. For this analysis, we considered the comprehensive resource of mammalian protein complexes (CORUM) protein complexes (42) consisting of at least four proteins.

We utilized DNA sequencing datasets of publicly available patient LGSOC cohorts from the American Association for Cancer Research (AACR) Project Genomics Evidence Neoplasia Information Exchange (GENIE; ref. 43). The dataset consisted of a total of 122 LGSOC tumors (both primary and metastatic tumors). We used somatic mutation and CNA profiles from the dataset. AACR Project GENIE data: version 5.0 was downloaded from (SynapseID: syn7222066).

The whole-exome and whole-transcriptome sequencing data from this study are available in the European Genome-phenome Archive, under accession number EGAS00001004724. The proteome data from MS proteome data are available in the PRIDE Archive, under accession number PXD019544.

To investigate the landscape of somatic mutation changes in LGSOC cell lines, we performed high-coverage whole-exome sequencing of 14 LGSOC cell lines and 10 matched normal samples from 11 independent patients. In 3 patients, two distinct cell lines were established from tumor/ascites samples taken at different timepoints during the patients' treatment course (Supplementary Table S1). We refer to this cohort of LGSOC cell lines as the VGH cohort. We achieved a mean coverage of 106 × for normal samples and 107 × for tumor samples (Supplementary Table S2). Altogether, we detected 1,893 highly confident nonsilent mutations in the protein coding regions in 14 cell lines (Supplementary Tables S3 and S4).

As described in our recent study (29), nonsynonymous mutations in KRAS were found in all four MEKi-sensitive LGSOC cell lines (CL-01, CL-02, CL-14, and CL-15; all derived from independent patients with LGSOC; Fig. 1A). Two of them (CL-01 and CL-14) harbored KRAS^G12D mutations, while the other two (CL-02 and CL-15) harbored KRAS^G12V mutations (Fig. 1B and C). The variant allele frequency (VAF) of these KRAS mutations was more than 70% in CL-14, CL-15, and CL-01 and that in CL-02 was 55% (Supplementary Fig. S2), indicating that these mutations are highly likely to be clonal in origin. Interestingly, all mutations detected in KRAS were at the glycine residue at position 12. NRAS was mutated in one MEKi-sensitive cell line and three MEKi-resistant cell lines. NRAS^Q61R mutation was specific to the three MEKi-resistant cell lines (CL-03, CL-04, and CL-10; CL-03 and CL-04 being derived from the same patient), whereas NRAS^C118Y mutation was present in one MEKi-sensitive cell line (CL-14) with a cooccurring KRAS mutation. The VAF of NRAS mutations ranged between 48% and 62%, also indicating that these mutations are highly likely to be clonal (Supplementary Fig. S2). A BRAF^D594G mutation was present in only one cell line (CL-09; Fig. 1B). We further evaluated somatic mutations in a large cohort of LGSOC tumors (n = 97) collected as part of the AACR Project GENIE (43). This cohort includes patients with both primary and metastatic LGSOC, however, information on tumor stage was not available. The most common mutations identified in the GENIE cohort were KRAS (40%), NRAS (11%), and BRAF (12%; Supplementary Fig. S3A and S3B). Figure 1C shows the specific mutations identified in the two cohorts. The NRAS^C118Y and BRAF^D594G mutations identified in the VGH cohort were not found in the GENIE cohort (Fig. 1C). Ten cell lines, including all MEKi-sensitive cases, had at least one mutated gene in the MAPK pathway (Supplementary Fig. S4A). In the GENIE cohort, RAS/RAF mutations were prevalent (63%) and were mutually exclusive (Supplementary Fig. S3A).

Apart from mutations in the MAPK signaling pathway axis, other sparsely mutated genes identified in the LGSOC cell lines were involved in Notch pathway, chromatin remodeling, and DNA repair (Fig. 1A; Supplementary Fig. S4A). Notably, TP53^R234H mutations were observed in two MEKi-resistant cell lines (CL-07 and CL-08) both with VAF above 96%, indicating their clonal origin (Fig 1A; Supplementary Fig. S2). While mutations in TP53 are rare in LGSOC tumors (8% frequency; refs. 10, 11), we derived two cell lines from sequential tumor samples obtained from the same patient with a BRCA1 mutation. This unusual case was described in a previous article (44), and was confirmed on pathology rereview. It is also noteworthy that 3% of LGSOCs in the GENIE cohort had TP53 mutations. Notch pathway genes were mutated in five cell lines; mostly in MEKi-resistant cases. Signaling pathways, such as PI3K, chromatin remodeling, MYC, TP53, and NRF2, were mostly mutated in MEKi-resistant cell lines (Supplementary Fig. S4A). Similarly, different subsets of LGSOC tumors in GENIE cohort also harbored mutations in DNA repair pathway, chromatin remodeling pathway, Notch pathway, and TP53 pathway (Supplementary Fig. S4B).

Using the software deconstructSigs (39), we evaluated the characteristic mutation patterns in our patient-derived cell lines against the mutational signatures obtained from the COSMIC mutational signature database (40). Analyses of the mutation patterns in the cell lines revealed strong enrichment of C>T transition nucleotide substitution mutations in MEKi-sensitive cell lines, whereas MEKi-resistant and MEKi-untested cell lines were enriched in both C>T (transition), as well as C>A (transversion) nucleotide substitution mutations (Fig. 1A; Supplementary Fig. S5). Moreover, GENIE cohort of LGSOC tumors also demonstrated strong enrichment of C>T and C>A nucleotide substitution mutations. We found mutational signature 1, predominantly related to age-related mutagenesis, to be operative in MEKi-sensitive cell lines (Fig. 1A). On the other hand, mutational signatures 4, 8, 24, and 29 that are associated with transcriptional strand bias for C>A mutations were operative in MEKi-resistant and MEKi-untested cell lines. Interestingly, mutational signatures 3, 6, and 15 that are associated with DNA mismatch repair were also operative in five MEKi-resistant and MEKi-untested cell lines. By contingency testing, signature 29 significantly discriminated between MEKi-sensitive and MEKi-resistant cell lines (P < 0.01; Fisher exact test).

Clinical annotation of these samples revealed that MEKi-sensitive cell lines were associated with primary tumors (rather than recurrent tumors) that had received a lower number of therapy regimens prior to collection than the MEKi-resistant associated cultures. Interestingly, MEKi-sensitive cell lines, despite being associated with older patients (average 65.5 vs. 51.7 years), their estimated overall survival rates were found to be longer (average 6.5 vs. 4.4 years; Supplementary Table S1). The presence of serous borderline or micropapillary tumor patterns, as well as the time to recurrence were not found to be associated with the two distinct MEKi response phenotypes observed in our LGSOC cell lines.

Finally, using the data available from the GENIE cohort, we compared the exonic tumor mutation burden (TMB) in the four major subtypes of ovarian cancer (HGSOC, LGSOC, clear-cell ovarian cancer, and endometrioid ovarian cancer) in the GENIE cohort. Compared with the other subtypes, LGSOC tumors had a low TMB (Fig. 1D).

CNAs were determined using the whole-exome sequence reads. We identified 6,480 CNAs across the 14 LGSOC cell lines analyzed (Supplementary Table S5). Copy-number changes for all cell lines are shown in Fig. 2A, and those for individual cell lines in Fig. 2B. As described previously, loss of chromosome 9p and gains in chromosomes 8, 12, and 20 are frequent (27). Overall, homozygous copy loss of MTAP was observed in 13 of 14 LGSOC cell lines and its LOH was observed in one cell line (Fig. 2C). Similarly, homozygous copy loss of CDKN2A/B was observed in 12 of 14 LGSOC cell lines and its LOH was observed in two cell lines. These copy-number changes were associated with loss of p16 expression as assessed by Western blotting (data available upon request).

Overall GISTIC analysis of focal amplifications and deletion revealed several recurrent events containing known oncogenic drivers (Supplementary Table S6). These include amplifications of PDGFRB (5q32), FGFR4 (5q35), MYC (8q24), and CCND2 (12p13) along with deletions of CDKN2A (9p21), MAPK1 (22q11), and SMARCB1 (22q11). The four MEKi-sensitive LGSOC cell lines, all containing KRAS mutations, also harbored copy-number gains affecting KRAS. CNAs were identified in different genes targeting oncogenic pathways, such as MAPK, NOTCH, PI3K, Hippo, Wnt, TGFβ signaling pathways, chromatin remodeling, DNA repair pathways, and cell-cycle process (Supplementary Fig. S6). In MEKi-resistant cell lines, we observed copy-number deletions of genes in chromatin remodeling, NOTCH, and Wnt signaling pathways. Interestingly, MEKi-sensitive cell lines were, for the most part, copy neutral on chromosomes 1, 9p, and 10, and tended to have chromosome 5, 8, and 12 copy-number gains.

Using our recently developed computational algorithm HIT'nDRIVE (41), we identified driver genes in these LGSOC cell lines. We first ran HIT'nDRIVE using SNV-mRNA expression data and CNA-mRNA expression data. HIT'nDRIVE analysis prioritized 17 unique driver genes in eight LGSOC cell lines for which matched genome and transcriptome data were available (Fig. 3A; Supplementary Table S7). Similarly, we also ran HIT'nDRIVE using SNV-protein expression data and CNA-protein expression data. This analysis identified 19 unique driver genes in seven LGSOC cell lines for which matched genome and proteome data were available (Fig. 3B; Supplementary Table S7).

Both HIT'nDRIVE analyses using mRNA and protein expression data identified genes in chromosome 9p21, MTAP, CDKN2A, and CDKN2B, as the most common driver genes among all LGSOC cell lines (Fig. 3A and B). For both CDKN2A and MTAP, the mRNA and protein expression profiles (available for seven cell lines) were found to correlate and to be expressed at very low levels. Analysis of CNA in the LGSOC tumors in the GENIE cohort revealed that loss of CDKN2A/B and MTAP was present in about 31.5% (17/54) of all tumors, with few homozygous deletions 7.4% (4/54; Supplementary Table S8).

Using HIT'nDRIVE CNA-mRNA expression analysis, copy-number gains that were shared by more than one cell line included ABCC5, XRN1, BRDT, and NFE2. Similarly, HIT'nDRIVE CNA-protein expression analysis revealed CNA gain of DIDO1, TPD52, DDX51, and SMUG1 in more than one cell line. Interestingly, NEF2, DDX51, and SMUG1 are located on chromosome 12q, indicating the potential role of this region in LGSOC progression. Furthermore, HIT'nDRIVE CNA-protein expression analysis prioritized mutation and loss of gene EIF4G3, specific to MEKi-resistant phenotype.

Next, we performed total RNA-seq (Illumina NextSeq500) on eight LGSOC cell lines with known sensitivity to MEKi treatment (27). We detected 14.3 × 10⁶ reads on an average and more than 27,700 uniquely expressed genes across all samples. Principal components analysis (PCA) using the RNA-seq profiles obtained from each sample revealed that, except for the sample CL-01, all other MEKi-sensitive samples were clustered together, whereas MEKi-resistant samples had high variance and thus, were spread along the principal component axes (Fig. 4A). This suggests that the MEKi-sensitive samples are transcriptionally very similar to each other, whereas the MEKi-resistant samples are very diverse.

Next, we performed proteome analysis by MS on seven LGSOC cell lines in triplicate and identified 5,110 unique proteins. The PCA on the MS protein expression data also revealed two distinct clusters for MEKi-sensitive and MEKi-resistant samples (Fig. 4B). This demonstrates that the protein expression profiles of LGSOC cell lines discriminate the MEKi-sensitive and MEKi-resistant phenotypes better than their mRNA expression profiles. The CL-01 sample was distinct from other MEKi-sensitive samples when we used either protein or mRNA data for the PCA. It should be noted that the protein samples were run in triplicate, although the results for individual samples were highly correlated (Supplementary Fig. S7). To ensure that the findings were not influenced by running these samples in triplicate, the PCA was repeated using individual samples and the results were the same.

We then identified differentially expressed transcripts and proteins between MEKi-sensitive and MEKi-resistant phenotypes (see Materials and Methods). A total of 1,231 and 1,202 differentially expressed genes (in mRNA expression profiles) and proteins (in protein expression profiles) were respectively identified, of which, 113 gene/protein pairs were detected in both datasets (Supplementary Fig. S8A; Supplementary Tables S9 and S10). Similar to the results obtained from the PCA, the expression profiles of the differentially expressed proteins discriminated the MEKi-sensitive and MEKi-resistant phenotypes better than that of the differentially expressed genes (mRNA; Supplementary Fig. S8B and S8C). The differentially expressed genes/proteins between MEKi-sensitive and MEKi-resistant LGSOC cell lines included EGFR, PRKCA, MAPK4, NF2, MTOR, MTAP, ATM, TGFB2, and BRCA2.

To identify signaling pathways dysregulated by the differentially expressed genes/proteins between MEKi-sensitive and MEKi-resistant cell lines, we performed overrepresentation analysis (see Materials and Methods; Supplementary Tables S11 and S12). Although, distinct sets of genes (in the mRNA expression profile) and proteins (in the protein expression profile) were found to be differentially expressed between the MEKi response phenotypes, intriguingly, we observed similar sets of differentially expressed signaling pathways enriched between the MEKi response phenotypes between the two datasets (Fig. 4C). Majority of the proteins in the signaling pathways, such as MAPK pathway, Wnt pathway, P53 pathway, DNA replication, DNA repair pathway, cell-cycle process, and amino acid metabolism, were found to be highly expressed in MEKi-resistant as compared with MEKi-sensitive cell lines. On the other hand, majority of the proteins in oxidative phosphorylation and tight junction pathways were expressed in the MEKi-sensitive lines. Similar to the observation made in Fig. 4A and B (Supplementary Figs. S8B and S8C and S9A–S9C), the pathway level analysis also indicated that the protein expression profiles separated the MEKi response phenotypes much more distinctly than that from mRNA expression data. Thus, subsequent analyses focused on differences between MEKi response phenotypes based on protein expression data.

As a novel discovery approach, to identify large protein complexes characterizing the differences in the MEKi response phenotypes, we leveraged a curated set of core protein complexes from the CORUM database (42). To identify protein complexes that characterize the MEKi response phenotypes, we performed Wilcoxon rank-sum test on the average protein expression profiles of CORUM protein complex members (see Materials and Methods). In MEKi-resistant cell lines, we identified protein complexes such as BRCC complex, p130–ER-alpha complex, VEGFR2 complex, Condensin complex, minichromosome maintenance (MCM) complex, and replication factor C (RFC) complex had high protein expression of their complex members (Fig. 5A and B; Supplementary Table S13). In MEKi-sensitive cell lines, we found F1F0-ATP synthase complex and RNA-induced silencing complex had high protein expression of their complex members.

Furthermore, to identify protein complexes with all of their members coexpressed at both transcript and protein levels, we calculated coexpression score of each protein complex (see Materials and Methods). We found 56 (of 343) protein complexes with correlated mRNA and protein coexpression (Fig. 5C; Supplementary Table S14). We also found 28 protein complexes with high mRNA coexpression score (⁠|{R_{mRNA}} > 0.5$|⁠), as well as high protein coexpression score (⁠|{R_{protein}} > 0.5$|⁠). This included the Condensin complex, which was found to have strong coexpression scores. MCM and RFC complexes, on the other hand, had strong protein coexpression, but weak mRNA coexpression.

With the knowledge that MEKi treatment shows activity in some patients with incurable LGSOC disease, we further investigated MAPK pathway activation in more detail. First, we correlated mRNA and protein expression fold change of genes/proteins involved in this pathway (Fig. 6A). In MEKi-resistant cell lines, genes such as PRKCA, HSPA2, TGFB2, FGF3, TP53, and FLNC were highly upregulated in both mRNA and protein expression data, and that these findings were highly correlated (Fig. 6B). On the other hand, MAPK13, HSPB1, and CACNA2D1 were upregulated in MEKi-sensitive cell lines. Furthermore, we found that the RAS–RAF–MEK–ERK and PI3K–AKT–mTOR pathways were highly altered in the MEKi-resistant cell lines as compared with the MEKi-sensitive cell lines (Supplementary Figs. S10 and S11).

Fatality rates in women with LGSOC are high due to a lack of effective treatments. Our study is the first to use a multiomics approach to more comprehensively characterize the genome, transcriptome, and proteome profiles of 14 LGSOC patient–derived cell lines. These cell line models closely recapitulate the genomic abnormalities depicted in LGSOC tumors from the GENIE cohort (43) and other previous studies, and highlight common somatic mutations known to involve the MAPK (RAS/RAF/MEK) pathway (6–8). LGSOC cell lines were enriched with C>T and C>A nucleotide substitutions and KRAS mutations. CNAs were characteristic of those described in previous reports of LGSOC (12, 13). As described previously by Jones and colleagues (5), both LGSOC cell lines and tumors were noted to have lower TMB than other ovarian cancer subtypes. Our findings also support earlier observations that BRAF (45, 46) and NRAS mutations are less frequent than KRAS mutations in advanced/recurrent LGSOC tumors, and NRAS mutations are less frequent than described previously (47).

Mutational signatures associated with LGSOC have not been utilized previously for drug prediction. Although these signatures have provided a remarkable basis for characterizing patterns of genomic disruption associated with different cancer etiologies (40), they may also serve as predictive biomarkers. Both LGSOC tumors and cell lines were found to be associated with transcriptional bias for C>A transversion nucleotide substitution mutations in LGSOC, indicating guanine damage that is being repaired by transcription-coupled nucleotide excision repair. Interestingly, mutational signatures also distinguished MEKi response phenotypes in LGSOC cell lines. In particular, mutational signatures 1, 5, and 12 were associated with MEKi-sensitive cell lines, while signatures 8, 4, and 29 were much more prevalent in MEKi-resistant cell lines. In particular, signature 29 significantly discriminated between the MEKi-responsive phenotypes. Mutational signature 1, which is associated with age, correlated with older patients. In the literature published to date, associations between age and response to MEKi have not been described previously. On the basis of our findings, the hypothesis that age is related to response should be evaluated in patient samples.

Loss of the chromosome 9p21 tumor suppressor locus (including CDKN2A and MTAP) is the most common deletion event across many human cancer types (48, 49). Both LGSOC cell lines and tumors had deletions affecting these genes. However, we noted an important difference in the frequency of these events. As described previously, LGSOC cell lines (86%; 12/14) had focal/homozygous deletions in chromosome 9p21, confirmed by a low CDKN2A/p16 expression. In contrast, LGSOC tumors from the GENIE project harbored heterozygous deletions of 9p21 in only 31.5% to 53% of cases. By way of explanation, there may well be a selective advantage for establishing primary cultures/cell lines of LGSOC with 9p loss, through other factors, such as the high frequency of samples from patients with recurrent disease play a role. Interestingly, absence of p16 expression had been associated with an unfavorable prognosis in LGSOC, suggesting a potential role for the use of cell-cycle inhibitors, such as those targeting CDK4/6 (13).

MTAP is a key enzyme in the methionine salvage pathway and is frequently codeleted with CDKN2A because of its chromosomal proximity (50). From large-scale short hairpin RNA–mediated screens across multiple human cancer cell lines, it is known that viability of MTAP-deficient cancer cells is impaired by depletion of the protein arginine methyltransferase PRMT5 (50). Inhibitors of PRMT5 may be a synthetically lethal strategy for this dysregulated metabolic state and are of interest in MTAP- and CDKN2A/B-deleted tumors. The potential therapeutic value of MTAP and CDKN2A deficiency in LGSOC should be further explored.

Our study describes several different techniques to evaluate multiomics data for the identification of disease drivers and MEKi response prediction. Recently, Schutte and colleagues showed integration of genomic and transcriptomic data significantly outperformed less comprehensive analysis in the identification of clinically relevant molecular alterations (51). The combination of DNA- and RNA-based analyses has shown to serve as mutual controls for verifying potential findings (26). As LGSOC is a rare cancer, model systems and patient samples are limited. Thus, efficient analytic tools and the use of biologically complimentary datasets may help reduce FDRs and/or improve prediction accuracy. HIT'nDRIVE analysis exploits the power of analyzing multiomics data and can also be used to enhance discovery in rare cancers, as we demonstrated previously (26). Using HIT'nDRIVE analysis, we confirmed potential driver genes (KRAS, MTAP, and CDKN2A/B) in our LGSOC cell lines. In the recently reported MILO trial (NCT01849874), the presence of a KRAS mutation appears to be a predictive biomarker (23). In a subgroup analysis, patients with advanced or recurrent LGSOC treated with binimetinib had a significantly better PFS (17.5 vs. 10.8 months; P < 0.01) compared with physician choice of chemotherapy. Furthermore, in a recent genotype-matched trial containing 13 patients with LGSOC, a response rate of 30% was reported for patients treated with MEKi-based–targeted combinations (52). While these results are promising, biomarkers that predict treatment efficacy in addition to KRAS are very much needed.

Here, using LGSOC cell lines, we demonstrated that proteomic profiles were better at discriminating distinct MEKi-response phenotypes than transcriptomic profiles, and found that protein alterations affecting TP53, cell cycle, and DNA repair pathways are common in MEKi-resistant cell lines. These pathways should be further studied as potential targets for therapeutic intervention. The addition of proteomic-based analysis facilitates the understanding of the biological implications of genomic change and can be used as discovery tools for both drug prediction and identification of novel therapeutic targets. We found considerable consistency in the identification of differential signaling pathway activity between MEKi response phenotypes using proteomic and transcriptomic data, strongly suggesting that these pathways play a role in MEKi resistance.

MEKi will soon be commonly prescribed as a treatment for women with recurrent LGSOC (53). To identify robust biomarkers that predict MEKi response, we used transcript/protein expression to evaluate MAPK pathway regulation. MAPK signaling feedback loops play a key role in MEKi resistance (54). In MEKi-resistant cell lines, both PI3K-AKT-mTOR and MAPK pathway activity was found to be altered. Recognizing that KRAS mutations may infer MEKi sensitivity, it may seem counterintuitive that MAPK pathway upregulation results in MEKi resistance. However, complex feedback resistance mechanisms have been described within the pathway and RAS mutations may uniquely impact pathway function and cellular signaling dependencies (oncogene addiction) that promote drug efficacy (55, 56). The most differentially expressed protein identified in the MAPK pathway analysis was protein kinase C alpha (PRKCA or PKCα). We identified previously both PKCα and EGFR as biomarkers of MEKi resistance using reverse phase protein array and EGFR was validated in vitro as a therapeutic target (29). Identifying the same two biomarkers using MS represents a form of orthogonal validation, thus lending support for the utility of this multiomics analytic approach for predictive biomarker discovery. The predictive biomarkers described in this study will hopefully guide future research for those groups/investigators who have access to clinical trial samples or additional cell line models to validate these findings.

While MEKi-sensitive LGSOC cell lines depend on MAPK pathway activation for survival (as shown by the lethal effects of MEKi treatment in vitro), future study of the intricacies of MAPK signaling should be explored as a means of improving treatment efficacy in both MEKi phenotypes. More accurate characterization of MAPK signaling, including a better understanding of pathway feedback loops, will be of importance. Therapeutic interventions that target multiple proteins in the pathway (combinatorial therapies) are being evaluated and alternatively, the simultaneous targeting of other pathways might be efficacious.

Finally, we used our proteomic data to identify protein complexes with potential therapeutic relevance in LGSOC. Using the CORUM database, we elucidated novel protein complexes differentially expressed between the MEKi response phenotypes. Previous studies have shown that some of these protein complexes regulate cell proliferation; including RFC in ovarian cancer (57) and MCM in other cancers (58). The MCM complex forms the core of the DNA replicative helicase and, interestingly, is associated with cisplatin resistance in ovarian cancer (59). Increased RFC expression has been shown to be associated with epithelial–mesenchymal transition and poor outcomes in both lung and breast cancer (60). Other complexes we identified, such as the BRCC and the Condensin complexes, have been linked to DNA repair in cancer (61, 62). Although MEKi represents a promising new therapy in LGSOC, MEKi treatment is only associated with a 6- to 8-month improvement in PFS. Increased expression of proteins within the complexes found in MEKi-resistant cell lines provides a strong rationale to study these complexes as potential therapeutic targets.

LGSOC is a rare and often lethal cancer when diagnosed at an advanced stage. Research models have been challenging to develop. As a result, our study is limited in the number of samples available for analysis. Data are not readily available to validate our findings in patients treated with MEKi. In this context, we have used multiomics analyses as a discovery tool to uncover potential driver and drug prediction biomarkers. Our cell line cohort had a higher frequency of p16 loss (chromosome 9p loss) than what is described in LGSOC tumors and thus, models that preserve p16 expression will be needed for further evaluation.

Using a multiomics approach, we have characterized the largest worldwide collection of LGSOC cell lines to improve their use as research model systems. These cell lines accurately reflect the molecular makeup of LGSOC tumor samples characterized to date, although LGSOC cell lines without 9p loss are not adequately represented in our sample. We have also described key molecular aberrations in our cell lines (including signaling pathway aberrations) that describe potential biomarkers of MEKi response and disease drivers. We found that proteomic data were more robust in differentiating MEKi drug sensitivity and identified several protein complexes that warrant further study in LGSOC. Multiomics analytic approaches have the potential to accelerate the discovery of clinically useful biomarkers and novel drug targets so we can leverage this knowledge to improve patient outcomes.

M. Llaurado Fernandez reported grants from Janet D. Cottrelle Foundation, Canadian Institutes of Health Research, Women's Health Research Institute, Canadian Foundation for Innovation, and Terry Fox Research Institute during the conduct of the study, as well as BC Cancer Foundation philanthropic support and Society of Gynecologic Oncology of Canada nonfinancial support. A. Dawson reports grants from Janet D. Cottrelle Foundation, Canadian Institutes of Health Research, Women's Health Research Institute, Canadian Foundation for Innovation, and Terry Fox Research Institute during the conduct of the study. G.B. Morin reports grants from Terry Fox Research Institute and British Columbia Cancer Foundation outside the submitted work. M.S. Carey reports grants from Janet D. Cottrelle Foundation, Canadian Institutes of Health Research, Women's Health Research Institute, Canadian Foundation for Innovation, and Terry Fox Research Institute during the conduct of the study; other funding from Bausch Health Companies Inc., Canopy Growth Corporation, Illumina Inc., and Exact Sciences Corporation outside the submitted work; and BC Cancer Foundation philanthropic support of Society of Gynecologic Oncology of Canada for nonfinancial support. No disclosures were reported by the other authors.

R. Shrestha: Conceptualization, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. M. Llaurado Fernandez: Conceptualization, resources, data curation, validation, investigation, methodology, writing–original draft, project administration, writing–review and editing. A. Dawson: Resources, data curation, validation, investigation. J. Hoenisch: Resources, data curation, software, formal analysis, methodology. S. Volik: Software, formal analysis, methodology. Y.-Y. Lin: Data curation, software. S. Anderson: Data curation, software, formal analysis. H. Kim: Resources, data curation, validation. A.M. Haegert: Resources, data curation, validation. S. Colborne: Resources, validation. N.K.Y. Wong: Data curation, investigation. B. McConeghy: Resources, data curation, software, validation, investigation. R.H. Bell: Data curation, software, formal analysis, investigation. S. Brahmbhatt: Data curation, software, validation. C.-H. Lee: Data curation, validation, investigation. G.E. DiMattia: Resources. S. Le Bihan: Resources, supervision, methodology, project administration. G.B. Morin: Supervision, methodology, writing–review and editing. C.C. Collins: Conceptualization, supervision, funding acquisition, project administration. M.S. Carey: Conceptualization, resources, supervision, funding acquisition, investigation, methodology, writing–original draft, project administration, writing–review and editing.

The results published here are, in part, based upon data generated by the AACR Project GENIE. This study was funded by the Canada Foundation for Innovation (CFI-IF 33440), Terry Fox Research Institute (TFRI NF PPG Project #1062; to C.C. Collins), British Columbia Cancer Foundation (to M.S. Carey), and the OvCaRe research program. The authors extend a special thanks to the MacKenzie, Lawler, Luther, MacRae, Ho, and Ludeman families; the Janet D. Cottrelle Foundation; and to all patients, families, and donors who supported this research. The authors also thank the Society of Gynecologic Oncology of Canada for support in advancing research and knowledge translation in LGSOC, the Canadian Institutes of Health Research for planning and disseminating grant support (NRF 152680), and the Women's Health Research Institute for Catalyst Grant support. The authors thank the members of the Carey and Collins laboratories for their contributions to this study.

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.