Background

Epithelial ovarian cancer encompasses five primary histological subtypes, with high-grade serous ovarian carcinoma (HGSOC) constituting about 75% of all cases (1). The standard treatment for HGSOC diagnosed at stage IIB and beyond involves a combination of surgery and chemotherapy, primarily using platinum-based compounds and taxanes (2,3). While initial chemotherapy results in tumor response in most patients with HGSOC, there is a very high recurrence rate (4). The addition of poly-ADP ribose polymerase and vascular endothelial growth factor A inhibitors to chemotherapy for subsets of patients currently results in a 5-year disease-specific overall survival (OS) rate of approximately 45% for patients with HGSOC. This rate has hardly improved in the last three decades (5–8). Therefore, new insights into the complex biology underlying HGSOC are urgently needed to develop more effective treatment strategies.

Previous studies using bulk transcriptomes of patients with HGSOC have identified expression-based molecular subtypes. However, these subtypes did not provide insights that have translated into novel drug targets (9–11). A common limitation of such studies is their reliance on bulk transcriptomes, containing both tumor cells and tumor microenvironment (TME) components, thus reflecting the average transcriptional patterns of the combination of all biological processes present in the tumors. This averaging often masks subtle transcriptional patterns pivotal to understanding HGSOC biology, especially when these are overshadowed by dominant patterns from other less relevant (non-)biological processes (12). Consensus-independent component analysis (c-ICA) offers an alternative by decomposing such bulk transcriptomes into statistically independent transcriptional patterns (i.e., transcriptional components; TCs)(13,14). This approach reveals both dominant and subtle patterns and provides a measure of TC activity for each sample (15).

In the present study, our aim was to utilize c-ICA to dissect HGSOC transcriptomes to identify as many TCs associated with patient OS as possible, which could reveal potential novel drug targets.

Methods

See Supplementary methods for the extended methods.

Data acquisition

Raw microarray bulk transcriptomes and clinicopathological details for patients with HGSOC, low-grade serous ovarian cancer (LGSOC), non-serous ovarian cancer, and benign ovarian tissues were sourced from the Gene Expression Omnibus (GEO)(16). We exclusively utilized transcriptomes generated from primary tumor samples. Our analysis was confined to samples on the Affymetrix HG-U133 Plus 2.0 platform (GEO accession identifier: GPL570) and excluded cell line samples. The datasets were pre-processed, and quality controlled as previously described (17). Furthermore, for comprehensive analyses, we incorporated transcriptomes from five distinct resources: the Cancer Cell Line Encyclopedia (CCLE, n = 969), Genomics of Drug Sensitivity in Cancer (GDSC, n = 959), Gene Expression Omnibus (GEO, n = 13,810), and The Cancer Genome Atlas (TCGA, n = 8,150), and spatially resolved transcriptomes from 10xGenomics (16,18–20).

Consensus-independent component analysis (c-ICA)

The bulk transcriptome input data underwent a whitening transformation to optimize the subsequent analysis. Following this pre-processing step, we employed the FastICA algorithm to execute the ICA, facilitating the derivation of estimated sources (ESs). We determined the appropriate number of ESs to extract based on the number of principal components that captured 90% of the total variance within the dataset. To rigorously evaluate the stability of the ESs, we conducted 25 separate ICA runs, each with uniquely randomized initial weight factors. ESs with an absolute Pearson correlation exceeding 0.9 were grouped. By analyzing the outcome of these runs, we derived consensus ESs, denoted as transcriptional components or TCs. This was achieved by computing the mean vector of weights from grouped ESs present in over half of the conducted runs. The gene weights within each TC provide insights into the strength and directionality by which an underlying latent biological process influences the related gene expression levels. Subsequently, these TCs are used to generate a consensus mixing matrix (MM). The coefficients within this MM depict the activity scores of the TCs across the samples.

Survival analysis

