Abstract Purpose Ovarian cancer is the leading cause of gynecologic cancer-related death worldwide. Early diagnosis of ovarian cancer can significantly improve patient prognosis. Hence, there is an urgent need to identify key diagnostic and prognostic biomarkers specific for ovarian cancer. Because high-grade serous ovarian cancer (HGSOC) is the most common type of ovarian cancer and accounts for the majority of deaths, we identified potential biomarkers for the early diagnosis and prognosis of HGSOC. Methods Six datasets ([38]GSE14001, [39]GSE18520, [40]GSE26712, [41]GSE27651, [42]GSE40595, and [43]GSE54388) were downloaded from the Gene Expression Omnibus database for analysis. Differentially expressed genes (DEGs) between HGSOC and normal ovarian surface epithelium samples were screened via integrated analysis. Hub genes were identified by analyzing protein–protein interaction (PPI) network data. The online Kaplan-Meier plotter was utilized to evaluate the prognostic roles of these hub genes. The expression of these hub genes was confirmed with Oncomine datasets and validated by quantitative real-time PCR and Western blotting. Results A total of 103 DEGs in patients with HGSOC—28 upregulated genes and 75 downregulated genes—were successfully screened. Enrichment analyses revealed that the upregulated genes were enriched in cell division and cell proliferation and that the downregulated genes mainly participated in the Wnt signaling pathway and various metabolic processes. Ten hub genes were associated with HGSOC pathogenesis. Seven overexpressed hub genes were partitioned into module 1 of the PPI network, which was enriched in the cell cycle and DNA replication pathways. Survival analysis revealed that MELK, CEP55 and KDR expression levels were significantly correlated with the overall survival of HGSOC patients (P < 0.05). The RNA and protein expression levels of these hub genes were validated experimentally. Conclusion Based on an integrated analysis, we propose the further investigation of MELK, CEP55 and KDR as promising diagnostic and prognostic biomarkers of HGSOC. Keywords: high-grade serous ovarian cancer, integrated analysis, bioinformatic analysis, differentially expressed genes, survival, biomarker Introduction Ovarian cancer is the leading cause of gynecologic cancer-related death and the fifth most common cause of cancer death in the United States. In 2018, approximately 22,240 new ovarian carcinoma cases and 14,070 associated deaths were expected in the United States, corresponding to almost 39 deaths per day.[44]^1 Because of the location of the ovaries and the lack of symptoms of early-stage disease, approximately 70% of ovarian cancer patients present with advanced disease (FIGO stage III/IV) and have a poor prognosis.[45]^1^–[46]^4 In contrast, if ovarian cancer could be diagnosed at a local stage, the 5-year relative survival rate would exceed 90%.[47]^1^,[48]^2 Hence, there is an urgent need to identify key early diagnostic and prognostic biomarkers specific for ovarian cancer that could significantly impact patient survival. Approximately ninety percent of ovarian cancers are epithelial ovarian cancers (EOCs), which are divided into four main histologic subtypes—serous, mucinous, endometrioid, and clear cell.[49]^3^,[50]^5^–[51]^7 The vast majority (70%) of EOCs are serous.[52]^8 Serous ovarian cancer can be categorized as high-grade serous ovarian cancer (HGSOC, grade 2 or 3) or low-grade serous ovarian cancer (LGSOC, grade 1), which have different pathogeneses and clinicopathologic features.[53]^9^–[54]^13 HGSOC is the most common type of EOC and has aggressive behavior, thereby accounting for the majority of deaths. Thus, we focused on HGSOC in this study. Recent advances in microarray-based profiling and high-throughput biological sequencing technologies have yielded considerable abilities to identify differentially expressed genes (DEGs) and discover potential biological mechanisms, thus allowing the identification of promising biomarkers for cancer diagnosis, treatment and prognosis. In this study, gene expression profiling data for human HGSOC and normal ovarian surface epithelium (OSE) samples were used to identify DEGs and analyzed their underlying biological functions by functional and pathway enrichment analyses. Then, protein-protein interaction (PPI) networks and Cytoscape were utilized to identify the hub genes, which were found to be closely related to the pathogenesis and progression of HGSOC; thus, these genes may play a vital role in the early diagnosis of HGSOC. Finally, the prognostic value of the hub genes in HGSOC was further confirmed by survival analysis, which can be used to identify predictors of poor clinical outcomes for patients with HGSOC. A flowchart of the analysis process is shown in [55]Figure 1. Figure 1. [56]Figure 1 [57]Open in a new tab Flowchart of the analysis process. Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction. Materials And Methods Gene Expression Profiling Data Gene expression microarray data ([58]GSE14001, [59]GSE18520, [60]GSE26712, [61]GSE27651, [62]GSE40595, and [63]GSE54388) were downloaded from the National Center for Biotechnology Information Gene Expression Omnibus (GEO) database ([64]http://www.ncbi.nlm.nih.gov/geo/). All included datasets met the following criteria: (1) they contained human HGSOC and normal OSE tissue samples; (2) the data were subjected to expression profiling by array analysis; and (3) they contained at least ten samples. These profile datasets were based on the Affymetrix [65]GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array), except for [66]GSE26712, which was based on the Affymetrix [67]GPL96 platform (Affymetrix Human Genome U133A Array). Six datasets including a total of 359 tissue samples (318 HGSOC samples and 41 normal OSE samples) were chosen for analysis. Integrated Analysis Of Microarray Datasets Raw data from [68]GSE14001, [69]GSE18520 and [70]GSE27651 were first subjected to base 2 logarithmic conversion for further analysis in R software. The BetweenArray normalization function in the limma package of R was applied to normalize the matrix data from each GEO dataset.[71]^14 After normalization, DEGs between HGSOC and OSE were screened by using empirical Bayes methods within the limma package.[72]^15 Based on the robust rank aggregation (RRA) method,[73]^16 the RobustRankAggreg package in R was utilized to integrate the DEGs identified in the six datasets. The following threshold values were established for screening DEGs: |log[2]FC| ≥ 2, where FC is the fold change; adjusted P < 0.05; and P < 0.05. Functional Enrichment Analysis Of DEGs We performed Gene Ontology (GO) enrichment analysis of the DEGs using the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8) ([74]http://string-db.org/).[75]^17 GO analysis annotates genes with respect to three independent ontologies—biological process (BP), cellular component (CC), and molecular function (MF).[76]^18 To elucidate potential pathways associated with the DEGs, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted in the clusterProfiler package[77]^19 to annotate genes with pathways.[78]^20 The false discovery rate (FDR) was obtained by adjusting the P value according to the Benjamini-Hochberg method.[79]^21 Genes with P < 0.05 and adjusted P < 0.05 or FDR < 0.05 were considered statistically significant. PPI Network Construction And Module Analysis The STRING online database (version 11.0) ([80]http://string-db.org/) was used to assess potential interactions among the DEGs.[81]^22 PPIs with an interaction score ≥ 0.4 (medium confidence) were utilized to construct the PPI network. Cytoscape software (version 3.7.0) was used to visualize and analyze the degree of connectivity to identify hub genes in the PPI networks.[82]^23 According to the degree of connectivity,[83]^24 we screened the top 10 hub genes. To detect the densely connected protein complexes in the PPI network, the Molecular Complex Detection (MCODE) app from Cytoscape was utilized with the default parameters to identify modules.[84]^25 Then, for significant module 1, we performed further GO and KEGG pathway enrichment analyses as previously described. Survival Analysis The Kaplan-Meier Plotter website ([85]www.kmplot.com/ovar) was utilized to validate the prognostic role of the ten hub genes in ovarian cancer patients. This website includes data for 2190 ovarian cancer samples on Affymetrix microarrays.[86]^26 We selected survival information for patients with HGSOC (grades 2 and 3) from multiple datasets (all available on the website) for the analysis.[87]^27 The patients were divided into two groups based on the best cutoff for gene expression (high vs low). The overall survival (OS) rates of the two groups were analyzed and Kaplan-Meier survival plots were then generated. Then, we performed a subgroup analysis considering stage, grade, TP53 mutation and treatment to understand how the expression of the identified hub genes impacts OS. Hazard ratios (HRs) with 95% confidence intervals (CIs) were calculated to identify protective (HR < 1) or risk genes (HR > 1), and the survival curves were plotted to visualize the relationships. A log rank P < 0.05 was set as the cutoff criterion. Validation Of Key Genes Hub genes were identified as the top 10 nodes in the PPI network. To confirm the reliability of these detected genes, we evaluated their expression in ovarian serous adenocarcinoma and normal ovarian tissues using datasets from Oncomine ([88]www.oncomine.org). In addition, the expression of the ten hub genes was experimentally validated by quantitative real-time PCR (qRT-PCR). Among the hub genes, EPCAM (epithelial cell adhesion molecule), ZWINT (ZW10-interacting kinetochore protein), DLGAP5 (DLG-associated protein 5) and KDR (kinase insert domain receptor) protein expression levels were confirmed by Western blotting. Clinical Samples With the approval of the Ethics Review Committee of Peking Union Medical College Hospital, Chinese Academy of Medical Sciences (ZS-1771), twenty-two HGSOC and twenty-two normal ovarian tissue samples were collected during initial operations between 2017 and 2018. All subjects gave written informed consent in accordance with the Declaration of Helsinki. HGSOC samples were obtained from primary ovarian cancer patients who had not previously received chemotherapy. Normal ovarian tissues were obtained from patients who underwent a total hysterectomy and bilateral salpingo-oophorectomy for benign uterine diseases (uterine prolapse or uterine leiomyoma) or precancerous lesions of the uterine cervix. [89]Table 1 summarizes the clinicopathological characteristics of the patients with HGSOC included in the study. All tissue specimens were pathologically confirmed before inclusion. All the fresh samples were frozen in liquid nitrogen and stored at −80°C. Table 1. Clinicopathological Characteristics Of The Patients With HGSOC Included In The Study Patient Age, y Histology Type Grade FIGO Stage Residual Lymphectomy Lymph Nodes Metastasis 1 41 Serous High IIIC Yes No NA 2 53 Serous High IIIC Yes No NA 3 55 Serous High IV Yes Yes + 4 51 Serous High III No Yes − 5 38 Serous High IA NA Yes − 6 55 Serous High IC No Yes − 7 63 Serous High IIIC Yes No NA 8 53 Serous High II Yes No NA 9 50 Serous High IV No Yes + 10 80 Serous High IIIC Yes No NA 11 40 Serous High III Yes Yes + 12 51 Serous High IIIC Yes Yes + 13 63 Serous High III Yes No NA 14 50 Serous High IIB No Yes − 15 49 Serous High IIIC Yes Yes + 16 48 Serous High IIIB Yes Yes − 17 63 Serous High III Yes No NA 18 59 Serous High IIIB Yes Yes − 19 65 Serous High IV No Yes − 20 52 Serous High III No Yes + 21 62 Serous High IV Yes No NA 22 56 Serous High III NA Yes + [90]Open in a new tab RNA Isolation And qRT-PCR Total RNA was extracted using TRIzol Reagent (Invitrogen) according to the manufacturer’s protocol. After the RNA concentration was measured with the Genova Nano Micro-volume Spectrophotometer (Jenway), cDNA was synthesized using GoScript^TM Reverse Transcriptase (Promega). qRT-PCR was performed using GoTaq^® qPCR (Promega) with the LightCycler^® 480 System (Roche) in accordance with the manufacturer’s instructions. The qRT-PCR primers are shown in [91]Supplementary Table 1. The conditions for all qRT-PCRs were as follows: Hot Start Taq activation for 1 min at 95°C, followed by 45 cycles of 95°C for 10 seconds, 58°C for 15 seconds and 72°C for 15 seconds, with a final step at 55°C for 1 min. The relative expression of target genes was calculated by the 2^−△△Ct method using GAPDH as the internal control. Western Blot Analysis Tissues were lysed in RIPA buffer with Halt^TM Protease and Phosphatase Inhibitor Cocktail (Thermo Scientific). Protein concentration was determined by the Bradford method using the Enhanced BCA Protein Assay Kit (Beyotime). Equal amounts of protein from each sample were separated in NuPAGE^TM 4–12% Bis-Tris Protein Gels (Invitrogen) and transferred to PVDF membranes (0.2 µm pore size) using iBlot^® 2 Transfer Stacks (Invitrogen). After being blocked with 5% nonfat milk in TBS with 0.1% Tween 20 for 1 hr at room temperature, the membranes were incubated overnight at 4°C with primary antibodies against the following proteins: EPCAM (1:1000; Cell Signaling #2929), ZWINT (1:2000; Abcam #71982), DLGAP5 (hepatoma upregulated protein, HURP; 1:2000, Abcam #70744), KDR (VEGFR2, 1:1000; Cell Signaling #2479) and β-actin (1:2000; ZSGB-BIO #TA-09). The membranes were subsequently incubated with a horseradish peroxidase (HRP)-conjugated secondary antibody at room temperature for 1 hr. Then, immunoreactive bands were detected with Immobilon Western Chemiluminescent HRP substrate (Millipore). Statistical Analysis GraphPad Prism 6.0 software was used to analyze the experimental data. The results are shown as the mean ± SEM. The statistical significance of differences between two groups was evaluated using Student’s t-test (two-tailed). P < 0.05 was considered to indicate statistical significance. Results Identification Of DEGs Information on the six GEO datasets included in the current study is provided in [92]Table 2.[93]^28^–[94]^34 The six datasets included a total of 318 HGSOC samples and 41 normal OSE samples. [95]Supplementary Table 2 presents detailed information on the samples in each included dataset. The data from each GEO dataset were normalized, and the results are presented in [96]Supplementary Figure 1. Using |log[2]FC| ≥ 2, adjusted P < 0.05 and P < 0.05 as cutoff criteria, we successfully screened a total of 103 DEGs via integrated analysis of the six GEO datasets ([97]Supplementary Table 3). Of the 103 DEGs, 28 were significantly upregulated and 75 were downregulated in HGSOC tissues compared to normal OSE tissues. The DEGs obtained from each GEO dataset are shown in [98]Figure 2. [99]Figure 3 presents the 103 DEGs identified by integrated analysis of the six GEO datasets based on the RRA method. [100]Supplementary Figure 2 shows the heat map of these 103 integrated DEGs in HGSOC and normal OSE tissues acquired from the GEO datasets. [101]Supplementary Table 4 shows the detailed expression data for each included sample. Table 2. Information On The Six GEO Datasets Included In The Current Study Dataset Reference Platform Number Of Samples (HGSOC/OSE) [102]GSE14001 Tung et al[103]^28 [104]GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array 13 (10/3) [105]GSE18520 Mok et al[106]^29 [107]GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array 63 (53/10) [108]GSE26712 Bonome et al[109]^30 Vathipadiekal et al[110]^31 [111]GPL96 [HG-U133A] Affymetrix Human Genome U133A Array 195 (185/10) [112]GSE27651 King et al[113]^32 [114]GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array 28 (22/6) [115]GSE40595 Yeung et al[116]^33 [117]GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array 38 (32/6) [118]GSE54388 Yeung et al[119]^34 [120]GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array 22 (16/6) [121]Open in a new tab Abbreviations: HGSOC, high-grade serous ovarian cancer; OSE, ovarian surface epithelium. Figure 2. [122]Figure 2 [123]Open in a new tab DEGs obtained from each GEO dataset ([124]GSE14001, [125]GSE18520, [126]GSE26712, [127]GSE27651, [128]GSE40595, and [129]GSE54388). The red points represent upregulated genes screened with |log[2]FC| ≥ 2 and adjusted P < 0.05. The green points represent downregulated genes screened with |log[2]FC| ≥ 2 and adjusted P < 0.05. The black points represent genes with no significant difference in expression. Abbreviations: DEGs, differentially expressed genes; GEO, Gene Expression Omnibus; FC, fold change. Figure 3. [130]Figure 3 [131]Open in a new tab Heat map of the 103 DEGs identified in the integrated microarray analysis. Each column represents one dataset, and each row represents one gene. The number in each rectangle is the log[2]FC value. The color gradient from blue to red represents the progression from down- to upregulation. Abbreviations: DEGs, differentially expressed genes; FC, fold change. Functional Enrichment Analysis Of The DEGs GO enrichment analysis for the three GO categories (BP, CC and MF) was performed in DAVID. As shown in [132]Figure 4A, a total of nine GO terms with P < 0.05 and FDR < 0.05 were enriched ([133]Supplementary Table 5). The 28 upregulated genes were significantly enriched in multiple BPs related to cell division and cell proliferation. The 75 downregulated genes were closely correlated with the ethanol oxidation term in the BP category; the extracellular region, proteinaceous extracellular matrix and extracellular exosome terms in the CC category; and the frizzled binding term in the MF category. Subsequently, we conducted KEGG pathway analysis the integrated 103 DEGs with the clusterProfiler package to explore potential enriched pathways. The downregulated genes mainly participated in the Wnt signaling pathway and diverse metabolism-associated signaling pathways, such as retinol metabolism, tyrosine metabolism, and drug metabolism/cytochrome P450 ([134]Figure 4B and [135]Supplementary Table 6). However, significantly enriched KEGG pathways were not identified for the upregulated genes. Figure 4. [136]Figure 4 [137]Open in a new tab Functional enrichment analysis of the 103 integrated DEGs. (A) GO annotation. The y-axis shows significantly enriched GO terms, and the x-axis shows the enrichment factor and -log[10]P values. (B) KEGG pathway enrichment analysis. The y-axis shows significantly enriched KEGG pathways, and the x-axis shows the enrichment factor and -log[10]P values. Note: The enrichment factor refers to the ratio of the number of DEGs enriched in a GO term/KEGG pathway to the total number of annotated genes enriched in the GO term/KEGG pathway. Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular component; MF, molecular function. PPI Network Construction And Module Analysis The STRING website was used to filter a total of 99 DEGs (28 upregulated and 71 downregulated) with a combined interaction score of > 0.4 and construct the PPI network, which comprised 65 nodes and 226 interactions ([138]Figure 5A and [139]Supplementary Table 7). The NetworkAnalyzer app in Cytoscape was used to calculate the node degrees to identify hub genes.[140]^23 After analyzing the data from STRING using the NetworkAnalyzer app, we screened the top 10 hub nodes according to node degree.[141]^24 [142]Table 3 presents detailed information for the 10 hub nodes—EPCAM, ALDH1A1 (aldehyde dehydrogenase 1 family member A1), ZWINT, BUB1B (BUB1 mitotic checkpoint serine/threonine kinase B), NEK2 (NIMA-related kinase 2), DLGAP5, MELK (maternal embryonic leucine zipper kinase), CEP55 (centrosomal protein 55), CKS2 (CDC28 protein kinase regulatory subunit 2), and KDR. Among these hub genes, only ALDH1A1 and KDR were downregulated. Figure 5. [143]Figure 5 [144]Open in a new tab PPI network construction and module analysis of DEGs associated with HGSOC. (A) PPI network of the integrated DEGs visualized by Cytoscape software. The node size represents the node degree (a larger node indicates a higher degree). The width of the edge indicates the combined score (a wider edge indicates a higher combined score). The color of the edge represents the EdgeBetweenness (a more orange edge indicates a higher EdgeBetweenness). (B-E) MCODE module screening for the DEGs, including module 1 (MCODE score = 7, nodes = 9, edges = 69), module 2 (MCODE score = 4, nodes = 5, edges = 28), module 3 (MCODE score = 3, nodes = 4, edges = 17), and module 4 (MCODE score = 2, nodes = 3, edges = 9). The red rectangle indicates the seed gene, and the blue circles indicate the clustered genes. Abbreviations: PPI, protein–protein interaction; DEGs, differentially expressed genes; HGSOC, high-grade serous ovarian cancer; MCODE, Molecular Complex Detection. Table 3. Information On The Ten Hub Nodes Gene Name Node Degree Betweenness Centrality Closeness Centrality Category EPCAM 10 0.36360281 0.34645669 Upregulated gene ALDH1A1 9 0.47501762 0.32592593 Downregulated gene ZWINT 8 0.00669486 0.22335025 Upregulated gene BUB1B 8 0.00669486 0.22335025 Upregulated gene NEK2 8 0.30443975 0.27160494 Upregulated gene DLGAP5 8 0.00669486 0.22335025 Upregulated gene MELK 8 0.00669486 0.22335025 Upregulated gene CEP55 8 0.00669486 0.22335025 Upregulated gene CKS2 8 0.00669486 0.22335025 Upregulated gene KDR 8 0.16213254 0.30985915 Downregulated gene [145]Open in a new tab Next, we performed module analysis in the MCODE app from Cytoscape to detect significant protein complexes in this PPI network and obtained four modules ([146]Figure 5B-[147]E, [148]Supplementary Table 8). Seven of the selected hub nodes were contained in module 1, suggesting that this module might play a key role in this PPI network. Thus, we performed further GO and KEGG enrichment analyses for module 1 ([149]Supplementary Tables 9 and [150]10). As shown in [151]Figure 6A, module 1 was closely correlated with a variety of GO terms, such as cell division, cell proliferation and mitotic nuclear division in the BP category; kinetochore and condensed chromosome kinetochore in the CC category; and protein binding in the MF category. The KEGG pathways enriched by the genes in module 1 included cell cycle and DNA replication ([152]Figure 6B). Figure 6. [153]Figure 6 [154]Open in a new tab Functional enrichment analysis of the genes partitioned into module 1. (A) GO annotation. The y-axis shows significantly enriched GO terms, and the x-axis shows the enrichment factor and -log[10]P values. (B) KEGG pathway enrichment analysis. The y-axis shows significantly enriched KEGG pathways, and the x-axis shows the enrichment factor and -log[10]P values. Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular component; MF, molecular function. Survival Analysis To explore the association between the aforementioned potential target genes and survival, the OS rates were calculated with the Kaplan-Meier Plotter website. A total of 1144 HGSOC patients from 9 datasets (TCGA, [155]GSE9891, [156]GSE26193, [157]GSE63885, [158]GSE18520, [159]GSE30161, [160]GSE14764, [161]GSE23554 and [162]GSE15622) were included in the OS analysis. Patients were stratified by low and high expression of each gene based on the best cutoff value. In the survival analysis ([163]Supplementary Table 11), the expression levels of MELK, CEP55 and KDR were significantly correlated with the OS rates of patients with HGSOC (P < 0.05). Additionally, OS analyses of prespecified subgroups based on stage, grade, TP53 mutation and treatment were also performed ([164]Supplementary Table 12). As shown in [165]Figure 7A–[166]C, a high MELK expression level was associated with a high OS rate for HGSOC patients, especially for those with advanced-stage disease (stages 3 and 4), grade 3 disease, and TP53 mutations and those administered platinum-containing chemotherapy regimens. High CEP55 expression afforded a poor OS, particularly for HGSOC patients in stage 3, without a TP53 mutation and treated with chemotherapy containing Taxol. Similarly, elevated KDR expression was significantly associated with a poor prognosis for HGSOC patients with grade 2 disease or under chemotherapy with Taxol ± platinum. Taken together, these findings indicate that MELK, CEP55 and KDR might represent important prognostic factors of survival for patients with HGSOC. Figure 7. [167]Figure 7 [168]Open in a new tab Survival curves and subgroup analyses of the OS rates of patients with HGSOC stratified by high and low expression of MELK (A), CEP55 (B) and KDR (C). Abbreviations: OS, overall survival; HGSOC, high-grade serous ovarian cancer. Validation Of Key Genes Oncomine was used to analyze the expression of the ten selected hub genes in ovarian serous adenocarcinoma and normal ovarian tissues[169]^35^–[170]^38 ([171]Table 4 and [172]Figure 8). Our findings regarding the expression of the ten selected genes were consistent with the expression data in the Oncomine datasets; P < 0.05 indicated statistical significance. Specifically, EPCAM, ZWINT, BUB1B, NEK2, DLGAP5, MELK, CEP55 and CKS2 were upregulated, and ALDH1A1 and KDR were downregulated. Table 4. Expression Of The Ten Selected Hub Genes In Five Oncomine Datasets Gene Name Welsh Ovarian, 2001[173]^35 Lu Ovarian, 2004[174]^36 Adib Ovarian, 2004[175]^37 Hendrix Ovarian, 2006[176]^38 TCGA Ovarian, 2013[177]^39 Fold Change P-value Fold Change P-value Fold Change P-value Fold Change P-value Fold Change P-value EPCAM 7.389 2.61E−09 7.44 4.25E−04 14.692 6.86E−05 2.867 4.44E−04 1.365 2.43E−04 ALDH1A1 −6.913 6.46E−12 −5.273 2.46E−04 −6.999 5.31E−04 −2.523 8.75E−13 −11.737 8.41E−08 ZWINT NA NA 2.662 1.61E−07 6.237 7.54E−04 1.541 9.37E−04 7.001 1.71E−05 BUB1B NA NA 2.312 3.07E−06 3.208 2.00E−03 1.514 2.20E−02 8.04 2.56E−07 NEK2 1.95 2.37E−06 1.186 1.10E−02 1.221 2.70E−02 1.22 3.05E−04 5.887 1.63E−06 DLGAP5 12.717 4.80E−02 1.321 2.00E−03 1.76 3.00E−03 1.479 3.70E−02 8.151 1.01E−06 MELK 23.125 1.90E−02 1.379 1.19E−04 1.331 5.05E−04 1.609 9.00E−03 10.6 2.98E−07 CEP55 NA NA 2.457 9.44E−07 NA NA 1.748 1.20E−02 8.075 1.50E−08 CKS2 13.417 4.20E−02 3.046 3.02E−08 4.612 8.00E−03 1.451 1.00E−02 5.956 3.85E−05 KDR 1.074 1.00E+00 −1.622 7.00E−03 −1.406 1.50E−02 −1.137 3.30E−02 −1.708 2.00E−03 [178]Open in a new tab Abbreviations: TCGA, the cancer genome atlas; NA, not available. Figure 8. [179]Figure 8 [180]Open in a new tab Heat map of the ten hub genes in the Oncomine datasets. Each column represents one dataset and each row represents one gene. The number in each rectangle is the log[2]FC value. The color gradient from blue to red represents the progression from down- to upregulation. Abbreviation: FC, fold change. We performed qRT-PCR experiments to validate the expression of these ten hub genes in 22 HGSOC samples and 22 normal tissues. Consistent with the Oncomine data, the mRNA expression of EPCAM, ZWINT, BUB1B, NEK2, DLGAP5, MELK, and CKS2 was significantly higher in HGSOC than in normal ovarian tissue, and the mRNA expression of ALDH1A1 and KDR was lower in HGSOC tissues, as mentioned above ([181]Figure 9). Although the mRNA expression level of CEP55 was elevated in HGSOC tissues compared with normal tissues, the difference was not statistically significant ([182]Figure 9). Figure 9. [183]Figure 9 [184]Open in a new tab qRT-PCR analysis of the ten hub genes in 22 HGSOC and 22 normal ovarian samples. *P < 0.05, **P < 0.01, ***P < 0.001 and ****P < 0.0001. Abbreviations: qRT-PCR, quantitative real-time PCR; HGSOC, high-grade serous ovarian cancer. Additionally, we investigated EPCAM, ZWINT, DLGAP5, and KDR protein expression in clinical tissue specimens collected from patients. As shown in [185]Figure 10, Western blotting confirmed the overexpression of both EPCAM and ZWINT and the downregulation of KDR in HGSOC tissues (n= 22) compared with normal ovarian tissues (n= 22). These results were consistent with the qRT-PCR results. However, DLGAP5 expression levels were not as expected. Although it is possible that complex physiological regulatory mechanisms lead to discrepant protein and RNA levels, more clinical samples are needed to verify the results for this gene. Figure 10. [186]Figure 10 [187]Open in a new tab The protein expression levels of EPCAM (A), ZWINT (B), DLGAP5 (C) and KDR (D) in 22 HGSOC and 22 normal ovarian samples were analyzed by Western blotting. Signal intensities were quantified by ImageJ and normalized to an internal control (β-actin). *P < 0.05, **P < 0.01, ***P < 0.001 and ****P < 0.0001. Abbreviation: HGSOC, high-grade serous ovarian cancer. Discussion EOC is a highly lethal malignancy due to the difficulty of early diagnosis and the lack of effective treatments for advanced-stage disease. HGSOC, the main type of EOC, accounts for the majority of deaths. Therefore, it is essential to identify key genes associated with the early diagnosis and prognosis of HGSOC to improve the survival rate. Recently, integrated bioinformatics approaches based on microarray technology have been widely used to discover promising biomarkers for cancer diagnosis, treatment and prognosis. In this study, we identified candidate biomarkers for HGSOC via integrated bioinformatic analysis. Six microarray datasets from GEO were integrated, and 103 DEGs between HGSOC and normal OSE samples were successfully identified, comprising 28 upregulated and 75 downregulated genes. Functional enrichment analysis showed that the upregulated genes were enriched in cell division and cell proliferation, which are critical for tumor growth, and that the downregulated genes mainly participated in various metabolic processes, including retinol metabolism, tyrosine metabolism, and drug metabolism/cytochrome P450 pathways, as well as the Wnt signaling pathway, which contributes to carcinogenesis.[188]^40^–[189]^42 A large number of studies have proposed that disorders in metabolism, including retinol metabolism[190]^43^,[191]^44 and tyrosine metabolism,[192]^45^,[193]^46 play important roles in carcinogenesis. PPI network analysis enabled screening of the top 10 hub nodes according to the degree of connectivity, revealing eight upregulated genes (EPCAM, ZWINT, BUB1B, NEK2, DLGAP5, MELK, CEP55 and CKS2) and two downregulated genes (ALDH1A1 and KDR) that might be important in HGSOC pathogenesis. These ten genes can be used in the early diagnosis of patients with HGSOC. The expression levels of these ten hub genes were validated in five Oncomine datasets. In addition, we performed module analysis by using the MCODE app. Seven of the selected hub nodes, namely, ZWINT, BUB1B, NEK2, DLGAP5, MELK, CEP55 and CKS2, partitioned into module 1; coincidentally, these nodes contained genes upregulated in HGSOC that are closely associated with cell cycle and DNA replication. Survival analysis revealed that MELK, CEP55 and KDR expression levels were significantly correlated with the OS rates of patients with HGSOC. MELK overexpression was associated with prolonged OS for HGSOC patients, whereas increased CEP55 and KDR expression levels in HGSOC tissues were significantly associated with decreased OS rates for patients with HGSOC. Similar results were found in some prespecified subgroup analyses restricted to patients at different stages or grades or based on TP53 mutations and chemotherapy regimens. These analyses indicated that MELK, CEP55 and KDR are potential predictors of HGSOC prognosis. Integrated bioinformatics methods are efficient tools for discovering promising biomarkers for cancer diagnosis, treatment and prognosis. However, further biological experiments are needed to validate these data analysis results. In this study, the significant expression levels of nine hub genes, namely, EPCAM, ALDH1A1, ZWINT, BUB1B, NEK2, DLGAP5, MELK, CKS2 and KDR, were confirmed by qRT-PCR at the RNA level. The Western blotting experiments showed that EPCAM, ZWINT, and KDR protein expression levels were consistent with the database information; only the DLGAP5 expression level seemed contrary to expectations. Overall, our experimental results verified the reliability of this research approach. The current study identified ten hub genes that might play important roles in HGSOC pathogenesis, and three hub genes were significantly associated with HGSOC prognosis. These genes have been widely studied in many other types of cancer. MELK is a serine/threonine kinase in the Snf1 (sucrose nonfermenting 1)/AMPK family of kinases.[194]^47 Multiple studies have reported that MELK is overexpressed and plays vital roles in various cancer types, including ovarian cancer,[195]^48 colorectal cancer,[196]^49 breast cancer,[197]^50^,[198]^51 small cell lung cancer,[199]^52 brain cancer,[200]^53 pancreatic cancer,[201]^54 prostate cancer,[202]^55 gastric cancer,[203]^56^,[204]^57 hepatocellular carcinoma,[205]^58 and melanoma.[206]^59 Pitner et al discussed the roles of MELK in cancer, including its ability to mediate proliferation, apoptosis, cancer stem cell phenotypes, epithelial-to-mesenchymal transition (EMT), metastasis, and therapeutic resistance, characterizing MELK as a promising therapeutic target.[207]^60 In addition, Kohler et al showed that MELK overexpression was associated with histological grade (P < 0.05) and progression-free survival (HR = 5.73, P < 0.01) in patients with ovarian cancer;[208]^48 notably, MELK was shown to participate in cell proliferation and tumor growth via cell cycle arrest at the G[2]-M phase.[209]^48 Gray et al reported the same results.[210]^49 Multiple studies have reported that MELK overexpression in patients with cancer is correlated not only with poor prognosis[211]^50^,[212]^53^,[213]^58 but also with chemoresistance and radioresistance.[214]^61^,[215]^62 Based on the above findings, MELK inhibition may be a novel strategy for the treatment aggressive malignancies and combining future MELK inhibitors with chemotherapy and/or radiotherapy may enhance the therapeutic effect of these traditional approaches. Jeffery et al comprehensively discussed the roles of CEP55 in regulating the PI3K/AKT pathway and stemness and in promoting tumorigenesis.[216]^63 CEP55 overexpression is common in cancer tissues and is correlated with an unfavorable prognosis in many human cancers, including non-small cell lung carcinoma,[217]^64 pancreatic cancer,[218]^65 osteosarcoma,[219]^66 and ER+ breast cancer.[220]^67 Shiraishi et al identified CEP55 as a marker to differentiate patients with prostate cancer recurrence following radical prostatectomy.[221]^68 Tao et al discovered that CEP55 expression was strongly elevated in gastric cancer and that CEP55 overexpression promoted the proliferation, colony formation and tumorigenesis of gastric cancer cells.[222]^69 KDR, one of the two VEGF receptors, is also known as vascular endothelial growth factor receptor 2. It plays an essential role in regulating VEGF-induced endothelial proliferation, migration and sprouting and thus promoting angiogenesis, which is essential for cancer growth and metastasis. Takahashi et al showed that VEGF and KDR were overexpressed in metastatic colon cancer, which correlated with the vascularity, metastasis, and proliferation of human colon cancer.[223]^70 An et al found that KDR mRNA expression levels were significantly higher in normal tissues than in lung cancer tissues. In addition, higher KDR expression levels in lung cancer tissues were associated with a shorter survival.[224]^71 The prognostic value of KDR was also identified in patients with stage I non-small cell lung cancer.[225]^72 The authors suggested that patients with KDR-positive tumors have a poor prognosis. Many studies have also reported that SNPs in KDR are associated with survival in patients with colorectal cancer.[226]^73^–[227]^75 EPCAM is an epithelium-specific intercellular adhesion molecule that mediates Ca^2+-independent homophilic cell-cell adhesion.[228]^76 EPCAM has been widely studied as a therapeutic target and a prognostic marker in various epithelial malignancies, including breast, ovarian, urothelial, and non-small cell lung carcinomas. Tayama et al found that high EPCAM levels correlated with chemotherapy resistance and poor prognosis in ovarian cancer patients; thus, targeting EPCAM is a potential promising approach to treat chemoresistant ovarian cancer.[229]^77 However, data from Woopen et al showed that EPCAM overexpression in EOC was significantly associated with improved OS (P = 0.015) and an increased rate of response to platinum-based chemotherapy (P = 0.048).[230]^78 Therefore, future large-scale and well-designed studies are warranted to validate these contradictory findings. In addition, Tomita et al demonstrated that EPCAM gene deletion attributed to Lynch syndrome (also known as hereditary nonpolyposis colorectal cancer, HNPCC) with MSH2 defects.[231]^79 Notably, women with Lynch syndrome have a high risk of ovarian cancer. Liu et al reported that elevated EPCAM expression is correlated with metastasis and cell adhesion in breast cancer tissue.[232]^80 In addition to being observed in epithelial malignancies, EPCAM overexpression is associated with larger tumor size, lymph node metastasis and a poorer prognosis in gastric cancer.[233]^81 However, Wen et al provided evidence strengthening the role of EPCAM in endometrial carcinoma progression and prognosis and reported that EPCAM overexpression favored survival.[234]^82 Furthermore, EPCAM was found to regulate EMT, which is important for cancer metastasis.[235]^83^,[236]^84 The ZWINT protein is involved in kinetochore function and cell proliferation and is thought to be an important regulatory protein for chromosome movement and mitotic spindle checkpoints.[237]^85^,[238]^86 Consistent with our results, Xu et al identified ZWINT overexpression in ovarian cancer by integrated bioinformatic analysis and reported its association with reduced patient OS.[239]^87 Urbanucci et al demonstrated that ZWINT expression was increased in patients with castration-resistant prostate cancer.[240]^88 Moreover, ZWINT expression has been proven to be significantly upregulated and associated with tumor progression and poor prognosis in various cancers, including hepatocellular carcinoma[241]^86 and lung cancer.[242]^85 In contrast, a study by Yang et al indicated that ZWINT expression is downregulated in hepatocellular carcinoma and that patients with low ZWINT expression have a shorter OS and time to recurrence than patients with high ZWINT expression.[243]^89 Further studies are required to define the detailed roles of ZWINT in cancer. NEK2 is a centrosomal serine/threonine kinase that plays a critical regulatory role in mitosis. Recent studies have shown pivotal roles of NEK2 in the development of several cancers. Wu et al showed that NEK2 overexpression promoted liver cancer cell growth, metastasis and angiogenesis.[244]^90 NEK2 induced chemotherapeutic resistance in patients with multiple myeloma,[245]^91 hepatocellular carcinoma,[246]^92 and nasopharyngeal carcinoma.[247]^93 Many other reports have shown that high NEK2 expression predicts poor prognosis in various cancer types.[248]^94^–[249]^97 DLGAP5, also known as HURP and DLG7 (disks large homolog 7), is a cell cycle-related protein that controls microtubule organization and is required for formation of the bipolar spindle.[250]^98^,[251]^99 Chen et al identified DLGAP5 as the direct target gene of NOTCH3 in ovarian cancer progression.[252]^100 Many studies have suggested that DLGAP5 overexpression plays a role in carcinogenesis,[253]^101 for example, in prostate cancer,[254]^102 lung cancer,[255]^103 gastric cancer,[256]^104 pancreatic carcinoma,[257]^105 and hepatocellular carcinoma.[258]^106 Moreover, elevated DLGAP5 expression levels were strongly associated with poor survival in patients with non-small cell lung cancer,[259]^107^,[260]^108 breast cancer,[261]^109 and adrenocortical carcinoma.[262]^110 Moreover, Gomez et al found that DLGAP5 was significantly overexpressed in prostate cancer and that higher DLGAP5 transcript levels were associated with advanced tumor progression and worse prognosis in this disease.[263]^111 Thus, DLGAP5 might be a novel biomarker of tumor progression and prognosis, but additional studies are needed to clarify its role in cancer. CKS2, a cyclin-dependent kinase-interacting protein, is also a cell cycle regulatory protein. Studies have demonstrated that CKS2 expression is elevated in multiple types of cancer, including gastric cancer,[264]^112^,[265]^113 prostate cancer,[266]^114 cholangiocarcinoma,[267]^115 esophageal carcinoma,[268]^116 and liver cancer.[269]^117 CKS2 overexpression contributes to tumor progression and predicts an unfavorable prognosis. Kita et al found that CKS2 protein expression was significantly correlated with several clinicopathologic parameters, including the depth of tumor invasion, clinical stage and poor five-year survival in esophageal squamous cell carcinoma.[270]^118 Yu et al reported that CKS2 was upregulated in colorectal cancer tissues and that CKS2 expression levels were significantly associated with tumor size and tumor stage. The author and his colleagues also observed that higher CKS2 expression predicted a poor outcome.[271]^119 Conclusion In conclusion, by performing an integrated bioinformatic analysis of six GEO datasets, we identified ten hub genes that might be involved in HGSOC pathogenesis. Module analysis placed seven of the selected hub nodes in module 1, in which cell cycle and DNA replication were identified as the crucial pathways enriched by the hub genes. Additionally, survival analysis revealed that three hub genes are associated with the OS rate in patients with HGSOC. Taken together, our findings reveal that MELK, CEP55 and KDR are potential biomarkers that could facilitate the early diagnosis and treatment of HGSOC and predict the prognosis of patients with HGSOC. However, further molecular biology experiments are needed to verify our findings and confirm the potential clinical value of these genes as biomarkers. Acknowledgments