Gene expression profile analysis identifies metastasis and chemoresistance-associated genes in epithelial ovarian carcinoma cells

The purpose of this study was to identify genes that associated with higher ability of metastasis and chemotherapic resistance in epithelial ovarian carcinoma (EOC) cells. An oligonucleotide microarray with probe sets complementary to 41,000+ unique human genes and transcripts was used to determine whether gene expression profile may differentiate three epithelial ovarian cell lines (RMG-I-C, COC1 and HO8910) from their sub-lines (RMG-I-H, COCI/DDP and HO8910/PM) with higher ability of metastasis and chemotherapic resistance. Quantitative real-time PCR and immunohistochemical staining validated the microarray results. Hierarchic cluster analysis of gene expression identified 49 genes that exhibited ≥2.0-fold change and P value ≤0.05. Highly differential expression of GCET2, NLRP4, FOXP1 and SNX29 genes was validated by quantitative PCR in all cell line samples. Finally, FOXP1 was validated at the protein level by immunohistochemistry in paraffin embedded ovarian tissues (i.e., for metastasis, 15 primary EOC and 10 omental metastasis [OM]; for chemoresistance, 13 sensitive and 13 resistant EOC). The identification of higher ability of metastasis and chemotherapic resistance-associated genes may provide a foundation for the development of new type-specific diagnostic strategies and treatment for metastasis and chemotherapic resistance in epithelial ovarian cancer.


Introduction
Epithelial ovarian cancers (EOC) are high aggressive tumors associated with high mortality and morbidity in gynecology. Although the 5-year survival rate is 90 % for women with early-stage ovarian cancer and postoperative introduction of paclitaxel drug have improved the 5-year survival rate for advanced-stage ovarian cancer, patients with this cancer have a 5-year survival rate of only 30 % [1]. Standard therapy includes cytoreductive surgery with first-line combination chemotherapy, 75 % of patients initially respond to conventional chemotherapy; however, 80 % of these women eventually relapse and die from chemotherapy resistant disease. Thus, to understand the molecular basis of epithelial ovarian cancer metastasis and chemotherapeutic resistance is of vital importance and may have the potential to improve significantly the development of more specific and effective treatment against EOC. Comprehensive, high-throughput technologies such as gene expression microarrays have provided powerful tools for this purpose.
In the present study, the whole human genome oligo microarray was used to investigate the differential expression genes (DEGs) in the human ovary cancer cell lines RMG-1-C, COC1, HO8910 and their high malignant and chemoresistant sub-cell lines RMG-I-H, COC1/DDP and HO8910/PM. Hierarchical clustering of genes by the expression level of the DEGs was performed. The potential functions of the DEGs were analyzed by Gene Ontology (GO) and pathway enrichment analyses. In addition, the interaction relationships between these DEGs were investigated by regulatory network. We hope these metastasis and chemotherapy resistance-associated genes may be used for early detection of epithelial ovarian cancer and for development of more specific chemotherapy drugs against EOC.

