NECTIN4 (PVRL4) as Putative Therapeutic Target for a Specific Subtype of High Grade Serous Ovarian Cancer—An Integrative Multi-Omics Approach

In high grade serous ovarian cancer patients with peritoneal involvement and unfavorable outcome would benefit from targeted therapies. The aim of this study was to find a druggable target against peritoneal metastasis. We constructed a planar—scale free small world—co-association gene expression network and searched for clusters with hub-genes associated to peritoneal spread. Protein expression and impact was validated via immunohistochemistry and correlations of deregulated pathways with comprehensive omics data were used for biological interpretation. A cluster up-regulated in miliary tumors with NECTIN4 as hub-gene was identified and impact on survival validated. High Nectin 4 protein expression was associated with unfavorable survival and (i) reduced expression of HLA genes (mainly MHC I); (ii) with reduced expression of genes from chromosome 22q11/12; (iii) higher BCAM in ascites and in a high-scoring expression cluster; (iv) higher Kallikrein gene and protein expressions; and (v) substantial immunologic differences; locally and systemically; e.g., reduced CD14 positive cells and reduction of different natural killer cell populations. Each three cell lines with high (miliary) or low NECTIN4 expression (non-miliary) were identified. An anti-Nectin 4 antibody with a linked antineoplastic drug–already under clinical investigation–could be a candidate for a targeted therapy in patients with extensive peritoneal involvement.


Introduction
High grade serous ovarian cancer (HGSOC) is a difficult to treat cancer entity, usually diagnosed at an advanced stage. The characteristic local peritoneal tumor spread is surgically limiting and followed by therapy resistant local (peritoneal) tumor recurrence as the predominant cause of cancer-associated death. Several molecular subclassification systems have been established [1][2][3][4][5] but did not lead to clinically relevant targeted therapies and over the last years little improvement in overall survival has been achieved.