To discern the relationship between TC activity and patient OS, a univariate Cox proportional hazards analysis was conducted on a select group of patients with available follow-up data (n = 541, Supplementary Table S1). In addition, a multivariate Cox proportional hazards analysis was carried out, including covariates such as age, stage, debulking status, and tumor grade. This latter analysis was based on a subset of patients with comprehensive clinicopathological data available (n = 373, Supplementary Table S1). We implemented a multivariate permutation framework encompassing 10,000 permutations to mitigate the risk of false discoveries. We established the acceptable false discovery rate (FDR) at 1%, maintaining an 80% confidence level, applicable for both the univariate and multivariate analyses.

Survival tree analysis

We performed a survival tree analysis to delineate groups of patients with HGSOC treated with platinum-based chemotherapy based on distinct transcriptional and clinicopathological attributes. The analysis utilized activities of TCs associated with OS (either from univariate or multivariate survival analysis as mentioned in supplementary methods) in conjunction with relevant clinicopathological factors, such as age, tumor stage, debulking status, and grade, as potential classifiers. We divided patients into two subsets using every plausible cut-off point for each classifier and compared the resulting survival curves employing the log-rank statistic. Consequently, the division was based on the most significant classifier at its optimal cut-off based on the smallest p-value of the log-rank test mentioned above. This divisional process was successfully reiterated on the derived subsets until any of the following stipulated conditions was satisfied: i) the total patient count across both subsets fell below 50, ii) the collective number of uncensored events in both subsets was < 25, or iii) one of the subsets contained < 17 patients. To gauge the stability of our classifiers, we performed 20,000 iterations, randomly selecting 80% of the patient group in each iteration. The significance-based ranks of classifiers in these iterations were correlated with those from the primary survival tree.

Associating the identified transcriptional components with biological processes

To discern the biological processes associated with the TCs, we adopted a multifaceted approach encompassing i) Transcriptional Adaptation to Copy Number Alterations (TACNA) profiling, targeting the identification of TCs that reflect the downstream implications of copy number alterations (CNAs) on gene expression levels(21); ii) Execution of gene set enrichment analysis (GSEA) for each TC, utilizing gene set collections (n = 16) from The Human Phenotype Ontology (The Monarch Initiative), the Mammalian Phenotypes (Mouse Genome Database), and the Molecular Signatures Database (MsigDB)(22,23); iii) The formation of co-functionality networks on the top and bottom genes of each TC, achieved using the GenetICA methodology, accessible via https://www.genetica-network.com(24). For clusters comprising ≥ 5 genes, the enrichment of the predicted functionality was quantified. This served as the foundation for determining the biological process associated with the TC being examined.

Cross-study transcriptional component projection

To determine whether a biological process captured by an identified TC is also active in other cancer types and to investigate if it is more active in tumor cells or in the TME, we collected raw expression profiles from multiple sources: the Cancer Cell Line Encyclopedia (CCLE, n = 969), Genomics of Drug Sensitivity in Cancer (GDSC, n = 959), Gene Expression Omnibus (GEO, n = 13,810), and The Cancer Genome Atlas (TCGA, n = 8,150)(16,22–24). While the CCLE and GDSC datasets comprise cell line profiles across many solid and hematologic malignancies, the GEO and TCGA datasets offer an extensive set of bulk transcriptomes derived from patient samples spanning 27 tumor types. We pre-processed the raw data as previously described (18). Next, we projected the TCs identified via c-ICA onto the cell line expression profiles from CCLE and GDSC and the patient-derived expression profiles from GEO and TCGA. This projection methodology has been described in more detail previously (18). To identify potential variations in the activity scores of the TCs, we compared the activity scores among cell lines and samples derived from patients within all four repositories. We used an absolute activity score threshold of 0.05 for each TC to pinpoint outlier cell lines and patient tumors with heightened activity.

Determination of spatial transcriptomic profiles’ significant activity locations for individual transcriptional components