Cell culture
The human ovarian cancer cell strain RMG-1 was the courtesy of Doctor Iwamori Masao of Kinki University in Japan. We transfected the gene of extrinsic a1,2-fucosyl transferase (a1,2-FT) into RMG-1 to create cell line RMG-1-H with high expression of Lewis (y) and a1,2-fucosyl transferase [2,3], and we discovered that compared with the empty plasmid vector transfected cell line RMG-1-C, cell line RMG-1-H showed enhanced cellular malignant biological behaviors, such as enhanced metastasis and proliferation [4], adhesion [5] and multiple drug resistance [6].
Human ovarian cancer cell lines HO-8910 and HO-8910 PM (a highly metastatic cell line derived from HO-8910) were purchased from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China) [7], human ovarian cancer cell lines COC1 and COC1/DDP (a platinum resistance cell line derived from COC1) were purchased from the China Center for Type Culture Collection (Wuhan, China). Cells were cultured in RPMI-1640 medium supplemented with 100 units/mL penicillin/streptomycin and 10 % fetal bovine serum (FBS) and maintained in an incubator at 37°C under a humidified atmosphere of 5 % CO 2 . COC1/DDP cells were cultured in RPMI-1640 medium containing 0.5 lg/mL cisplatin (Sigma, St. Louis, MO, USA) to maintain the drug resistant phenotype. The cell lines and labels in this experiment are listed in Table 1.   Total RNA extraction and gene chip hybridization   Total RNA was extracted from all 6 cell line samples with  TRIZOL reagent (Life Technologies, Inc, Carlsbad, CA) and further purified with RNeasy Min-elute Clean-up Columns (Qiagen, Valencia, CA), as described by the manufacturers. Optical density for each sample of RNA was measured at OD 260 nm and OD 280 nm using NanoDrop ND-1000 (NanoDrop Technologies, Wilmington, DE, USA). All RNA samples isolated OD260/280 ratio should be close to 2.0 for pure RNA (ratios between 1.8 and 2.1 are acceptable). The OD A260/A230 ratio should be more than 1.8. Each isolated RNA sample was subjected to further quality check to ensure integrity of RNA with Agilent RNA 6000 Nano LabChip using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). All RNA samples were verified to be intact with distinct 28S and 18S RNA bands at a ratio of approximately 2:1 and a RNA integrity number (RIN) [ 7.
The samples were amplified and labeled using the Agilent Quick Amp labeling kit and hybridized with Agilent whole genome oligo microarray in Agilent's SureHyb Hybridization Chambers in accordance with the manufacturer's instructions. This array contains 41,000 ? unique human genes and transcripts represented, all with public domain annotations, content sourced from RefSeq, Goldenpath Ensembl Unigene Human Genome (Build 33) and GenBank databases, over 70 % of the represented probes are validated by Agilent's laboratory validation process, 4 9 44 K slide formats printed using Agilent's 60-mer SurePrint technology. After hybridization and washing, the processed slides were scanned with the Agilent DNA microarray scanner (part number G2505B) using settings recommended by Agilent Technologies.

Data analysis and clustering
Agilent Feature Extraction Software (version 10.5.1.1) was used to extract the signal intensity values from each gene chip, and the resulting text files were imported into the Agilent GeneSpring GX software (version11.0) for further analysis. The 6-microarray data sets were normalized in GeneSpring GX using the Agilent FE one-color scenario (mainly median normalization), and genes marked present or marginal in all samples were chosen for data analysis. DEGs were identified through fold-change screening comparing between cell lines group A, B, C and cell lines group 1, 2 and 3. The threshold used to screen up or downregulated genes is fold change C2.0 and P value B0.05. A scatter plot was made to visualizingly assess the variation Real-time polymerase chain reaction (RT-PCR) was performed in triplicate with primer sets and probes that were specific for 4 selected genes that were found to be significantly differentially expressed. These 4 genes were 2 upregulated genes: GCET2, CFTR and 2 down-regulated genes: FOXP1, GARS. cDNA was synthesized using random primers (hexamers) and Oligo 18dT and Superscript II Reverse Transcriptase (Invitrogen  [8,9], and we randomly selected 40 samples in sensitive and 30 samples in resistant group for FOXP1 staining. There was no statistical difference between these two groups of ovarian samples in age, pathological subtype, lymph node metastasis or residual tumor size (data not shown). FOXP1 staining was performed using JC12 mouse anti-human monoclonal antibodies (diluted 1:40, JC12 was kindly provided by Alison H. Banham, University of Oxford, UK [10]) using the Envision detection kit (Maixin. Bio China). Positive myoepithelial cell staining and negative stromal cell staining in breast carcinoma were used as internal positive and negative controls, respectively. FOXP1 nuclear expression was scored using the following system: negative = 0; weak/focal = 1; strong focal/widespread moderate staining = 2; or strong/widespread staining = 3. Tumors that scored 2 or 3 were considered positive for FOXP1 nuclear staining. Survival analysis was performed on those patients, and the overall survival (OS) time was defined from the date of surgery (earliest was in July, 2004) to the date of death or the last follow-up (Jun, 2014).

Enrichment analysis of DEGs
Gene Set Analysis Toolkit (Gestalt) tool was used to do enrichment analysis on the DEGs in three sections of function, biological process and molecular composition. Gestalt is a suite rich of analysis of biologically relevant content collecting eight species, including human, rat, mouse and other data from various different public data resources, such as NCBI, Ensemble, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG).

Construction of gene regulatory network
The gene regulatory network was visualized by Cytoscape [11]. Proteins in the network served as the ''nodes,'' and each pairwise protein interaction (referred to as edge) was represented by an undirected link. The property of the network was analyzed with the plug-in network analysis.

Statistics
Statistical analysis was performed using Graphpad Prism 6.0e Software for Mac OS X (GraphPad Software, La Jolla California USA, www.graphpad.com). Student's t test was employed for comparison between two groups and one-way ANOVA with Tukey's post hoc test was used for comparison between more than two groups. As to the analysis of quantitative RT-PCR result, data were expressed as mean ± SEM to compare on mRNA expression between different groups. The Chi square and Kaplan-Meier survival analysis were applied to analyze the nuclear expression of FOXP1. For these tests, a P value of\0.05 was considered statistically significant.

Gene expression analysis and clustering
The expression profiles of all the samples passed the microarray quality control (Table 3); a scatter plot was constructed with a two-dimensional rectangular coordinate plane ( Fig. 1).
Using hierarchical clustering map analysis with probe sets, the DEGs were identified in visualization, which readily distinguished the 2 groups as shown in Fig. 2. The volcano plot of DEGs revealed a total of 49 probe sets that showed a C2.0-fold change and P value B0.05, as shown in Fig. 3. Of 49 genes, 14 genes were found to be up-regulated and 35 genes down-regulated ( Table 4).
Validation of gene expression results by using quantitative RT-PCR Four highly differentially expressed genes (i.e., GCET2, CFTR, FOXP1 and SNX29) were selected for quantitative RT-PCR analysis as shown in Fig. 4. These results were in good agreement with the microarray data, confirming the reliability of the microarray results.

Validation of protein expression by immunohistochemical staining
To confirm gene expression results at the protein level, immunohistochemistry for FOXP1 was carried out on all paraffin embedded samples. For metastasis, as shown in Table 5, FOXP1 nuclear positive staining in EOC was detected in 12 of 29 EOC samples (41.4 %), while only 4 of 25 OM samples (16.0 %) showed positive nuclear staining for FOXP1 (P = 0.042). For chemotherapeutic resistance, as shown in Table 6 and depicted in Fig. 5,  GO function analysis and Signal pathway result of differential genes Significant bioprocesses of the DEGs, gene expression, biopolymer biosynthetic process, macromolecule biosynthetic process, cAMP-mediated signaling, nucleic acid metabolic process, transcription and so on (Table 7). A total of 176 KEGG pathways were enriched for the 49 DEGs, including 20 significantly enriched pathways (Table 8), such as P450 hydroxylations, HIF-1-alpha transcription factor network, mechanism of acetaminophen activity and toxicity, cystic fibrosis transmembrane conductance regulator (CFTR) and beta 2 adrenergic receptor (B2AR) pathway, alpha6beta4integrin, negative regulation of the PI3K/AKT network.

Establishment of regulatory network for the DEGs
In order to further investigate the global expression occurring and to define how individual up-or downregulated genes interact with each other to have a coordinated role, we identified potential networks for these DEGs (Fig. 7). Among the 49 DEGs, 20 were involved in the establishment of regulation network, of which 4 were up-regulated and 16 were down-regulated. A total of 21 transcription factors (TFs) were predicted, in which UBC and EP200 are the most connected predicted hub genes.

Discussion
Ovarian cancer has the highest mortality rate of all gynecologic cancers, and 75 % of patients diagnosed with ovarian cancer are already at an advanced stage.  Chemotherapy is important in treating and preventing the recurrence of ovarian cancer; however, resistance is an obstacle to overcome and finding new treatment strategies has become increasingly valuable. High-throughput technologies for assaying gene expression, such as high-density oligonucleotide and cDNA microarrays, may offer the potential to identify clinically relevant genes highly differentially expressed between different cell lines. Thus, this study showed the first communication of an investigation that involved the genome-wide examination of differences in gene expression between ovarian cancer cell lines and their sub-lines with enhanced metastasis and chemotherapic resistance. We identified 49 genes that were expressed differentially between Group A, B, C and Group 1, 2, 3, and the average change in expression level between the two groups was at least twofold. The known functions of some of these genes can provide insights with the highly metastasis and chemotherapic resistance for ovarian cancer, although others are still useful for a further research. GCET2 is found to be the most up-regulated genes in the more enhanced metastasis and chemoresistance cell group, and it is also known as human germinal center associated lymphoma (HGAL) gene, is specifically expressed in germinal center B-lymphocytes and germinal center-derived B cell lymphomas [12], but its function is largely unknown. The GCET2 gene is located on chromosome 3q13 and encodes a 178-amino acid (aa) protein with 51 % identity and 62 % similarity to the murine M17 protein [13]. GCET2 is a cytoplasmatic protein that may also associate with cell membrane. GCET2 expression is associated with improved survival in diffuse large B cell lymphoma (DLBCL) and classic Hodgkin lymphoma patients [14]. In vitro studies in human lymphocytes demonstrated that HGAL increased the binding of myosin to F-actin and inhibits the ability of myosin to translocate actin by reducing the maximal velocity of myosin head/actin movement [15]. In vitro HGAL enhances BCR signaling by binding and increasing Syk activation, in vivo older HGAL transgenic animals progressively developed polyclonal lymphoid hyperplasia and reactive AA amyloidosis [16], these finding suggests that GCET2 may play a role in humoral immune responses. No articles about expression and function of GCET2 on ovarian tissue have been published until now, and our findings make a possible insight of this gene in the study of ovarian cancer, especially about the aspects of metastasis and chemoresistance.
Cystic fibrosis transmembrane conductance regulator (CFTR, ABC35 or ABCC7) was found among the most upregulated genes in more metastasis and chemoresistance cell line group, and it participates in the beta-adrenergicdependent CFTR expression pathway. Loss of function mutations of this gene causes the autosomal recessive lethal disease cystic fibrosis (CF) and congenital bilateral aplasia of the vas deferens. There is an increasing interest in the association of cancer incidence with the genetic variations in the CFTR gene. Large cohort studies in North American and European patients with CF found that there was a marked increase in the risk of malignancies affecting the gastrointestinal tract, even to 17 times higher risk of digestive cancer with most cases arising in the bowel [17]. Meanwhile, mutations and low-penetrance polymorphisms in the CFTR gene have been found in patients with various cancers, including pancreatic cancer [18], breast cancer [19], cervical cancer [20], melanoma [21], prostate cancer [22] and lung cancer [23,24]. On the other hand, CFTR has been suggested to interact with various cancer-related kinases [25]. It encodes a member of the ATP-binding cassette (ABC) transporter superfamily. ABC proteins transport various molecules across extracellular and intracellular membranes. ABC genes are divided into seven distinct subfamilies (ABC1, MDR/TAP, MRP, ALD, OABP, GCN20, White). Meanwhile, CFTR is a member of the MRP subfamily that is involved in multidrug resistance. The encoded protein is a cAMP-activated Clchannel lining the luminal/apical surfaces of epithelial cells in airway, gut, and exocrine glands, and there is a functional coupling between CFTR and MRP2 that may be mediated by PDZ protein [24]. Taken together, our gene expression profile that show a significant up-regulated result of CFTR in more metastasis and chemoresistant ovarian cancer cell lines are consistent with previous findings, a further research on the mechanism of CFTR on ovarian cancer or   selective inhibition of CFTR or its pathway may give a insight in therapeutic effects against metastatic and chemoresistant of ovarian cancer. RBMX gene, also known as HnRNP G, is a member of heterogeneous nuclear ribonucleoprotein (hnRNP) family and can collaborate with hTra2-beta1 (human transformer-2-beta1) as sequence-specific transacting factors to exert antagonistic effects on alternative splicing which is recognized as a pivotal mechanism in regulation of gene expression and associated to tumorigenesis and metastasis of a wide variety of human cancers [26]. It is proposed that the ratio of hnRNP G/hTra2-beta1 influenced cellular splicing preference [27,28]. Some researches revealed that hnRNP G-protein showed as tumor suppressor in endometrial carcinoma [27] and oral squamous cancer [29], its activity was elicited by transactivating tumor suppressor Txnip gene [30]. A recent research showed that high frequency of hnRNP G-protein reduction and loss of expression in precancerous and human oral squamous cell carcinoma tissue specimens, suggesting that reduction in hnRNP G may play an important role in the early pathogenesis of oral squamous cell carcinomas [31].
GPRC6A encodes an orphan G-protein-coupled receptor, mediates the non-genomic effects of testosterone and other anabolic steroids in multiple tissue, and it is a potential target for developing antagonists and agonists that would have broad applications in health and disease [32], including cancer. A genome-wide association study on prostate cancer identified GPRC6A was one of the five novel genetic loci associated with prostate cancer in Japanese and Chinese Han population [33,34], and the same result was verified by a genome-wide testing of putative functional exonic variants in a multiethnic population [35]. GPRC6A is expressed at higher levels in human prostate cancer cells and prostate cancer tissues and small interfering RNA knockdown of GPRC6A attenuates these response in human prostate cancer cell lines [36]. GPRC6A is also coupled to signaling pathways, such as phosphatidylinositol 3-kinase that are known to be deregulated in prostate cancer [32]. On all accounts, nearly all researches about GPRC6A associated with cancer focused on prostate cancer, no investigation about this gene on ovarian cancer has been published.
FOXP1, as one of the down-regulated genes, drew our attention for a further research. FOXP1 is a FOX family member consisting of the winged-helix DNA-binding domain and the N-terminal transcriptional repression domain, and it is widely expressed and plays a key role in the development of various human tissues [37,38]. FOXP1 represses its target genes by forming homodimers or Kaplan-Meier survival analysis shows that the positive nuclear staining of FOXP1 is an independent risk factor in ovarian cancer patients and strongly correlates with good prognosis heterodimers with FOXP2 and FOXP4 [39], it has been suggested to be both a tumor suppressor candidate and potential oncogene, because of its differential expression levels in distinctive types of tumors, including B cell lymphomas [40], breast cancer [41,42], endometrial cancer [43], prostate cancer [44], non-small cell lung cancer [38] and renal cell carcinoma [45], the loss of FOXP1 in breast cancer has been associated with shorter survival times [42]. Until now, no article about FOXP1 expression in ovarian cancer has been published, and we made the first investigation of FOXP1 protein expression in ovarian tissue and found that nuclear staining of FOXP1 decreased as the metastasis increased, a significant decrease in FOXP1 expression in the resistance group, nuclear FOXP1 expression were independent risk factors strongly correlated with prognosis of ovarian cancer, above all, FOXP1 may serve as a good marker for late stage ovarian cancer and chemoresistance EOC patients, high expression of FOXP1 in nucleus is associated with improved survival in patients with ovarian cancer.
There are some pseudogenes which shows significant expression difference in enhanced metastasis and chemoresistant ovarian cancer cell lines, in which BC031676 and BC113708 are up-regulated, and RPL28P1, RPL23A, RPL13AP3, LOC341412, LOC641784, LOC391560, RPS16P9, LOC732186, RPL13AP23, RPLP1P7, RPL31P10, LOC648361 are down-regulated, and most of them are ribosomal protein pseudogenes. Pseudogenes are DNA sequences similar to genes encoding functional proteins but are presumed to be nonfunctional due to mutations and truncation by premature stop codons [46]. Ribosomal protein (RP) pseudogenes constitute the largest family of pseudogenes (approximately 2000 RP processed pseudogenes), and they are constitutively expressed at reasonably stable levels and are very highly conserved [47]. Although pseudogenes have long been considered as nonfunctional genomic sequences, during recent two decades, especially with the broad applications of next-generation sequencing technologies, emerging evidences have confirmed that some pseudogenes have acquired diverse functions in regulating development and diseases, especially in cancers [48]. Some pseudogenes are specifically expressed in certain cancers or diseases. It has been shown that the pseudogene of PTEN, PTENP1, was selectively lost in some human cancer cells, resulting in decreased expression of PTEN and abnormal proliferation of cancer cells [49]. The expression of MY-LKP1, a duplicated pseudogene of MYLK, can decrease the stability of MYLK mRNA at the posttranscriptional level and stimulate cell proliferation [50]. Recently, a study provided a systematic approach to analyze expressed pseudogenes, enabling comparisons of cancer versus benign tissues in multiple solid tumors, which overcome the limitations of previous analyses of pseudogene expression. They observed 218 pseudogenes expressed only in cancer samples, of which   7 Interaction network of the differentially expressed gene. Genes with more links are shown in bigger size. Proteins shown in red are encoded by upregulated genes, while those in green are encoded by downregulated genes, the gray represents the predicted genes. Arrow line represents definite control relationship, dotted line represents predicted control relationship, solid line represents inhibition 178 were observed in multiple cancers, and 40 were found to have highly specific expression in a single cancer type only, finally they described ATP8A2-J and CXADR-J pseudogenes preferentially associated with distinct subsets of breast cancer and prostate cancer patients, respectively [51]. Besides cancers, pseudogenes also involve in the development of other diseases, such as HMGA1-p [52]. Although the regulatory functions of pseudogenes seem to be striking, the functional studies of pseudogenes are still in its early stage.
The pseudogenes in our study should not be useless, their functions and relationships with ovarian cancer, especially the enhanced metastasis and chemoresistance should be investigated in near further.
In conclusion, this study has identified potential DEGs responsible for enhanced metastasis and chemoresistance in ovarian cancer cell lines. Among the 49 DEGs, 14 genes were up-regulated and 35 genes were down-regulated. Prospective investigations using a combination of genomic and proteomic approaches are required to validate the functionality of these targets identified.