Identification of Chemoresistance-Associated Key Genes and Pathways in High-Grade Serous Ovarian Cancer by Bioinformatics Analyses

Introduction

Epithelial ovarian cancer has the highest mortality rate of any gynecological cancer, with patients with high-grade serous ovarian cancer (HGSOC) having particularly poor outcomes.^1,2 Typical treatments for HGSOC include surgical resection in combination with postoperative chemotherapy using cisplatin and paclitaxel.³ However, the vast majority of patients with advanced disease relapse within 5 years, often owing to metastasis and drug resistance of ovarian cancer cells.^4–7 Therefore, there is a need to find more efficient therapeutic targets, including key genes and signaling pathways that drive therapy resistance.

Some progress has been made in determining the mechanisms of chemoresistance in HGSOC.^5,8 For example, Zhang et al⁹ isolated CD44+/CD117+ ovarian cancer cell stem cells (CSCs) and found that they exhibited enhanced chemoresistance to the ovarian cancer chemotherapeutics cisplatin or paclitaxel. Liu et al¹⁰ found that C/EBPβ-mediated reprogramming of gene expression triggered a broad signaling network that synergistically promoted cisplatin resistance in HGSOC. Luo et al¹¹ showed that loss of ARID1A in HGSOC led to multiple drug resistance through the upregulation of MRP2. Recently, tumor metabolism has also been implicated in chemoresistance; chemotherapeutic drugs combined with metabolic targeting appears to be a promising approach to overcoming chemoresistance.^12,13 Nevertheless, these studies are only the tip of the iceberg regarding the mechanisms of platinum resistance in HGSOC. Stronger links between molecular profiles and drug resistance are needed.

Bioinformatics analysis is an effective and practical method to predict key genes and pathways in tumorigenesis or other pathological processes.¹⁴ In this study, we analyzed microarray profiles of 12 platinum-resistant HGSOC samples and 16 platinum-sensitive control samples from the Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) were screened, and gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted. We also constructed a protein–protein interaction (PPI) network and screened out the hub genes involved in the development of platinum-based chemotherapy resistance in HGSOC patients. Finally, we evaluated the significance of these hub genes with respect to prognosis, and chose one marker to validate its expression. This study may provide novel insights into HGSOC chemoresistance and identify potentially important therapeutic targets.

Materials and Methods

Data Sets

Microarray gene expression profiles from the GSE51373 data set¹⁵ were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). This data set is based on the GPL570 Affymetrix Human Genome U133 Plus 2.0 Array platform (Affymetrix Inc., Santa Clara, CA, USA) and contained 28 patients, who were divided into two groups (chemoresistant and chemosensitive groups). The inclusion criterion was as follows: patients with high-grade serous HGSOC, treated with the same standard platinum-based chemotherapy (Carboplatin/paclitaxel). Twelve patients demonstrating relative resistance to platinum chemotherapy corresponding to shorter progression-free survival (PFS < 8 months) were compared with 16 platinum-sensitive patients (PFS > 18 months).

Data Pre‑processing and Differential Expression Analysis

RStudio software (version 1.1.447) and various R packages were used to analyze the original array data. In short, we first assessed the quality of the raw data by plotting the normalized unscaled standard error (NUSE) boxplot and plot of residuals. Then, background correction and quantile normalization were performed on the raw data using the robust multi-array average algorithm in the R affy package. Subsequently, the DEGs between the chemosensitive and chemoresistant group samples were sorted by paired t-tests using the limma package in R. Multiple comparisons were corrected by the Benjamini–Hochberg method to obtain adjusted P-values. Finally, genes with adjusted P<0.05 and |log₂ fold change (FC)|>1 were considered to be significant.

Functional Enrichment Analysis

GO enrichment and KEGG pathway analyses were performed for gene annotation and functional enrichment using the online tool Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/).¹⁶ The resulting GO terms and KEGG pathways with P<0.05 were considered to be significantly enriched in the obtained DEGs.

Construction of PPI Network