To further assess whether a biological process captured by an identified TC is more active in tumor cells or in the TME, we collected publicly available spatial resolved transcriptomic profiles of ovarian cancer samples. Eight were sourced from GEO (study ID GSE211956), and three were sourced from the public dataset repository of 10xGenomics (see supplementary methods for details) generated using the 10xGenomics Visium platform. The samples were from patients with HGSOC, serous papillary, and endometrioid ovarian cancer. Activity for each TC across every location within the spatial samples was ascertained through the cross-study projection methodology referred to in the previous method section (21). We incorporated a permutation-driven approach to discern the markedly active areas within the spatial samples for each TC. We derived a null distribution of activities for each TC-location pairing by performing 3,000 permutations and subsequent projections. The p-value of each observed TC activity quantifies the significance of the deviation of the TC activity at a given location from its baseline null distribution. After this, we visualized the z-transformed p-values using a heatmap, followed by obtaining colocalization scores for each combination of TCs in the spatial transcriptomic profiles for each ovarian cancer sample (25). This visualization aided in highlighting the areas with notable activity aligned against the stained representation of the tissue sample.

Results

An integrated data set containing 1,125 bulk transcriptomes from ovarian tissues

We curated 1,193 bulk transcriptomes from the GEO, including patients with HGSOC, low-grade serous ovarian cancer (LGSOC), non-serous ovarian cancer, and benign ovarian tissues (16). Pre-processing, which included removing duplicates and quality checks, culminated in a refined dataset of 1,125 samples (17). These were extracted from 32 distinct studies (Supplementary Table S2) and represented the entire spectrum of ovarian cancer types, stages, and grades, and included 43 samples from non-malignant ovarian tissue (Supplementary Table S1). The ovarian cancer dataset comprised bulk transcriptomes of patients with HGSOC (n = 678), other serous (n = 110), endometrioid (n = 110), and clear-cell ovarian cancer samples (n = 96). Additionally, for 541 patients, comprehensive survival data was available, as well as additional clinicopathologic information, including age, grade, stage, subtype, treatment schedule, and debulking status for 373 patients (Fig. 1).

Figure 1.

Consensus-independent component analysis identifies 374 transcriptional components

c-ICA on the 1,125 bulk transcriptomes revealed 374 independent TCs. Notably, 135 TCs captured the impact of copy number alterations on gene expression levels. Each TC displayed enrichment for at least one gene set from the 16 gene set collections, with an absolute Z-score of more than two. For example, the number of enriched gene sets from the Hallmark gene set collection in an individual TC ranged from zero to 28 enriched gene sets (interquartile range 3 – 7). The median top Z score for Hallmark gene sets was 3.21 (range 1.55 – 37.54, interquartile range 2.6 – 4.25). A comprehensive database, including all TCs and GSEA outcomes, has been made accessible at http://transcriptional-landscape-ovarian.opendatainscience.net.

Using consensus-independent component analysis (c-ICA), we decomposed the 1,125 bulk transcriptomes, revealing 374 independent transcriptional components (TCs). Each TC displayed enrichment for at least one gene set from the 16 gene set collections, with an absolute Z-score of more than two. The activities of 13 TCs were associated with patient overall survival (OS) in a univariate analysis, with an additional TC (TC166) identified in a multivariate analysis accounting for age, stage, debulking, and tumor grade. Combined, these 14 OS-associated TCs were enriched for gene sets associated with diverse biological processes and clinicopathological characteristics, with four TCs capturing the effects of copy number alterations on gene expression levels.

The activities of six transcriptional components are associated with patient overall survival