Co-Association Gene Expression Clusters Associated with Tumor Spread Types
The aim of this analysis was to find a relevant over-expressed gene in patients with miliary tumor spread (associated with a significantly worse outcome) as defined in Auer et al. [6]. Therefore, our own RNA-sequencing data from tumors with known peritoneal tumor spread characteristics (miliary or non-miliary) were correlated with functional gene expression sub-clusters derived from a TCGA high grade serous ovarian cancer co-expression network. Clusters with up-regulated hub-genes already targeted by therapeutics in clinical studies were preferred, mainly to speed up potential clinical usage. As a first validation step the candidate sub-cluster was investigated via gene-signature analyses for its impact on overall survival in several other HGSOC cohorts with available gene-expression, clinicopathologic, and outcome data.
Whole transcriptome RNA-sequencing data from isolated tumor cells (i.e., enzymatically singularized and enriched for EpCAM protein abundance) of solid tissues and from ascites of patients with either miliary or non-miliary tumor spread were characterized according to three gene signatures: a 272 miliary-versus-non-miliary gene signature [6], a 211 ovarian-versus-tubal origin gene signature [12] (usually highly correlated with each other), and a 212 epithelial-mesenchymal (EM) gene signature [17]. Eight partly immortalized "normal" ovarian surface (four, three of them immortalized) and tubal secretory (four immortalized) epithelial cells were also analyzed accordingly. In Figure 1B results of these analyses are shown in triplots, clearly showing a discrimination in all three parameters of samples from patients with clinically evident miliary or non-miliary tumor spread (i.e., unambiguously classified intrasurgically [6]), most pronounced in ascites tumor cells. In solid tumor tissues no difference in the EM-status was revealed, similar to the non-immortalized OSE cells compared to the immortalized OSE cells ( Figure 1B, left triplot). This probably is due to the growth condition as solid tumor tissues (or as epithelia on the surface of ovaries), which forces cells to a more epithelial status. In Figure 1A typical representations of both peritoneal spread types, miliary and non-miliary, are shown.  [6], a 211 ovarian-versus-tubal origin gene signature (left edge) [12], and a 212 epithelial-mesenchymal (EM) gene signature (bottom edge) [17] and colored according tumor spread characteristic (green and blue, miliary; cyan and violet, non-miliary). For details see Material and Methods.  [6], a 211 ovarian-versus-tubal origin gene signature (left edge) [12], and a 212 epithelial-mesenchymal (EM) gene signature (bottom edge) [17] and colored according tumor spread characteristic (green and blue, miliary; cyan and violet, non-miliary). For details see Material and Methods. In Figure 2 the complete process and the results are shown, described briefly in the following sections. Whole transcriptome data of 303 HGSOC samples were normalized and used for building of a co-association network with about 4.3 million edges between the 20,501 genes with mutual information criterion as measure for the association (Figure 2A). Planar filtering of this network revealed a network of 47,581 edges connecting all genes in one large scale-free (degree distribution follows a power law) small world network ( Figure 2B). This network was used for multiscale clustering analysis to get sub-networks (clusters) significantly more associated internally compared to the total network. Hub-genes of these clusters were also identified ( Figure 2C). Finally, 653 hierarchically nested clusters were defined, with node numbers larger than ten and smaller than halve of the entire network ( Figure 2D). Zero to 245 hub-genes were identified for each cluster. These 653 clusters were used for gene-set analyses with samples from ascites (A/S; cf. Figure 1B) and solid tissues (P/M; cf. Figure 1B) between patients with miliary and non-miliary tumor spread ( Figure 2E; Table S1). Finally, significant values of both comparisons were combined and corrected for multiple testing by the False Discovery Rate (FDR) approach. Interestingly, 149 out of 201 significant clusters below an FDR of 0.01 showed a matching direction of global deregulation in both comparisons (A/S and P/M), i.e., up or down in both, compared to only 400 out of all 653 clusters (enrichment p-value < 0.001). As for therapeutic targeting overexpressed genes are much easier druggable, the clusters up-regulated in miliary tumors were searched for hub-genes already used as targets in clinical studies in decreasing order of significance (Table 1). Static plots (with regulation information, red, up in miliary, green, down in miliary, each with gene names labelled if significant and darker colored and always name-labelled, if hub genes) and interactive network-representations with genes (nodes) linked to GeneCards (http://www.genecards.org/) of the first three Clusters are linked in Table S1. Bold, selected cluster with Nectin 4 as hub-gene.
The predictive value was dichotomized at the 25 percentile, taking into account our observation that approximately 75% of FIGO III/IV HGSOC patients present with miliary tumor spread and according to published data of frequencies of putative ovarian and tubal origin of HGSOCs [15,20]. This gene signature was applied to six publicly available gene expression data sets of ovarian cancer samples [4,[21][22][23][24] meeting the following criteria: the clinicopathologic parameters grade, histology, FIGO stage, age at diagnosis, and residual tumor mass after debulking surgery had to be present; the first three parameters were used for selection of samples: Grade ≥ 2, serous histology, and FIGO stage ≥ III (in words: high grade late stage serous) the last three were used in multiple Cox-regression models to correct for these predictive clinicopathologic factors, together with the dichotomized Nectin 4 cluster gene signature value. A meta-analysis over these six independent data sets was significant for the dichotomized Nectin 4 cluster signature in an univariate setting ( Figure 2H left; p = 0.003), and a trend in multiple Cox regression models together with clinicopathologic parameters FIGO stage, age, and residual tumor ( Figure 2H right; p = 0.063, only one data set predicts favorable OS for Nectin 4 high samples in the multiple Cox-regression), showing an unfavorable outcome of patients with high Nectin 4 cluster expression.  The plot shows the node-degrees for all nodes (sorted according node-degree) in a log 10 -log 10 -scale. (B) Planar filtered-small world scale free-network (MEGENA, R-package). The plots shows the same as in

HGSOC Cell Line Models
Forty-eight ovarian cancer cell lines [25] were also characterized with above mentioned gene signatures for tumor spread and origin and the NECTIN4 gene expression ( Figure S1). Several high grade serous ovarian cancer cell lines with putative different spread and origin characteristics and appropriate Nectin 4 expressions (high in miliary of tubal origin, Caov-3, OVCAR-3, and OVKATE, and low in non-miliary of ovarian origin, OV-90, ES-2, and TYK-nu, with mean 36-fold higher Nectin 4 expression in the former) were identified, which could be used in further functional in vitro and in vivo (mouse) experiments.

Protein Expression of Nectin 4 (Assessed by Immunohistochemistry) and Impact on Overall Survival
Aim of this section was to validate the independent impact on overall survival of the candidate over-expressed hub-gene at the protein level (via immunohistochemistry staining of FFPE tissues) to strengthen the rational for using this gene product as target for a targeted therapy in selected patients with HGSOC.
Nectin 4 protein expression was assessed on 90 late stage (FIGO III/IV) HGSOC tissues by immunohistochemistry and scored according to percentage of positive tumor cells and staining intensities (H-score; 0, 1+, 2+, and 3+) (Figures 3 and S2). The Cohen's Kappa coefficient for the H-score between both raters was κ = 0.78 with ±1 as largest discrepancy. Samples with discordant H-scores were assessed together upon agreement. Percentage positive cells correlated with R = 0.91 and were averaged. Correlation of Nectin 4 protein expression to tumor spread types, miliary and non-miliary, was not possible as this information is not available from this cohort of patients. But gene expression analysis of all Nectin 4 cluster genes (c1_143, see Figure 2F) revealed clearly an up-regulation of the complete cluster in association to high Nectin 4 protein expression ( Figure 4B,C, the latter showing the log 2 fold-change enrichment barcode plot of the cluster). As there was no information available from literature concerning Nectin 4 impact on overall survival (OS) in HGSOC, we exploratively assessed the optimal cutoff by non-linear modeling of the Nectin 4 impact on OS by fractional polynomials Cox regression estimation [26] correcting for the known clinicopathologic factors age, FIGO stage, grade, and residual tumor mass after debulking surgery. In Figure 4A the corrected relative hazard against the percentage of Nectin 4 positive tumor cells is shown indicating a negative impact on OS in tumors with >50% and therefore this cutoff was used for outcome analyses. For all molecular biologic analyses (cf. Section 2.4. and Figure 5) a trichotomized Nectin 4 protein expression score was used (0, negative; 1, ≤50%; and 2, >50%) to increase statistical power. Patient characteristics related to Nectin 4 expression are presented in Table 2. Incomplete surgical resection was significantly associated with a high Nectin 4 score (p-trend = 0.020; Table 2). In addition, a high Nectin 4 score was significantly associated with a higher frequency of neoadjuvant chemotherapy mode (overall 13.8%), indicating a higher tumor burden at diagnosis, but not with pre-operative CA-125 or ascites volume ( Figure S3). All patients received a systemic platinum-based chemotherapy, either only adjuvantly of both, neo-adjuvantly and adjuvantly, the latter only if primary debulking surgery was not possible due to high tumor burden. The Nectin 4 score was not significantly associated to platinum-based chemotherapy resistance (defined as recurrent or progressive disease within six months after debulking surgery or neo-adjuvant chemotherapy start), with only a slightly higher frequency in the >50% Nectin 4 group (25.0% compared to each 16.7% for the 0% and ≤50% Nectin 4 groups).      Univariate Cox regression analyses (Table 3) and a Kaplan-Meier estimate indicate a negative impact of >50% Nectin 4 abundance on overall survival (HR 3.03, CI 95 1.37-6.68; p = 0.006; Figure 4D). Even after correction for age, FIGO stage, grade, and residual tumor mass in a multiple Cox regression analysis (Table 3), high Nectin 4 abundance showed an independent negative impact on overall survival (HR 3.62, CI 95 1.52-8.63; p = 0.004), illustrated as survival curves of the Cox regression model in Figure 4E. Adding CA-125 and chemotherapy mode to the multiple Cox regression model increased the HR for OS of high Nectin 4 expression to 4.64 (CI 95 1.52-14.13) (Table S4). CA-125 showed no impact on OS and chemotherapy mode only univariately (HR for neodajuvant chemotherapy 2.46, CI 95 1.01-6.04, p = 0.049). Similar analyses for progression free survival (Table S3) revealed a significant univariate impact (HR 2.60, CI 95 1.33-5.08; p = 0.005) which did not hold in the multiple Cox regression model (HR 1.92, CI 95 0.94-3.95; p = 0.074).  3 The optimal cutoff was assessed by non-linear modeling of the Nectin 4 impact on OS by fractional polynomials Cox regression estimation correcting for known clinicopathologic factors age, FIGO stage, grade, and residual tumor mass after debulking surgery. In Figure 4A the corrected relative hazard against the percentage of Nectin 4 positive tumor cells is shown indicating a negative impact on OS in tumors with >50%, and therefore this cutoff, >50% versus ≤50%, was used for outcome analyses; Bold, statistically significant.

Network of Associations of the Nectin 4 Driven Expression Clusters with Omics Data and Biological Interpretation
Aim of this section was to examine the biological impact of the candidate hub-gene on tumor tissues, the microenvironment (i.e., peritoneal), and systemically, mainly on the immune reactions. Therefore Nectin 4 protein expression was correlated to omics and medium dimensional data derived from tumor tissues, the ascites, and the serum of the patients by various statistical approaches, like pathway, network, and gene-set analyses. Finally, the complex associations between these data and results-correlated to the Nectin 4 expression-were presented in a richly annotated network representation.
To examine the impact of high Nectin 4 expression on tumor cells, the local intraperitoneal microenvironment and systemic characteristics of patients, the trichotomized Nectin 4 score (0%, ≤50%, and >50%, respectively) was correlated to (i) gene expression data of isolated solid tumor cells (ovarian and peritoneal tumor masses) and free floating ascitic tumor cells (single and spheroids) [6,12], to (ii) immune cell compositions in tumors, ascites, and serum of patients (FACS and immunofluorescence staining data) [7], to (iii) cyto-and chemokines in ascites and serum [7], and to (iv) untargeted proteomics data from cell free ascites. In Table 4 these data and results are summarized.
Ad (i) Differential gene expression analysis revealed 908 significantly deregulated genes (FDR<5%; thereof 643 up-regulated in Nectin 4 high expressing tumors). Using 2003 deregulated genes (FDR < 10%) for a pathway analysis applying enrichment and topological information with 199 KEGG pathways revealed eight pathways significantly deregulated (FDR < 5%; all inhibited; Figure S4 and Table S5). The major histocompatibility complexes I (MHC I) and II (MHC II) were down-regulated in many of these inhibited pathways (Figures S5-S12), therefore we assessed the expression of the single HLA genes in detail. In Figure S13 all expressed MHC I and II HLA genes are shown, correlated to the Nectin 4 score. Five of seven MHC I and six of 13 MHC II associated genes were down-regulated (FDR < 20%) in Nectin 4 high expressing tumors. Whereas only two of 15 non-coding (pseudo-or antisense-) HLA genes were down-regulated ( Figure S14). Table 4. Significant genes, proteins, and immune cell populations associated with the Nectin 4 score and significant correlations between clusters and these factors; and significant genes and pathways associated with NECTIN4 expression in HGSOC cell lines (n.a., not applicable). Using a gene-set enrichment analysis with the 653 TCGA HGSOC expression data derived-clusters (cf. above and Figure 2) revealed 20 significantly deregulated clusters (FDR < 5%; thereof 12 up-regulated in Nectin 4 high expressing tumors; Table S5). A further gene-set enrichment analysis with the gene-sets v6.2 from the Broad Institute [27] revealed 1828 of total 17,810 gene-sets significantly deregulated (FDR < 5%; 133 up-regulated in Nectin 4 high expressing tumors). The gene-sets "Chr22q12" and "Chr22q11" were five orders of magnitude more significant compared to all other gene-sets, indicating a chromosomal aberration at this region associated with Nectin 4 expression. To further examine this finding the log 2 fold changes associated to the Nectin 4 score of all genes mapped to chromosome 22 were plotted ordered relative to their positions. And indeed, significantly reduced log 2 fold changes can be seen between 22q11.21 and 22q12.3 chromosomal bands ( Figure S15). Another evidence for copy loss of a chromosomal region correlated with the Nectin 4 score arose from the analysis of cluster c1_496 (cf. Table 5 and Figure S26). Twenty out of 23 genes from this cluster map to chromosomal band 6q21 and two further genes to 6q16 and 6q22. An indeed, 6q is the most affected chromosomal loss region in ovarian cancer [4], indicating LOH as driving event during the co-association network generation from expression data and sub-cluster identification of c1_496. The expression of the cluster c1_496 is negatively correlated to the Nectin 4 score indicating a positive correlation of 6q21 LOH frequency with Nectin 4 expression. genes from this cluster map to chromosomal band 6q21 and two further genes to 6q16 and 6q22. An indeed, 6q is the most affected chromosomal loss region in ovarian cancer [4], indicating LOH as driving event during the co-association network generation from expression data and sub-cluster identification of c1_496. The expression of the cluster c1_496 is negatively correlated to the Nectin 4 score indicating a positive correlation of 6q21 LOH frequency with Nectin 4 expression.

Figure 5.
A zoomable version of this network is provided as pdf in the Supplement as Figure S39. Network of significant Nectin 4 correlated gene-expression co-association clusters (cf. Table 5) correlated to other omics (proteomics), FACS (immune cell populations) and Luminex (cyto/chemokines) data (cf. Table 4  . The type and source of the analytes/cell-populations is coded in the node-label and node border color, respectively. Analytes and cell-populations are given as node-labels (for complete lists of analyzed factors cf. Table 4 and Table S5). A comprehensive legend is shown on the left border of the figure. Table 5. List of Nectin 4 associated significantly deregulated HGSOC specific gene co-association clusters (derived from gene expression data of 303 HGSOC patients (TCGA), cf. Figure 2) with corresponding hub-genes and putative functions. The correlations with the Nectin 4 expression are indicated by UP/DOWN and colors red/green, respectively. These clusters were used for multi-omics correlation analyses shown in Figure 5 (cf. Table 4).   Furthermore, significantly deregulated genes were mapped onto four networks, (i) the TCGA expression derived co-association network ("TCGAnet", a HGSOC specific co-expression network, cf. above), (ii) the STRING v10.0 network ("STRING", a direct and functional (literature) driven protein-protein interaction network, [28], (iii) an experimentally verified protein-protein interaction network ("PPI"; CCSB Human Interactome database, HI-III, preliminary release 2.5, downloaded from http://interactome.dfci.harvard.edu), and (iv) a network of all connected KEGG pathways (KEGG super-pathway, "KEGG", [29]), and the maximum scoring sub-networks with at least 30 nodes (20 for the PPI network) determined. In Supplementary Figures S16-S19 these four sub-networks and in Supplementary Figures S20-S30 eleven of the smaller (<30 nodes) significantly deregulated TCGA clusters are shown and in Table 5 some hub genes and the putative functions of these high-scoring sub-networks and clusters are summarized. Interestingly, seven out of the 110 genes present in the four high-scoring sub-networks were from the 263 genes mapping to bands q11 and q12 on chromosome 22 ( Figure S31), which is a highly significant enrichment (p < 0.0001).

Network
To correlate these sub-networks and clusters (represented by the first principal component (PC1) of the expressions of the network genes) with other omics and medium dimensional data from tumor tissues, the microenvironment (i.e., ascites), and systemically (i.e., blood) derived data (cf. Table 4), significant correlations between PC1s of the high-scoring sub-networks and clusters with all analytes were determined (FDR < 10%) and are shown as comprehensive network in Figure 5. The validity of this approach is exemplified by the revealed correlations of the PC1 of the c1_358 TCGA cluster with-inter alia-eleven Kallikrein genes (including KLK6) with the soluble KLK6 abundance in ascites (bottom middle Figure 5) or of the PC1 of the PPI network (including BCAM) with the soluble BCAM abundance in ascites. Using this approach, sub-networks and clusters with similar functions were arranged close together even if the overlap of the sub-networks and clusters was rather small ( Figure S31), like clusters c1_747 and c1_782 (Figures S29 and S30) or cluster c1_315 and the TCGAand KEGG-sub-networks (Figures S16, S19 and S21, respectively).
Among the immune cell populations, mainly CD14 positive cells (e.g., monocytes) and several different natural killer (NK) cell populations (especially NK1 and NK2 [7,30]) in the tumor tissue, blood, and ascites were underrepresented in Nectin 4 positive tumors and were significantly correlated to many sub-networks and clusters (cf. Figure 5).

NECTIN4 Dependent Expression Differences in HGSOC Cell Lines
To estimate functional expression differences between each three NECTIN4-high (Caov-3, OVCAR-3, and OVKATE) and NECTIN4-low (OV-90, ES-2, and TYK-nu) expressing HGSOC cell lines, publicly available RNA-sequencing data were used for gene expression analyses [25]. Of about 18k expressed genes, 53 genes were higher and four lower expressed (FDR < 5%) in NECTIN4-high cell lines (or 718 higher and 194 lower, respectively, with a 20% FDR cutoff) (cf. Table 4, Table S5 "sign_RNAseq_CLs" and Figure S32). Again, genes from the NECTIN4-cluster (c1_143) introduced in Figures 2F and 4B were significantly enriched in high ranks, indicating a significant up-regulation of the complete cluster ( Figure 6B). L1CAM was the most up-regulated single gene, known to be a negative prognostic factor in ovarian cancer [31,32] (Figure 6A), therefore fitting to the worse outcome of patients with high Nectin 4 expressing tumors in our study.

Discussion
In search of new targets for which putative therapeutics have already been developed, a comprehensive analysis including various validation steps revealed new evidence for Nectin 4 as potential target for selected ovarian cancer patients. A clinically functional antibody against Nectin 4 has already been developed (enfortumab) and its relevance for epithelial cancers hypothesized [33],  Figure 4B). (C) Gene Ontology enrichment analysis with down-(G1) and up-regulated (G2) genes. (D) SPIA evidence plot of significantly deregulated KEGG pathways (blue line: cutoff FDR < 20%; for details cf. Figure S4). (E) KEGG pathway "Arrhythmogenic right ventricular cardiomyopathy" significantly activated in NECTIN4 high expressing cell lines. Node colors represent the log 2 fold changes (Table S5 "sign_KEGG_CLs"; for the other activated KEGG pathways cf. Figures S33-S38).

Discussion
In search of new targets for which putative therapeutics have already been developed, a comprehensive analysis including various validation steps revealed new evidence for Nectin 4 as potential target for selected ovarian cancer patients. A clinically functional antibody against Nectin 4 has already been developed (enfortumab) and its relevance for epithelial cancers hypothesized [33], especially as antibody-drug combination (monomethyl auristatin E, MMAE, a small molecule microtubule-disrupting agent), known as enfortumab vedotin or ASG-22ME. This antibody-drug combination is already under investigation in a phase I clinical trial (clinicaltrials.gov identifier ASG-22CE-13-2) and first results were presented at the European Society for Medical Oncology (ESMO) 2016 congress. So far, in ovarian cancer only little is known about Nectin 4 expression and its impact on outcome. Derycke et al. showed that Nectin 4 is overexpressed in 51.5% of high grade serous ovarian cancer tissues but found no impact on survival in a heterogeneous population comprising different histologies (serous, mucinous, endometrioid, and clear cell), grades, and FIGO-stages in Kaplan-Meier estimates [34]. Aberrant expression of Nectin 4 was shown for several cancer entities, including ovarian cancer with comparable expression frequencies in primary serous tumors, i.e., 47%, not discriminating between high and low grade, and higher frequencies in metastases, i.e., 79%, in serous ovarian cancer metastases (of unknown origin, peritoneal local or "real" distant) [33]. In a study evaluating NECTIN4 expression by qPCR and ELISA in 39 ovarian cancer patients, 21 subjects with benign ovarian pathologies and 25 healthy controls, Nabih et al. [35] demonstrated increased NECTIN4 mRNA expression in 97.4% of the ovarian cancer samples. In human ovarian cancer cells the significance of Nectin 4 in cell-cell adhesion, proliferation, and migration has been demonstrated by over-expression and silencing experiments [36]. Artificial higher NECTIN4 expression compared to biological or silenced NECTIN4 expression in cell line models was associated with a more epithelial phenotype [36], fitting to our results that cell lines with biological higher NECTIN4 expression are from the miliary type [12], shown to exhibit high epithelial characteristics [6].
Further, growth of orthotopically implanted breast cancer cells in nude mice can be inhibited by blocking PVRL4 (NECTIN4)-driven cell-to-cell attachment with monoclonal antibodies against PVRL4 [37].
Further information about expression or impact on outcome in (high grade serous) ovarian cancer is not available. In tumor specimens of 5673 triple negative breast cancer patients Nectin 4 has been identified as cell surface biomarker [38]. Nectin 4 was overexpressed in triple negative breast cancer and basal breast cancer patients. Further, these findings were strongly correlated to NECTIN4 mRNA expression. In multivariate analysis high NECTIN4 mRNA expression has been identified as independent poor prognosis factor for metastasis-free survival. To test Nectin 4 as therapeutic target an monoclonal anti-Nectin 4 antibody drug (monomethyl auristatin-E; MMAE) conjugate was investigated in vitro and in vivo. In both, localized and metastatic triple negative and Nectin 4-positive breast cancer this therapy led to rapid and long-lasting regression in mice xenograft models. In hepatocellular carcinoma (HCC) Nectin 4 was over-expressed in 68% HCC tissues and positive Nectin 4 expression was significantly correlated with tumor size, status of metastasis, vascular invasion and tumor-node-metastasis stage. Nectin 4 expression was also associated with worse recurrence-free and overall survival, univariately and multivariately [39]. In gastric cancer higher Nectin 4 expression was found in cancerous compared to normal gastric tissues and was associated to unfavorable outcome [40].
In the data presented above, NECTIN4 was identified as hub-gene in a highly significant co-association cluster between samples from patients with miliary and non-miliary peritoneal tumor spread behavior. Subsequent protein expression analyses revealed positive Nectin 4 expression in about 53% of late stage HGSOC tissues but impact on survival was only seen in those 13% of cases with >50% of positive tumor cells (Table 3 and Figure 4). A comprehensive network based multi-omics integration revealed several interesting biological and immunological characteristics of high Nectin 4 expressing tumors ( Figure 5): (i) expression of HLA genes-mainly from the MHC I complex-were negatively correlated to Nectin 4 expression, indicating an immunologic silencing strategy of these tumors, similar as we revealed in PD-L1 negative HGSOC tumors [41]; (ii) genes from chromosomal region 22q11-12 were highly over-represented and down-regulated in Nectin 4 high tumors (Figure S14), and this even seems to drive Nectin 4 associated high-scoring sub-networks (gene co-associations (TCGA clusters), functional (STRING) and experimentally verified PPI, and KEGG-super-pathway), i.e., seven out of 110 genes in these pathways were from this region, which is significantly over-represented ( Figure S30; p < 0001). It is known that loss of chromosomal region 22q with 79% is the second most common chromosomal loss in HGSOC (after 16q with 80% [4], which seems also correlated to Nectin 4 expression, indicated by cluster c1_496, which is comprised of only genes mapping to 6q16-22 and expression of this gene-cluster was significantly negatively correlated to Nectin 4 expression); (iii) a cluster with 11 Kallikrein genes (c1_358; Figure S22) was positively correlated to Nectin 4 expression and also soluble Kallikrein 6 in cell free ascites (untargeted proteomics). Kallikreins are a subgroup of serine proteases capable of cleaving peptide bonds in proteins, and are associated with an increased metastatic behavior and unfavorable survival in HGSOC [42] and seem to be attractive targets for ovarian cancer treatment [43]; (iv) higher soluble protein abundance in ascites and gene expression of basal cell adhesion molecule (BCAM) was also identified, a receptor for the extracellular matrix protein laminin α5 (LAMA5), known to be expressed in ovarian cancer tissues [44]. The function of BCAM in cancer is largely unknown, only in KRAS-mutated colorectal cancer it was shown, that the axis BCAM-LAMA5 mediates the recognition between tumor cells and the endothelium in the metastatic spread [45], and (v) several immune cell population changes related to the clusters, systematically in blood and locally in ascites and the tumor, i.e., reduced CD14 positive cell concentrations (e.g., monocytes) in ascites (negatively correlated to clusters c1_180 and c1_363), reduced concentrations of several natural killer cell populations in blood and the tumor [7] (negatively correlated to clusters c1_747, c1_782, c1_376, c1_335, and the high-scoring sub-networks PPI and STRING), and the cluster c1_143 (with Nectin 4 as hub-gene) correlated with several immune modulating cytokines, e.g., in serum negatively with IL6, the soluble interleukin-6 receptor (sIL6Ra), MIF, and CCL20 and in ascites negatively with, CCL25, and prolactin. The immune-suppressive IL10 was positively correlated with the Nectin 4 cluster (c1_143). For further correlations see Figure 5.
Interestingly, ATG5 as member of the negatively to the Nectin 4 expression correlated c1_496 cluster of 6q21 genes is a known key regulator of autophagosome formation. Next to autophagy activation, ATG5 is also involved in cell death progression [46]. Autophagy is abnormally activated during cancer cell metastasis and plays an important role in the maintenance of cancer cell viability [47]. Recently, in a mouse model ATG5 deficiency has been shown to be involved in tumorigenesis [48]. Downregulation of ATG5 transcription further leads to repression of macroautophagy activity, which promotes breast cancer cell metastasis [49]. In colorectal cancer tissues, ATG5 expression was also downregulated, which led to the assumption that ATG5 may function as a tumor suppressor [50].
Given the putative function of Nectin 4-Ca 2+ -independent cellular adhesion-the anti-Nectin 4 antibody might have additional effects in ovarian cancer besides being "the Trojan horse for an antineoplastic substance", i.e., potential direct effect through inhibition of the attachment of tumor cells on the peritoneal wall, which might inhibit local peritoneal tumor spread directly. First underlying data can be derived from cell line models [36] or in vivo mouse models (breast cancer [37]).
Further functional analyses uncovering the direct impact of Nectin 4 and the other-mostly surface-proteins from the identified cluster, probably involved in cell-cell or cell-matrix attachment, are necessary to enlighten the functional impact in HGSOC, especially for peritoneal tumor spreading.

Planar Filtered Gene Co-Association Network and Differential Cluster Expression Analyses
All analyses were conducted in R v3.5.0 [51] and CRAN or Bioconductor packages [52] (if not otherwise stated). Normalized read counts from RNA-sequencing data of 303 HGSOC tissues from The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/) were downloaded via R-package TCGA2STAT v1.2 and filtered for genes with counts per million > 0.1 in halve of samples (20,501 genes). Counts were transformed and normalized with the function voom ("transforms count data to log 2 -counts per million (logCPM) [with offset 0.5], estimates the mean-variance relationship and uses this to compute appropriate observational-level weights. The data are then ready for linear modelling") and the cyclic loess method ("Normalizes the columns of a matrix, cyclicly applying loess normalization to normalize each pair of columns to each other". Cyclic loess is a non-linear but less aggressive normalization method compared to quantile normalization) implemented in the Bioconductor R-package limma v3.30.7 [53]. The gene co-association network was built with Tool for Inferring Network of Genes (TINGe v1.062) using mutual information criterion as association measure, a list of known transcription factors (n = 1672), and following parameters (converting final mutual information to correlation values): "−C c −w −v −p 1e−3 −e 0.2" [54]. The planar filtered (thus scale free small world) network was built from the TINGe network using the calculate.PFN function from R-package Multiscale Clustering of Geometrical Network (MEGENA v1.3.5-2) [55]. Multiscale clustering analysis (MCA) and multiscale hub analysis (MHA) was performed with do.MEGENA function from the same R-package, allowing clusters of sizes larger than ten and smaller than halve of nodes in the network (i.e., 10,250) [55]. The final list of clusters (n = 653) was used for gene-set analyses comparing miliary and non-miliary samples (as defined in [6]) from ascites (n = 22; A, single ascites tumor cells; S, spheroidal ascites tumor cell aggregates) and solid tumor tissues (n = 20; P, primary, i.e., ovarian, tumor mass; M, metastasis, i.e., peritoneal tumor mass) of 21 patients (twelve miliary and nine non-miliary) [6,12]. Gene-set analyses were performed with mroast [56] function for AS-samples and with camera [57] function for PM-samples (limma) using gene-and sample-weights from the voomWithQualityWeights function (limma) [58]. Raw p-values of both comparisons were combined by the Fisher's method and corrected for multiple testing by the Benjamini-Hochberg (False Discovery Rate, FDR) method (R-package MetaDE v1.0.5). The same gene expression data and additional gene expression data of eight partly immortalized "normal" ovarian and tubal cell(-line)s were used for triplots (R-package ggtern v2.2.0) in Figure 1B, with calculated 272 miliary-versus-non-miliary gene signature (right edge) [6], 211 ovarian-versus-tubal origin gene signature (left edge) [12], and 212 epithelial-mesenchymal (EM) gene signature (bottom edge) [17] scores (always median first type up, i.e., miliary, ovarian, or epithelial, genes minus median second type up, i.e., non-miliary, tubal, or mesenchymal, genes). Values were scaled between zero and one before plotting with ggtern. Usually a false discovery rate (FDR) cutoff of 5% was chosen, but 10% (e.g., heatmap of differentially expressed genes, correlations of associations between different omics-results in Figure 4, or for some results in Table 4) or even 20% (e.g., pathway analyses for differences between Nectin 4 high and low cell lines) were also used in specific cases. The used FDR-cutoff is always indicated.

Cell Line Analyses
Ovarian cancer cell lines were analyzed as described in Auer et al. [12] and triplots plotted with Nectin 4 gene expression as third value (bottom edge), instead of the 212 epithelial-mesenchymal (EM) gene signature score. Cell lines were classified as of ovarian or tubal origin by the non-negative matrix factorization approach (NMF) using the 211 origin gene signature genes as described [12].

Overall Survival Meta-Analyses
Overall survival meta-analyses were performed with the curated ovarian cancer data sets from the Bioconductor R-package curatedOvarianData v1.12.0 [59]. In a first step, only samples of high grade (2/3) and late stage (FIGO stage III/IV) with serous histology were selected with further information for age and residual tumor mass after debulking surgery and used for univariate and multiple Cox-regression analyses (with age, FIGO-stage, and residual tumor mass as correcting factors) together with the dichotomized Nectin 4 cluster gene signature score (calculated as follows: median expression of all genes up in miliary, n = 31, minus median expression of all genes down in miliary, n = 4; dichotomized at the 25 percentile; Table S2). Results are shown in forest plots and combined hazard ratios with 95% confidence intervals and combined p-values.

Immunohistochemistry of Nectin 4 and Impact on Overall Survival
Whole tissue sections from ovarian tumor tissue were employed for immunohistochemistry (IHC). Staining was performed as described previously [60] using a polyclonal antibody against Nectin 4 (dilution 1:500, cat no. Ab192033; Abcam, Cambridge, MA, USA) that had already been implemented successfully in tissue samples of solid tumors. Antigen heat retrieval was performed by microwaving the slides in EDTA (0.01 M, pH 8.0). The sections were incubated at 4 • C overnight with primary antibodies. Placental tissue sections were used as positive controls. For negative control ovarian and placenta sections were incubated in absence of the primary antibody. All quantitative assessments (H-scoring, 0, 1+, 2+, and 3+, and percentage positive cells) were done by visual scoring independently by two trained IHC analysts specialized in gynecology. The Cohen's kappa coefficient (κ) was calculated for the inter-rater agreement according the H-score and discordant samples were looked at together to resolve the discrepancy. Percentage positive cells were averaged between both raters. Approval for this study was obtained by the ethical review board (Ethics Committee of the Medical University of Vienna: nos. 366/2003, 793/2011, and 1076/2018) and all patients signed an informed consent. Table 2 shows the characteristics of the 90 patients included in the IHC and survival analyses ( Table 3). Only patients with advanced stage, high grade serous EOC were included. Mean age at time of cytoreductive surgery was 60.0 years (SD 11.9 years). The median observation period for progression free survival was 56.8 months (IQR 27.3-62.9 months) and for overall survival 47.4 months (IQR 24.9-65.5 months), estimated from patients without events, respectively. Within the observation period, 40 (44.4%) patients died and 71 patients (78.9%) experienced tumor recurrence or progression.
Non-linear modelling of percent Nectin 4-positive cells with fractional polynomials Cox regression was performed with R-package mfp v1.5.2 [26] using fp(percNECTIN4, df = 4) together with age, FIGO stage, grade, and residual tumor. After dichotomization at >50% (cf. Figure 3A) progression free and overall survival was estimated by univariate and multiple Cox regression (coxph()-function from R). Kaplan-Meier estimates of univariate Cox regressions and survival curves of the multiple Cox regression model were plotted with R-function plot(survfit(..)), for the latter with all parameters averaged except the dichotomous Nectin 4 variable. p-values below 0.05 were considered as statistically significant.

Untargeted Proteomics of Cell Free Ascites
Cell free ascites samples were depleted using Pierce™ Top 12 Abundant Protein Depletion Spin Columns (Thermo Fisher Scientific, Waltham, MA, USA). Samples were depleted according to the protocol of the manufacturer. In-solution digestion of depleted samples was performed on 3 kDa MW cutoff filters (Nanosep with Omega membrane, Pall Austria Filter GmbH, Vienna, Austria). Total protein amount of 20 µg was used for digestion, which consists of protein reduction (dithiothreitol, Gerbu Biotechnik GmbH, Heidelberg, Germany), protein alkylation (2-iodoacetamide, Sigma-Aldrich, St. Louis, MO, USA) and protein digestion with a Trypsin/Lys-C Mix (MS grade; Promega Corporation, Madison, WI, USA) [61]. Untargeted proteomics analysis of digested peptide was conducted on a QExactive Orbitrap mass spectrometer (Thermo Fischer Scientific) coupled with a UltiMate 3000 RSLC nano System (Pre-column: Acclaim PepMap 100, C18 100 µm × 2 cm; Analytical column: Acclaim PepMap RSLC C18 75 µm × 50 cm; Dionex, Sunnyvale, CA, USA) [61]. A 54 min gradient from 1% to 40% solvent B was applied for chromatographic peptide separation. Solvent composition: A-7.9% water, 2% acetonitrile, 0.1% formic acid and B-97.9% acetonitrile, 2% water and 0.1% formic acid. The resolution on Orbitrap mass spectrometer was set to 70,000 for full MS scan and the MS2 scan was acquired at a resolution of 17,500, both at m/z 200. The MaxQuant 1.5.2.8 [62] software packet was used for protein identification and MS1 based label-free protein quantification (LFQ). The same criteria as described recently were applied for protein identification and MS1 based protein qualification [61]. Briefly, search criteria included a maximum of two missed cleavages and a maximal mass deviation of 5 ppm for peptide ions and of 20 ppm for fragment ions. Further, a minimum of two peptide identifications per protein (including one unique) was requested and an FDR of less than 0.01 was applied at both peptide and protein level.

Correlations to Omics and Other Data and Network Construction for Biological Interpretation
Differentially expressed genes associated with the Nectin 4 protein score (0, negative; 1, ≤50%; and 2, >50%) were assessed from raw RNA-sequencing reads [6,12] after cyclic loess normalization with function voomWithQualityWeights from R-package limma. Differentially abundance analyses of omics and other data (cf. Table 4) were performed also with limma after logarithmization (base 2) of the corresponding (normalized) read-outs (zeros were replaced by 0.1 before logarithmization). Differentially regulated KEGG-pathways were determined by a combined gene over-representation and perturbation analysis using Signaling Pathway Impact Analysis (SPIA), implemented in the Bioconductor R-package SPIA v2.32.0 [63], with 10,000 permutations and the normal inversion "norminv" p-value combination method. Correction for multiple testing was done by the False Discovery Rate (FDR) method (Benjamini-Hochberg). Significant pathways were illustrated with the Bioconductor R-package pathview v1.20.0 [64]. Over-representation analyses for clusters was performed with Quantitative Set Analysis for Gene Expression (QuSAGE), implemented in the R-package qusage v2.14.0 [65]. Gene-set analyses was performed with the Molecular Signatures Database (MSigDB) v6.2 [27] and the function camera [57] from R-package limma (independent variable: the trichotomized Nectin 4 score). Finally, the R-package dnet v1.1.4 [66] was used to find high-scoring sub-networks using the p-values from the differentially gene expression analysis above, either using the STRING v10 database of known (including literature-based) and predicted protein-protein interactions ("STRING"; [67]) or the more selective-only experimentally verified physical protein-protein interactions-CCSB Human Interactome-database ("PPI"; HI-III, preliminary release 2.5, downloaded from http://interactome.dfci.harvard.edu), a network of all connected KEGG pathways (KEGG super-pathway, "KEGG", [29]), and the complete scale-free small-world co-association gene expression network ("TCGA"; cf. Figure 2). p-value cutoff for the heuristically finding of the high-scoring sub-networks were optimized to get at least 30 (for the PPI network, 20) nodes.
For the correlations shown in Figure 4, the four high-scoring subnetworks (STRING, PPI, KEGG, and TCGA) and all significant with Nectin 4 associated TCGA-clusters with less than 50 nodes (n = 12) were used. All genes with available gene expression data from these clusters were summarized over all samples using the first principal component (PC1) of a principal component analysis (PCA) using the prcomp(center = TRUE, scale = TRUE) function from R. In Figure 4 (right border) a boxplot of all 16 percentage of explained variations (PEVs) of all PC1s are shown. Using these PC1s as independent variables, all omics and other data sets were screened for significant correlations (using limma) and all significant correlations (FDR < 10%) were collected (including their orientation, either positively or negatively correlated). These correlations were than plotted in a network as colored edges (blue, negatively; red, positively correlated) between all PC1s of the four high-scoring sub-networks and twelve significant clusters and all analytes (immune cell populations, chemo/cytokines, and cell-free proteins). In Table 4 these data and results are summarized. In addition, the correlations between the 16 high-scoring sub-networks and clusters were assesses by graphical Gaussian modelling (GGM) using R-package GGMselect v0.1-12.1) [68]. The principle behind GGM is to use partial correlations as a measure of independence of any two PC1s which allows for distinguishing direct from indirect correlations. Raw RNA-sequencing data are available at ArrayExpress with accession no. E-MTAB-7674.

Conclusions
This study identified Nectin 4 as correlated to extensive, miliary, peritoneal tumor spread in HGSOC patients and as significantly associated with survival. A clinically functional antibody (antibody-drug combination) against Nectin 4 has already been developed, its relevance for epithelial cancers hypothesized and first encouraging antitumor activities presented in metastatic urothelial cancers. Our study provides comprehensive biological data highlighting that this targeted therapy might be an interesting approach in selected patients with high grade serous ovarian cancer.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-6694/11/5/698/s1, Table S1: Table of sub-clusters differentially expressed in tumor cells of patients with miliary peritoneal tumor spread compared to patients with non-miliary peritoneal tumor spread (AS, ascites tumor cells; PM, solid tumor tissues), Table S2: Genes of sub-cluster c1_143 with information of different expression between miliary and non-miliary solid (PM) or ascites derived floating (AS) tumor cells and usage in the gene signature for survival validation as shown in Figure 2H, Table S3: Univariate and multiple Cox regression analysis for progression free survival in 90 patients with late stage (FIGO III/IV) high grade serous ovarian cancer, Table S4: Univariate and multiple Cox regression analysis for overall survival in 90 patients with late stage (FIGO III/IV) high grade serous ovarian cancer including pre-operative CA 125 levels and chemotherapy mode, Table S5: Tables of analytes (cf.  Table 4) with their associations with the protein Nectin 4 score (0, negative; 1, ≤50%; 2, >50%), Figure S1: 48 ovarian cancer cell lines characterized according the miliary-versus-non-miliary and the origin gene signatures (right and left edges, respectively) and Nectin 4 gene expression, Figure S2: Correlation of staining intensities and percentages of Nectin 4 positive cells on tumor tissues of 90 HGSOC patients, Figure S3: Correlation of CA 125 (log 10 U mL −1 ), amount of ascites (log 10 mL), and chemotherapy mode (ratio neoadjuvant to all (neoadj.+adjuvant)) with the Nectin 4 score (X-axis). p-values are not corrected for multiple testing, Figure S4: SPIA evidence plot of significantly deregulated KEGG pathways. Each pathway is represented by one dot, Figures S5-S12: Significantly downregulated KEGG pathways. Colors represent log 2 fold changes according correlations to the Nectin 4 score. White nodes were not represented in the RNA-sequencing data, either because they were not (reliably) expressed or could not be mapped, Figure S13: Coding HLA genes (major histocompatibility complexes I (MHCI) and II (MHCII)) correlated to the Nectin 4 score. Headers in red indicate significance at the FDR 20% level, Figure  S14: Non-coding (pseudo-or antisens-) HLA genes correlated to the Nectin 4 score. Headers in red indicate significance at the FDR 20% level, Figure S15: Log 2 fold-changes associated to the Nectin 4 score of all expressed genes on chromosome 22. Genes were ordered according their relative positions on chromosome 22. The length of the indicated chromosomal bands correspond not to the physical length of the bands on the chromosome. Red labeled are the seven genes from the high-scoring sub-networks, significantly associated to the Nectin 4 score (cf. Figure S31), Figure S16. Nectin 4 score associated high-scoring sub-network of the TCGA co-association gene expression network TCGAnet. Color: log 2 fold changes according Nectin 4 score associations; cave: all log 2 fold change values are negative, therefore the color bar is from −2 to 0, Figure S17: Nectin 4 score associated high-scoring sub-network of the functional protein-protein interaction network STRING v10. Color: log 2 fold changes according Nectin 4 score associations, Figure S18: Nectin 4 score associated high-scoring sub-network of the experimentally verified protein-protein interaction network (CCSB Human Interactome database, HI-III, preliminary release 2.5). Color: log 2 fold changes according Nectin 4 score associations, Figure S19: Nectin 4 score associated high-scoring sub-network of the KEGG pathway super-network (all KEGG pathways connected). Color: log 2 fold changes according Nectin 4 score associations, Figures S20-S30: TCGA co-association clusters significantly associated with the Nectin 4 score. Only clusters with less than 50 nodes are shown, and these clusters are together with cluster c1_143 (cf. Figure 3) and the high-scoring networks (TCGAnet, STRING, PPI, and KEGG; cf. Figures S16-S19) the basis of the integrative network shown in Figure 4. Color: log 2 fold changes according Nectin 4 score associations. If the association is significant (FDR < 5%), the node is framed by a black line, Figure  S31: Overlap of 110 genes from the high-scoring networks (cf. Figures S16-S19) with genes from chromosomal bands 22q11-q12 (cf. Figure S15), Figure S32: Heatmap of differentially expressed genes between three NECTIN4 high (Caov-3, OVCAR-3, and OVKATE) and three NECTIN4 low (OV-90, ES-2, and TYK-nu) expressing HGSOC cell lines (cutoff: FDR < 10% and log 2 fold change > 2), Figures S33-S38: Significantly (20% FDR) activated KEGG pathways. Colors represent log 2 fold changes according correlations to NECTIN4 expression. White nodes were not represented in the RNA-sequencing data, either because they were not (reliably) expressed or could not be mapped, Figure S39