To evaluate the interactive relationships among DEGs, we mapped the DEGs to the STRING database (http://string-db.org). Only the interactions with a combined score >0.15 were considered significant. Hub genes were then selected from the PPI network and a score was calculated for each gene (the number of genes directly interacting with it). The top ten hub genes in the network were identified based on their scores.

Evaluation of the Prognostic Significance of Top Ten Hub Genes

Kaplan–Meier Plotter, an online survival analysis tool (http://kmplot.com/analysis/),^17,18 was used to evaluate the prognostic significance of the top ten hub genes in HGSOC. According to their median expression value, the patient samples were divided into high- and low-expression groups. Hazard ratios with 95% confidence intervals and log-rank P-values were calculated.

Immunohistochemistry

Immunohistochemical analyses were performed using the Envision System with diaminobenzidine (Gene Tech Co., Ltd., Shanghai) according to the manufacturer’s protocol. Chemoresistant and chemosensitive samples were collected from Fudan University Shanghai Cancer Center. In brief, specimens were incubated first with an anti-MGP antibody (10734-1-AP; 1:400, Proteintech, China) overnight at 4°C and then with a biotinylated secondary antibody (1:100, goat anti-rabbit IgG) for 30 mins at 37°C.

Ethics Statement

The tissue samples from patients in this study were used with approval from the Ethics Committee of Fudan University Shanghai Cancer Center and with informed consent from all patients. All procedures were performed in accordance with the Declaration of Helsinki and relevant policies in China.

Results

Quality Control of Data Sets

Good quality control is essential for successful microarray data analysis. In the present study, we used NUSE boxplots and residual plots to assess the quality of the data sets. As described in the literature, if a microarray data set is very reliable, its NUSE values will be very close to one. Our results (Figure 1A) show that this was the case for our data set. The residuals were distributed uniformly (shown in red and blue in Figure 1B), which further confirmed the quality of the data set.

Figure 1 Evaluation of the quality of the microarray data set. (A) Normalized unscaled standard error (NUSE) boxplot for all samples. (B) Plot of residuals also confirmed the quality of the microarray data set.

Screening of DEGs

After quality assessment, a set of 20,460 genes were mapped to probes for each CEL file. Then, we defined 109 DEGs between chemoresistant and chemosensitive samples, of which 64 (58.7%) were upregulated and 45 (41.3%) were downregulated (Table S1). Genes with |logFC| >1 and adjusted P<0.05 are marked in red in the volcano plot in Figure 2A. A heat map of these DEGs is presented in Figure 2B, and the top ten DEGs are listed in Table 1.

Table 1 Top Ten Differentially Expressed Genes Between Chemosensitive and Chemoresistant Samples

Figure 2 Expression levels of genes and distributions in microarray data set. (A) Volcano plots of microarray data. Y-axis represents log2 FC. X-axis represents adjusted P-value (chemosensitive vs chemoresistant samples). Red dots indicate DEGs. (B) Heatmap of DEG clustering. Green represents downregulation and red represents upregulation.

Significant Functions and Pathway Enrichment Analysis

Potential biological functions associated with these DEGs may imply intrinsic chemoresistance mechanisms. We used the online software DAVID to identify the representative GO categories and KEGG pathways. The top ten most enriched GO terms according to P-values and numbers of enriched pathways are shown in Figure 3 and Table 2. These GO terms included response to stimulus, cell communication, positive regulation of biological process, and cellular component organization or biogenesis. Many of these terms are closely related to chemoresistance and tumorigenesis. The four KEGG pathways that were enriched are presented in Figure 3 and Table 2. These included the PI3K-AKT signaling pathway and the Rap1 signaling pathway, which may have important roles in chemoresistance in HGSOC patients.

Table 2 Enriched GO Terms and KEGG Pathways for the Identified DEGs

Figure 3 Representative GO categories and KEGG pathways obtained using DAVID.

PPI Network Analysis and Hub Genes Screening

To determine the interactions among the identified genes, we constructed a PPI network using the online tool STRING. Only genes with score >0.15 were included in the network. As shown in Figure 3A, we obtained a network containing a total of 213 PPI relationships and 97 nodes after removing the effects of free protein pairs, accounting for 89% of all DEGs.

In PPI networks, hub genes, which have strong interactions with many other genes, can exert huge effects on tumorigenesis. Owing to their key positions, hub genes are potential drivers of disease. In order to identify the key genes in chemoresistance of HGSOC, the cytoHubba plugin for Cytoscape was used to screen the hub genes. As shown in Figure 3B and Table 3, we obtained ten hub genes (nodes colored red or orange in the figure) according to their scores; these were proto-oncogene c-Kit (KIT), forkhead box protein M1 (FOXM1), fibroblast growth factor 2 (FGF2), H4 clustered histone 4 (HIST1H4D), zinc finger protein FOG family member 2 (ZFPM2), IFN-induced protein with tetratricopeptide repeats 2 (IFIT2), cyclin O (CCNO), matrix Gla protein (MGP), rho-related BTB domain-containing 3 (RHOBTB3), and cell division cycle 7 (CDC7). The hub genes identified in the PPI network analysis may be key players in the development of cancer chemoresistance. Table 4 provides a brief overview of the functions of these hub genes. In addition, we tested the expression of hub genes using a data set from The Cancer Genome Atlas (Figure S1); however, the expression of some hub genes in this data set, including KIT, did not match the results of the current bioinformatics analyses. This was probably due to the small data set itself or the changes in the genome during the development of platinum resistance in ovarian cancer patients.

Table 3 Identified Hub Genes in the PPI Network

Table 4 Overview of the Functions of Top Ten Hub Genes

Evaluation of Hub Genes as Prognostic Markers for HGSOC Patients

Given their vital roles in the process of tumor chemoresistance, we explored the prognostic significance of the hub genes in HGSOC patients. Survival analysis was performed using Kaplan–Meier Plotter, which contains a large number of HGSOC microarray data sets. As shown in Figure 4, we found that five of the ten hub genes could predict the prognosis of HGSOC patients (Figure 5, P<0.05). High expression levels of ZFPM2, MGP, and RHOBTB3 were correlated with decreased overall survival, whereas high expression levels of another two genes, HIST1H4D and CCNO, were associated with positive prognosis.

Figure 4 PPI network of DEGs. (A) A total of 213 PPI relationships and 97 nodes were obtained; PPIs with a score >0.15 were selected. (B) Sub-network of ten hub genes. Color gradient from orange to red indicates increasing score.

Validation of the Expression of MGP in Patients

To further validate the findings, we performed immunostaining with antibodies on chemoresistant (platinum-resistant) and chemosensitive (platinum-sensitive) HGSOC samples collected from our cancer center. Given that no antibodies specific to HIST1H4D, RHOBTB3, CCNO, or ZFPM2 were commercially available, we detected the protein expression of MGP only. Positive immunostaining signals for MGP were observed in the chemoresistant samples (n=3, Figure 6).

Figure 5 Exploration of prognostic value of hub genes. Survival data were analyzed using Kaplan–Meier Plotter, and overall survival curves were plotted based on data of all ovarian cancer patients (n=1656). Red line: patients with expression above the median; black line: patients with expression below the median. Abbreviation: HR, hazard ratio.

Figure 6 Expression of MGP detected in platinum-resistant HGSOC samples compared with platinum-sensitive tissues.

Discussion

Platinum resistance is a major obstacle that limits the efficacy of platinum-based chemotherapy in ovarian cancer.⁷ Therefore, new therapeutic targets are urgently needed. In the present study, the GSE51373 data set was analyzed and all DEGs between chemosensitive and chemoresistant samples were identified. Quality control of the test arrays before analysis confirmed that the samples were of good quality. A total of 109 DEGs were identified, comprising 64 upregulated and 45 downregulated genes.

To obtain a more in-depth understanding of these DEGs, GO analysis was performed to identify related biological processes. The results showed that positive regulation of biological process, cellular response to stimulus, cellular component organization, response to stimulus, and single organism signaling were the top five abundant functions. In addition, we conducted KEGG analysis to investigate the pathways involved. The enriched KEGG pathways for the DEGs included the PI3K-AKT and Rap1 signaling pathways. The PI3K-AKT signaling pathway, mediated by PI3K, is one of the most important signaling cascades, with pivotal roles in various cellular biological processes including survival, growth, proliferation, differentiation, and apoptosis. Dysregulation of this signaling cascade has commonly been observed in a broad range of human cancers.^19–23 A recent study showed that alterations of the PI3K pathway are also part of the chemoresistance profile in ovarian cancer patients.²³ PI3K, a frequently mutated kinase that is known to increase chemoresistance in ovarian cancer, usually acts as an essential lipid second messenger in the translocation of protein kinase B (AKT) to the plasma membrane.²⁴ Downstream of PI3K, AKT is phosphorylated at two sites, T308 in the kinase domain and S473 in the regulatory tail. Phosphorylation at T308 and S473 is essential for maximal AKT activation. Phosphorylated AKT regulates the functions of a broad array of intracellular proteins involved in chemoresistance against platinum in ovarian cancer.²⁵ The Rap1 signaling pathway was also highly enriched. Rap1, which belongs to the Ras family, regulates two important cellular processes: B-RAF/MEK/ERK activation, and integrin-mediated cell adhesion and migration. Many studies have shown that Rap1 and its upstream regulators play an important part in the development of human malignancies.^26–29 However, few studies have investigated its relationship with platinum resistance in ovarian cancer.

A network of interactions between proteins helps to dig up the hub gene, which may be a potential therapeutic target in HGSOC. A series of hub genes in this study were discovered, the top three of which were KIT, FOXM1, and FGF2. KIT, a type III receptor tyrosine kinase, has a crucial role in cancer biology.^30,31 Currently, c-Kit is mainly considered as a stem cell factor. In addition, early studies on c-Kit identified it as an oncogene.^32,33 Chau et al³⁴ first reported a relationship between c-Kit and chemoresistance of ovarian cancer, and showed that c-Kit mediates chemoresistance and the tumor-initiating capacity of ovarian cancer cells through activation of Wnt/β-catenin-ATP-binding cassette G2 signaling. The second hub gene, FOXM1, a transcription factor of the forkhead family, has been demonstrated to be critical for proliferation, apoptosis, migration, and invasion of human cancers.³⁵ FOXM1 has been shown to be overexpressed in cancer cells with resistance to DNA-damaging agents, and its expression can confer resistance to genotoxic agents.^36–38 For example, FOXM1 modulates sensitivity to the DNA-damaging agents cisplatin and doxorubicin through regulating Exo1 in ovarian and breast cancer, respectively.^39,40 The third hub gene, FGF2, a pleiotropic angiogenesis inducer belonging to the heparin-binding FGF growth factor family, is expressed by numerous tumor types and exerts its proangiogenic activity by interacting with tyrosine kinase receptors, heparan-sulfate proteoglycans, and integrins.^41–43 Moreover, accumulating evidence indicates that FGF2 contributes to cytoprotection against drug toxicity via regulating apoptotic mediators and thus has an important role in cancer drug resistance.^44–46 In summary, these hub genes may represent potential therapeutic targets for overcoming platinum resistance in ovarian cancer. However, the role of these genes has not been experimentally confirmed, and there may be some discrepancies owing to sample size, etc. Generally, our results suggest a possible role of some genes in platinum resistance of ovarian cancer from a bioinformatics perspective.

There were several limitations to the current study. First, this study was based on microarray data mining. Post-transcriptional modifications may exist, and protein expression levels may be different from mRNA levels. Therefore, further studies are necessary to confirm these results. Second, the results were only analyzed using bioinformatics; experimental verification is required to confirm our findings. Third, the threshold we set for constructing the PPI network may not have been stringent enough, because we screened out a relatively small number of DEGs. Therefore, PPI networks may not be constructed under a more stringent threshold.

Conclusions

In the present study, we analyzed a data set of microarray profiles from chemoresistant HGSOC patients and control chemosensitive samples using bioinformatics tools. Enriched GO terms and KEGG pathways were identified. Moreover, a PPI network was constructed and network modules were analyzed. Finally, we found that five of the top ten hub genes were associated with the prognosis of HGSOC patients. Collectively, these results may provide a basis for the screening of drug candidates and prognostic factors for HGSOC.