For a selected subset of 541 patients—including HGSOC, LGSOC, and non-serous ovarian cancer—with available OS information (Supplementary Table S1), 13 TC activities displayed an association with OS univariately (false discovery rate 5%, confidence level 80% in permutation-based multiple testing framework Supplementary Table S3; Fig. 2). For patients with serous ovarian cancer, treated with platinum-based therapy (n = 301, Supplementary Table S1), lower activity of one additional TC (TC166) was associated with worse OS independent of age, stage, debulking, and tumor grade (Supplementary Table S4). Combined, these 14 OS-associated TCs were enriched for gene sets associated with diverse biological processes and clinicopathological characteristics. Four of these TCs captured the downstream effects of CNAs on gene expression levels (Fig. 2, Supplementary Fig. S1-S3). Survival tree analysis identified ten groups of patients with platinum-treated HGSOC based on the activity of six OS-associated TCs and the presence of two clinicopathological characteristics, namely age and stage (Fig. 3, Supplementary Table S5, median robustness statistic of survival tree = 0.52, interquartile range = 0.36 - 0.69). The survival tree demonstrated good classification power (concordance statistic = 0.72, standard error = 0.021). As expected, patients were divided into separate survival groups based on stage (1/2 vs. 3/4) and age (<53.7 vs. ≥53.7 years). The most significant difference in OS was observed between the cohorts with low and high TC121 activity (Supplementary Table S5). Patients with high TC121 tumor activity exhibited the shortest OS, also observed for the subset of patients with advanced-stage HGSOC (Supplementary Fig. S4, Supplementary Table S6). These robust associations with OS for TC121 in these two subsets of patients indicate the relevance of TC121, irrespective of stage.

Figure 2.

Figure 3.

Distinct biological processes show enrichment in the transcriptional components associated with overall survival

Three of the six TCs associated with OS—TC166, TC247, and TC76—captured the effects of CNAs on the expression levels of genes mapping to chromosome regions 13q12-q14, 11q13-q14, and 9p13-p21, respectively (Supplementary Fig. S5, Supplementary Table S7) (6). The higher activity of TC166 was associated with better OS, whereas the higher activities of TC121, TC247, TC250, TC76, and TC146, were associated with worse OS. The top genes from TC166 were enriched for genes involved in replication and apoptosis. The chromosomal region 13q12-q14 linked to the TC166, contains the tumor suppressor genes retinoblastoma 1 (RB1) and Breast Cancer Type 2 Susceptibility Protein (BRCA2). Loss of heterozygosity of this chromosomal region is frequently observed in both sporadic and hereditary serous ovarian cancers (26,27). The top genes from TC247 were enriched for genes involved in proliferation and immune cell activation, TC76 in replication stress, TC250 in extracellular matrix (ECM) interactions, and TC146 in neurotransmitter signaling.

