Abstract Ovarian epithelial cancer (OEC) is an often fatal disease with poor prognosis in women with high-stage disease. In contrast, ovarian low malignant potential (LMP) tumors with favorable prognosis behaves as a disease between benign and malignant tumors. The involved genes and pathways between benign-like LMP and aggressive OEC are largely unknown. This study integrated two cohorts profile datasets to investigate the potential key candidate genes and pathways associated with OEC. Gene expression in two datasets ([29]GSE9891 and [30]GSE12172), including 327 OECs and 48 LMP tumors, were analyzed. 559 differentially expressed genes were found to overlap, 251 up-regulated and 308 down-regulated. Subsequently, analysis of gene ontology, signaling pathway enrichment and protein-protein interaction (PPI) network was performed. Gene ontology analysis clustered the up-regulated and down-regulated genes based on significant enrichment. 282 nodes/ differentially expressed genes (DEGs) were identified from DEGs PPI network complex, and two most significant k-clique modules were identified from PPI. In a summary, using integrated bioinformatics analysis, we are able to identify biomarkers potentially significant in the pathogenesis of OEC, which can improve our understanding of the cause and molecular events. These candidate genes and pathways could be used for further confirmation, and lead to better disease diagnose and therapy. Keywords: ovarian epithelial cancer, low malignant potential tumor, chemokine, neoplasm invasiveness, bioinformatics Introduction Ovarian epithelial cancer (OEC) is the fourth leading cause of female cancer death in the developed world [31]^1, with 50% of all cases occurring in women older than 65 years [32]^2. Ovarian epithelial cancers account for 75% of all ovarian tumors, and 90-95% of ovarian malignancies. The American Cancer Society estimates a total of 22,240 new cases and attributes 14,070 deaths to the disease in 2018. Ovarian cancer burden in China is relatively stable due to the increased aging population [33]^3. The disease has a low survival rate due to the fact that the majority of patients when first diagnosed are already at an advanced stage of ovarian cancer since symptoms are often not apparent during the early development of it. Surgery followed with platinum-based chemotherapy is the standard treatment. In spite of the great improvement in the current therapeutic approach, the survival rate at 5 years for the whole population of ovarian cancer patients is still low (46.5%), and is even worse (29%) in women with late-stage distant disease diffusion [34]^4. High mortality in ovarian cancer underscores the high demand to reveal underlying molecular mechanisms, and to discover molecular biomarkers for early diagnosis, prevention and target therapy. Different from ovarian epithelial cancer, low malignant potential (LMP) tumors are a distinct subset of epithelial tumor with behavior characteristics in between benign and malignant tumors, and make up approximately 15% of all epithelial ovarian tumors [35]^5. As abnormal cells form in the tissue covering the ovary and do not usually grow into the stroma, LMP tumors are also called borderline malignant ovarian cancer. LMP tumors, which most often affect younger women, display atypical nuclear structure and metastatic behavior, but are considered noninvasive with 5-year survival rates greater than 95% in contrast to a <45% survival for advanced high-grade OEC [36]^6^, [37]^7. In the past two decades, about twenty microarray-based gene profiles were performed to seek new insights for biomarkers of ovarian cancers [38]^8. They focused on the prognostic value of gene expression signatures and made advances in disease stratification and prognosis prediction. However, the results were generated from an individual cohort, and were always inconsistent in independent studies. The genetic mechanisms of OEC are far from being understood. An integrated bioinformatics method combined with gene expression profiling will be important to discover more reliable biomarkers for ovarian epithelial cancer. In this study, the screening of differentially expressed genes between ovarian low malignant potential tumors and ovarian malignant tumors in two individual gene expression datasets were performed, and subsequent analysis of gene ontology, signaling pathway enrichment and protein-protein interaction (PPI) network were carried out. We proposed novel biomarkers for further study. Materials and Methods Gene Datasets Two gene expression datasets ([39]GSE9891, [40]GSE12172) were obtained from NCBI Gene Expression Omnibus (GEO) database, available at [41]http://www.ncbi.nlm.nih.gov/geo/. [42]GSE9891 comprises molecular profiling from 285 ovarian samples, which include 18 ovarian low malignant potential tumors and 267 malignant ovarian cancers. [43]GSE12172 includes the expression profile of 30 low malignant potential tumors and 60 ovarian malignant serous tumors that originated from ovary epithelial tissue. The experiments of both datasets were performed in Affymetrix Human Genome U133 Plus 2.0 Array microarray platform. Data Preprocessing and Differentially Expressed Genes (DEGs) Analysis By using the robust multi array average (RMA) [44]^9 algorithm in R affy package, the raw array data was converted into expression values, and subsequently background correction, quintile normalization and probe summarization were performed. Differentially expressed genes between low malignant potential tumors and malignant ovarian cancers were analyzed by paired t-test based on the limma package in R language. The adjustment of raw p-value to false discovery rate (FDR) was carried out by the Benjamini & Hochberg method [45]^10. FDR < 0.01 and |log[2]FC| > 1.5 were considered as the cutoff value for DEGs screening. Gene ontology and signaling pathway enrichment analysis of DEGs Gene ontology analysis (GO) and functional enrichment of the DEGs in the molecular function, biological process and cellular component categories were performed with DAVID, Panther ([46]http://www.pantherdb.org/) [47]^11 and GO ([48]http://geneontology.org/) [49]^12 online database. KEGG PATHWAY ([50]http://www.genome.jp/kegg) and Reactome ([51]http://www.reactome.org) databases were used to perform signaling pathway enrichment analysis with p < 0.05 as a cut-off criterion. Construction of the PPI network and module analysis The interactions between the proteins translated from the identified DEGs were searched by STRING Database ([52]http://www.string-db.org/, version 10.5) [53]^13, and a confidence score > 0.4 was used as cut-off criterion. Then cytoscape software ([54]http://www.cytoscape.org/) [55]^14 was used to visualize the PPI network. Cluster analysis of the PPI network was performed by CFinder ([56]http://www.cfinder.org/) [57]^15. CFinder is based on the Clique Percolation Method algorithm to locate the k-clique communities of networks, where k refers to the number of nodes in the subgraph [58]^16. k is the size of the complete subgraphs whose large scale organizations are analytically and numerically investigated. A k-cliques value of > 10 was selected as the cut-off criterion. Pathway enrichment analysis of two selected modules was performed with a cut-off of p < 0.05. Results Identification of DEGs 795 DEGs were screened in [59]GSE9891 database, which include 404 up-regulated genes and 391 down-regulated ones.724 DEGs were screened in [60]GSE12172 database, which include 340 up-regulated genes and 384 down-regulated ones. Among these two databases, there are total 559 genes (251 up-regulated and 308 down-regulated) that overlapped (Figure [61]1A, B). Figure 1. [62]Figure 1 [63]Open in a new tab Identification of 559 commonly changes DEGs from the two profile datasets ([64]GSE9891, [65]GSE12172). Different color areas represented different datasets. The cross areas meant the commonly changed DEGs. DEGs were identified with t-test, statistically significant DEGs were defined with p<0.05 and [log[2]FC]>1.5 as the cut-off criterion. DEGs Gene Ontology Analysis Three groups of DEGs: molecular function, biological process and cellular component were classified by gene ontology (GO) analysis (Figure [66]2). Among these GO functions, cellular process (GO:0009987), cell part (GO:0044464), binding (GO:0005488), metabolic process (GO:0008152) and biological regulation (GO:0065007) were the top five ones involved in ovarian cancer. Up-regulated genes were mainly enriched in CXCR3 chemokine receptor binding, immune response, chemokine-mediated signaling pathway, chemokine activity, positive regulation of cAMP metabolic process, transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding and sequence-specific DNA binding (Table [67]1). Down-regulated genes were mainly enriched in microtubule, microtubule motor activity, dynein complex, nucleoside kinase activity and nucleoside diphosphate kinase activity (Table [68]2). Figure 2. [69]Figure 2 [70]Open in a new tab Gene Ontology analysis of DEGs in ovarian cancer. GO analysis classified the DEGs into 3 groups: molecular function, biological process and cellular component. Table 1. Significant enriched GO terms of up-regulated DEGs in ovarian cancer Term name Description Count p-value Enrichment Score: 2.46 GO:0048248 CXCR3 chemokine receptor binding 4 1.49E-05 GO:0006955 immune response 17 2.97E-05 GO:0070098 chemokine-mediated signaling pathway 7 1.64E-04 GO:0008009 chemokine activity 6 2.47E-04 GO:0030816 positive regulation of cAMP metabolic process 3 0.00191 GO:0043950 positive regulation of cAMP-mediated signaling 3 0.00804 GO:0007267 cell-cell signaling 9 0.00918 GO:0006935 chemotaxis 6 0.01338 GO:0002690 positive regulation of leukocyte chemotaxis 3 0.01781 GO:0051281 positive regulation of release of sequestered calcium ion into cytosol 3 0.03821 GO:0032496 response to lipopolysaccharide 6 0.04135 Enrichment Score: 2.06 GO:0001077 transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding 12 9.65E-05 GO:0043565 sequence-specific DNA binding 17 3.48E-04 GO:0045944 positive regulation of transcription from RNA polymerase II promoter 22 0.00447 GO:0000978 RNA polymerase II core promoter proximal region sequence-specific DNA binding 11 0.00836 GO:0006366 transcription from RNA polymerase II promoter 13 0.01569 GO:0003700 transcription factor activity, sequence-specific DNA binding 17 0.09053 [71]Open in a new tab Table 2. Significant enriched GO terms of down-regulated DEGs in ovarian cancer Term name Description Count p-value Enrichment Score: 3.67 GO:0005874 microtubule 17 1.12E-06 GO:0003777 microtubule motor activity 9 3.62E-06 GO:0030286 dynein complex 6 5.22E-06 GO:0007018 microtubule-based movement 7 3.48E-04 GO:0016887 ATPase activity 7 0.01885 Enrichment Score: 1.80 GO:0019206 nucleoside kinase activity 3 0.00550 GO:0004550 nucleoside diphosphate kinase activity 3 0.02158 GO:0015949 nucleobase-containing small molecule interconversion 3 0.03346 [72]Open in a new tab Signaling Pathway Enrichment Analysis Based on KEGG and Reactome databases, the top five enriched pathways of up-regulated genes were related to Chemokine receptors bind chemokines, G alpha (i) signalling events, Chemokine signaling pathway, Cytokine-cytokine receptor interaction and Mitotic Prometaphase (Table [73]3). ACKR2, CCL11, CDCA8, CENPE, CXCL10, CXCL11, CXCL13, CXCL14, CXCL9, GNGT1, GPR17, NPW, NUF2, PDYN, PTGER3, SKA1, SPC25 and TNFRSF17 were the genes involved in these five signaling pathways. The most significant enriched pathway in down-regulated genes was Interconversion of nucleotide di- and triphosphates, involving AK7, AK9 and AK8 (Table [74]3). Table 3. Signaling pathway enrichment analysis of DEGs function in ovarian cancer Pathway Name Count p-value Genes Up-regulated DEG R-HSA-380108 Chemokine receptors bind chemokines 6 5.26E-04 CCL11, CXCL13, CXCL9, ACKR2, CXCL11, CXCL10 R-HSA-418594 G alpha (i) signaling events 9 0.00675 GNGT1, PTGER3, CXCL13, NPW, CXCL9, GPR17, PDYN, CXCL11, CXCL10 KEGG: hsa04062 Chemokine signaling pathway 7 0.01032 CCL11, GNGT1, CXCL14, CXCL13, CXCL9, CXCL11, CXCL10 KEGG: hsa04060 Cytokine-cytokine receptor interaction 7 0.02668 CCL11, CXCL14, CXCL13, CXCL9, TNFRSF17, CXCL11, CXCL10 R-HSA-68877 Mitotic Prometaphase 5 0.03419 SPC25, CDCA8, NUF2, CENPE, SKA1 R-HSA-2500257 Resolution of Sister Chromatid Cohesion 5 0.04855 SPC25, CDCA8, NUF2, CENPE, SKA1 Down-regulated DEG R-HSA-499943 Interconversion of nucleotide di- and triphosphates 3 0.01646 AK7, AK9, AK8 [75]Open in a new tab Analysis of PPI network and modules A total of 282 DEGs of the 559 commonly changed DEGs were screened into the DEGs PPI network complex, containing 282 nodes and 652 edges (Figure [76]3). Using CFinder with a k-cliques value of > 10, two modules, including module 1 (Figure [77]4A) and module 2 (Figure [78]4B) were extracted from the constructed PPI network. Pathway enrichment analysis showed that Module 1 consisted of 12 nodes and 66 edges (Figure [79]4A, Table [80]4), which are mainly associated with G-protein coupled receptor signaling pathway, CXCR3 chemokine receptor binding and CXC chemokine, and that Module 2 consisted of 10 nodes and 45 edges (Figure [81]4B, Table [82]5), which are mainly associated with mitotic cell cycle process, nuclear division and chromosome segregation. Figure 3. [83]Figure 3 [84]Open in a new tab DEGs protein-protein interaction network complex. Total of 282 DEGs were screened into, containing 282 nodes and 652 edges. Figure 4. [85]Figure 4 [86]Open in a new tab Two most significant k-clique modules in the PPI network. (A) Module 1 consisted of 12 nodes and 66 edges, which are mainly associated with G-protein coupled receptor signaling pathway, CXCR3 chemokine receptor binding and CXC chemokine; (B) Module 2 consisted of 10 nodes and 45 edges, which are mainly associated with mitotic cell cycle process, nuclear division and chromosome segregation. Genes in red represent as upregulation and genes in green represent as downregulation. Table 4. Top 10 of pathway enrichment analysis of Module 1 genes function Category Term name Description Count p-value GO Process GO:0007186 G-protein coupled receptor signaling pathway 11 7.03E-10 GO Function GO:0048248 CXCR3 chemokine receptor binding 4 1.40E-09 InterPro IPR001089 CXC chemokine 4 2.73E-08 InterPro IPR018048 CXC chemokine, conserved site 4 2.73E-08 GO Process GO:0070098 chemokine-mediated signaling pathway 5 1.36E-07 GO Process GO:0048247 lymphocyte chemotaxis 4 9.50E-07 InterPro IPR001811 Chemokine interleukin-8-like domain 4 1.05E-06 Pfam PF00048 Small cytokines (intecrine/chemokine), interleukin-8 like 4 1.48E-06 GO Process GO:0032496 response to lipopolysaccharide 6 9.96E-06 KEGG Pathways 04062 Chemokine signaling pathway 5 1.02E-05 [87]Open in a new tab Table 5. Top 10 of pathway enrichment analysis of Module 2 genes function Category Term name Description Count p-value GO Process GO:1903047 mitotic cell cycle process 9 3.82E-09 GO Process GO:0007049 cell cycle 10 3.82E-09 GO Process GO:0000278 mitotic cell cycle 9 6.37E-09 GO Process GO:0007067 mitotic nuclear division 7 7.95E-08 GO Process GO:0007059 chromosome segregation 6 2.14E-07 GO Process GO:0051301 cell division 7 4.07E-07 GO Component GO:0000793 condensed chromosome 5 8.13E-06 GO Component GO:0000775 chromosome, centromeric region 5 8.13E-06 GO Component GO:0000777 condensed chromosome kinetochore 4 4.33E-05 GO Component GO:0000779 condensed chromosome, centromeric region 4 4.33E-05 [88]Open in a new tab Discussion Invasiveness is one of the aggressive features of ovarian epithelial cancer, in particular the advanced high-grade disease. In contrast, the low malignant potential tumor behaves as a located lesion within the tissue covering the ovary. Characterizing the molecular difference between LMP tumor and OEC will shed light on the mechanism controlling the invasive trait of OEC malignances. In the present study, we integrated two individual cohorts of profile datasets, and used multiple bioinformatics tools to identify significant genes and pathways between LMP tumors and ovarian epithelia cancers. A total of 559 differentially expressed genes were screened in two public available GEO datasets, including 251 up-regulated and 308 down-regulated genes. Based on the gene ontology analysis, a greater majority of the DEGs were involved in biological process, specifically the cellular process, than in cellular component or molecular function (Figure [89]2). GO function analysis revealed that more significant enrichments in the up-regulated DEGs were CXCR3 chemokine, immune response and chemokine-mediated signaling pathway. Through literature survey, we found that the increased expression and release of pro-inflammatory chemokines have been associated with progression of ovarian cancer, which may induce tumor cell proliferation, survival, migration, and angiogenesis [90]^17^-[91]^19. The enrichments of chemokine and immune response pathway in this study highly validated their important roles in OEC malignancies. Consistent with Gene Ontology analysis, signaling pathway enrichment analysis displayed the same signaling pathways in the up-regulated DEGs: Chemokine receptors bind chemokines, G alpha (i) signaling events, Chemokine signaling pathway and Cytokine-cytokine receptor interaction. Chemokines are a family of cytokines that induce chemotaxis of target cells and bind to the G protein coupled chemokine receptors. Other than the function of inducing leukocyte migration (including dendritic cells, macrophages, and neutrophils) into the infected or injured sites [92]^20, they can promote cancer progression [93]^21^-[94]^25. GO analysis in down-regulated DEGs showed that more significant enrichments were microtubule, microtubule motor activity, dynein complex and microtubule-based movement. All of these enrichments were associated with microtubule function. Microtubules are highly dynamic structures that play an important role in cellular growth, vesicular transport and mitosis [95]^26. They are composed of α/β-tubulin heterodimers [96]^27. Alterations in specific β-tubulin isotypes in epithelial cancers are associated with resistance to tubulin-binding agent chemotherapy and more aggressive disease [97]^28. Microtubules were linked in tumor cell migration and metastasis [98]^26. More interestingly, 12 genes were selected in the most significant k-clique module 1 of the PPI network: APLNR, CNR1, CXCL9, CXCL10, CXCL11, CXCL13, GNGT1, GPR17, GRM7, NPW, PDYN, PTGER3. Among them, most genes were chemokines and involved cancer immunity. APLNR is a member of G-protein coupled receptor signaling pathway, which was significantly mutated in human cancers [99]^29. APELA, the ligand of APLNR, can promote ovarian cancer cell growth and migration [100]^30. CNR1 gene expression was down-regulated in endometrial carcinomas, another gynecologic cancer which was validated by qRT-PCR [101]^31. Inversely, Messalli et al. found that expression of CB1R (alias of CNR1) increased from benign and borderline to malignant tumors by immunohistochemical quantification [102]^32. A possible reason for the discrepancy between gene and protein expression is that other post-transcription and/or post-translation mechanisms affect CNR1 expression. CXCL9 and CXCL10 expression was associated with improved patient survival, and these two chemokines were synergistically induced by inflammatory cytokines [103]^33. CXCL11 promoted proliferation and migration of ovarian cancer cells via the chemokine receptor CXCR3, thus CXCL11-CXCR3 signaling represented therapeutic targets in ovarian cancer [104]^34. In independent cohorts of ovarian cancer patients, high CXCL13 correlated strongly with better prognosis [105]^35. Based on the function of uncovered genes, we hypothesized, GNGT1, GPR17, GRM7, NPW, PDYN and PTGER3 can influence OEC development via chemokine signaling. In module 2, ASPM, CDCA8, CENPE, EXO1, FOXM1, MCM10, NCAPH, NUF2, SKA1 and SPC25 were displayed and mainly associated with mitotic cell cycle. Most of these genes were studied and identified in ovarian cancers. Deregulation of ASPM in OEC correlates with tumor progression, grade and survival [106]^36^, [107]^37. CDCA8 was significantly higher in ovarian cancer cells compared with ovarian epithelial cells by using quantitative PCR with reverse transcription analysis [108]^38. CENPE, a cell cycle regulating gene, was up-regulated in chemo-resistant ovary tumors [109]^39. CENPE proteins were significantly up-regulated in the fibroblasts co-cultured with ovarian cancer cells [110]^40. Attenuating EXO1 expression by small interfering RNA augments the chemotherapy efficacy against ovarian cancer [111]^41. FOXM1 acted as a transcriptional activator involved in cell proliferation to promote cell cycle progression in OEC cells [112]^42. NUF2 was significantly aberrantly overexpressed in ovary serous adenocarcinomas, and silencing of NUF2 induced increased apoptosis [113]^43. It is speculated that the undiscovered four genes (MCM10, NCAPH, SKA1, SPC25) involved in cell cycle may exert effects in OEC. Taken together, chemokine-related signaling, mitotic cell cycle and microtubule pathways play important roles in pathogenesis and aggressiveness of ovarian epithelial cancer. Nine novel up-regulated genes (GNGT1, GPR17, NPW, PDYN, PTGER3, MCM10, NCAPH, SKA1, SPC25) and one down-regulated gene (GRM7) are linked with OEC invasiveness. This is the first time the identified genes from this study have been proposed, but these results will require confirmation through further studies. Acknowledgments