Intriguingly, the top 100 genes in TC121 revealed a co-functional cluster enriched for genes involved in synaptic signaling, with the corresponding proteins reported to localize to the synaptic membrane of neurons (Fig. 4). Among these were pre-synaptic protein neurexin-1 (NRXN1) and its post-synaptic ligand leucine-rich repeat transmembrane protein 2 (LRRTM2), which regulates excitatory synapse formation (top 20 genes are described in Supplementary Table S8, for more details: http://transcriptional-landscape-ovarian.opendatainscience.net) (28,29). Furthermore, this co-functional cluster included neuron-specific synaptic structure proteins, neurofilament light, and medium chain. Moreover, genes encoding for potassium ion channel proteins integral to membrane repolarization during synapse signal transduction carried high weights in TC121. These genes included KCNC1, KCNN2, and KCNIP1(30–32). Several genes in TC121 encoded proteins related to glutamate receptor signaling, including GRIN2C and SLC7A10 (33). In line with this proposed function, high activity of TC121 was observed in neuroblastoma cell lines but not in ovarian or central nervous system cancer cell lines in the GDSC and CCLE datasets (Fig. 5A, Supplementary Fig. S6, S7). In the GEO and TCGA datasets, high activity of TC121 was observed in glioblastoma multiforme and lower-grade glioma but not in ovarian cancer patient samples (Fig. 5B).

Figure 4.

Figure 5.

Distinct spatial profiles of transcriptional components associated with overall survival

Cross-study TC projection onto spatial transcriptomic profiles from 11 ovarian cancer samples revealed that TC121 was highly active in profiles from the tumor region of the 11 ovarian cancer samples (Fig. 6A; Supplementary Fig. S8). Distinct regions with high activity of the copy number TCs (TC166, TC247, and TC76) in the HGSOC sample overlapped with the region containing cancer cells, as expected. TC250, enriched for extracellular matrix interactions, was also active in the stromal region. The strongest inverse colocalization (colocalization score -2.43) was observed between the activity scores of TC146, enriched for neurotransmitter signaling, and TC76, captured the effect of copy number alterations at chromosome 9p13-p21, at the serous ovarian cancer sample (Fig. 6B, Supplementary Table S9).

Figure 6.

Discussion

In this study, we identified 374 TCs, each enriched for gene sets representing various biological processes in HGSOC samples. Six could stratify patients with HGSOC who had received platinum-based treatment into ten distinct OS groups.

The most significant TC in the survival tree analysis, TC121, captured a clinically relevant subtle transcriptional pattern linked to synaptic signaling not previously recognized in HGSOC. In the survival tree, TC121 identified 12% of the HGSOC patients with the shortest OS and, based on spatially resolved transcriptomic analyzed samples, is active in tumor regions. This observation supports the emerging role of neurons and neuronal projections as cancer hallmark-inducing constituents of the TME (34–36).

Further investigation on whether the activity of TC121 originated from tumor cells or in the TME revealed that the TC121 signal is coming from cells within the TME. The high activity of TC121 in low-grade glioma and glioblastoma multiforme patient samples (Fig. 5B) is in agreement with the presence of neurons in large numbers within the TME of gliomas, where they form functional synapses with tumor cells (37,38). Moreover, TC121 activity was lower in non-brain cancers, such as ovarian cancers, which contain fewer neurons and synapses in the TME compared to brain cancers. We expected TC121 activity to be low in the bulk transcriptomes of all cell lines, since they lack TME. TC121 activity in most cell lines, which includes glioblastoma and ovarian cancer cell lines, was indeed low. Neuroblastoma cell lines, however, exhibited high TC121 activity, which is likely due to retained synaptic formation capacity originating from neuroblast cells (39,40). Lastly, TC121’s high activity observed in small, scattered regions within the tumor of spatially resolved transcriptomic ovarian cancer samples also supports TC121’s role in the TME.

TC121’s significant association with OS underscores the potential significance of synaptic signaling in HGSOC biology. Yet, the neuronal subtype and the molecular mechanisms associated with TC121 remain to be elucidated. A study in human ovarian cancer-bearing mice demonstrated that sympathetic innervation in HGSOC involves adrenergic signaling: norepinephrine released by sympathetic neurons binds to beta-adrenergic receptors on the cancer cells (41–43). This binding triggers the tumor cells to release brain-derived neurotrophic factor (BDNF), which enhances cancer innervation via activation of host neurotrophic receptor tyrosine kinase B receptors, thereby establishing a feed-forward loop of sustained signaling. BDNF and the nerve marker neurofilament protein expression were examined in 108 human ovarian tumors (41). This study revealed that increased intratumoral nerve presence strongly correlates with elevated BDNF and norepinephrine levels, advanced tumor stage, and shorter OS in patients with ovarian cancer. Our analysis indicated that BDNF is a prominent gene (with an absolute weight > 3) in 10 TCs but not in TC121, suggesting that TC121 may indicate a distinct process unrelated to BDNF.

The significance of sensory innervation in HGSOC was evidenced by the co-localization of TRPV1, a marker for sensory neurons, and β-III tubulin, a general neuronal marker, in immunofluorescent staining of histological sections from 75 patients (44). Additionally, a murine model study employing neural tracing identified sensory neurons originating from local dorsal root ganglia and jugular–nodose ganglia, with axons extending into the TME (44). A transgenic murine model lacking nociceptors demonstrated that this specific subtype of sensory neurons was involved in tumor progression (45). Another study showed that reducing the release of Calcitonin gene-related peptide from tumor-innervating nociceptors could be a strategy to alleviate this effect of nociceptors by improving anti-tumour immunity of cytotoxic CD8+ T cells in a melanoma model bearing mice (46). This indicates that the signal from TC121 may represent an indirect influence on tumor cells via interactions with immune cells and the promotion of an immune suppressive TME. Furthermore, in cell lines derived from Trp53^−/− Pten^−/− murine HGSOC, the influence of nociceptors was characterized by the release of substance P (SP), their primary neuropeptide. SP is an alternative splicing product of the preprotachykinin A gene (TAC1) and binds to the receptor neurokinin 1 (NK1R), encoded by the TACR1 gene. NK1R expression was confirmed in murine HGSOC cell line, and SP enhanced cellular proliferation in NK1R-positive murine HGSOC cancer cells in vitro (45). Our analysis identified TAC1 and TACR1 as prominent genes in 15 and 2 TCs, respectively, yet not in TC121, and none of these TCs were associated with patient survival. Currently, there are no drugs specifically targeting tumor innervation in (ovarian) cancer (47). Interestingly, the NK1R antagonist aprepitant effectively inhibited the metastasis-promoting effects of neural substance P in human breast and mammary cancer bearing mice (48), demonstrating the feasibility of such an approach. Therefore, it is essential that the mechanisms driving this nerve growth, the specifics of how nerves within the TME interact with ovarian cancer cells, and how they impact the survival of patients with HGSOC are further elucidated.

Altogether, the present study uncovered a clinically relevant TC linked to synaptic signaling not previously identified in HGSOC. This TC may represent a novel cancer cell-extrinsic mechanism within the TME, illustrating how cancer cells and nerve cells interact to promote enhanced proliferation. A deeper understanding of the molecular aspects of tumor innervation could pave the way for novel drug targets for patients with HGSOC.

Data availability

Microarray expression data was collected from three public data repositories: Gene Expression Omnibus with accession number GPL570 (generated with Affymetrix HG-U133 Plus 2.0), CCLE (generated with Affymetrix HG-U133 Plus 2.0, file CCLE_Expression.Arrays_2013-03-18.tar.gz) available at https://portals.broadinstitute.org/ccle/data and GDSC (generated with Affymetrix HG-U219) available at https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3610/. Pre-processed and normalized RNA-seq data was collected from TCGA using the Broad GDAC Firehose portal (https://gdac.broadinstitute.org/). Spatially resolved samples were sourced from 10xGenomics and GEO. The datasets generated during and/or analyzed during the current study are available in the website: http://transcriptional-landscape-ovarian.opendatainscience.net.

Additional information

Funding

This research was supported by a Hanarth Fund grant, the Netherlands (2019N1552 to R.S.N.F).

Author contribution

A. Bhattacharya: Data curation, methodology, data analysis, data interpretation, writing. T. S. Stutvoet: Data curation, methodology, data analysis, data interpretation, writing. M. Perla: Data interpretation, writing. S. Loipfinger: Data analysis, data interpretation, writing. M. Jalving: Data interpretation, writing. A. K. L. Reyners: Data interpretation, writing. P. D. Vermeer: Data interpretation, writing. R. Drapkin: Data interpretation, writing. M. de Bruyn: Data interpretation, writing. E. G. E. de Vries: Data interpretation, writing. S. de Jong: Conceptualization, data interpretation, writing. R. S. N. Fehrmann: Conceptualization, data curation, methodology, data analysis, data interpretation, writing. All were involved in the final decision to submit the manuscript.

Disclosure

The authors have declared no conflicts of interest.