Integrated Bioinformatical Analysis for Identification the invasion-related Genes in ovarian cancer

doi:10.21203/rs.3.rs-4925442/v1

Download PDF

Article

Integrated Bioinformatical Analysis for Identification the invasion-related Genes in ovarian cancer

https://doi.org/10.21203/rs.3.rs-4925442/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Ovarian cancer (OC) is the most prevalent cancer among females and its prognosis is closely related to the tumor microenvironment and the progression of the disease. This study aimed to explore the relationship between invasion-related genes (IRGs) and the prognosis of OC. Using the Cancer Genome Atlas (TCGA) database, 65 IRGs were identified and a risk model was constructed based on 4 of these IRGs. This model was able to differentiate OC patients into two risk groups with significant differences in prognosis, immune cell infiltration, and phenotypes. The model was validated using multiple methods and datasets, and the expressions of the four hub genes were confirmed using qPCR and immunohistochemistry staining assays. The results of this study suggest that the risk model could play a significant role in predicting the prognosis of OC. The findings also contribute to a better understanding of OC invasion and could assist in the development of more personalized and precise therapeutic strategies.

Biological sciences/Cancer

Biological sciences/Cancer/Cancer screening

Biological sciences/Cancer/Gynaecological cancer

ovarian cancer

invasion-related genes

bioinformatical analysis

risk model.

Ovarian cancer (OC) is the fifth leading cause of cancer mortality in women and the first cause of gynecologic tumor mortality in North America¹. Although ovariancancer mortality accelerated declined from 2017 through 2020, itwillwitness an increase of 19,710 cases in the US by 2023, accountingfor approximately 2.176% of all cancers and 13,270 deathswith a five-year survival estimate merely 47%^1,2.The patients with OC lack typical clinical symptoms in the early stages, and 75% of OC patients are diagnosed advanced disease with multiple organs invaded.Therefore, effective prediction of the invasion and metastasis of OC patients is of great significance for the clinical treatment of OC. As to FIGO, in stage I cancer, tumor is limited to ovary;in stage II, tumor extends to or is implanted on uterus or fallopian tubes; stage Ⅲ cancer shows positive retroperitoneal lymph nodes;stage Ⅳ cancer indicates pleural effusion or metastasis to extraabdominal organs such as hepatic or splenic parenchymal metastasis³. Therefore, we chose stage II~Ⅳ OC as invasion group and stage I OC as control group.

Tumor invasion is a complex multistep cell-biological process termed the invasionmetastasis cascade, which involves dissemination of cancer cells to anatomically distant organsites and their subsequent adaptation to foreign tissue microenvironments. Metastases represent the end products of a process. Each of these eventsis driven by the acquisition of genetic and/or epigenetic alterations within tumor cells and theco-option of nonneoplastic stromal cells, which together endow incipient metastatic cells with traitsneeded to generate macroscopic metastases⁴.There are many factors affecting tumor invasion in patients with OC, such as cell adhesion molecules (integrins), neurotrophic factors, cytotoxicity and cytokines (interleukin 6, insulin-like growth factor, angiogenesis, andvascular endothelial growth factor)⁵. Some monoclonal antibodies that target surface antigens often delivering toxic payloads or cells and finally there are antibodies and biologics that seek to overcome the immunosuppression caused by elements within the tumor microenvironment (TME)⁶.Therefore, advances provide provocative insights intothese cell-biological and molecular changes, which have implications regarding the steps of theinvasion-metastasis cascade that appear amenable to therapeutic targeting.

In recent years, numerous multi-center genomics discoveries have identified a group of cancer relatedgenes, which provide insight into the molecular mechanism of diseases progression and impact to fields⁷. Global views of analyzing the magnitude of information are necessary and traditional approaches to lab work have steadily been changing towards a bioinformatics era⁸. In the present study, we performed crossover analysis with gene expression datasets from The Cancer Genome Atlas (TCGA) and GeneExpression Omnibus (GEO) databases, identifiedinvasion-related genes (IRGs) in OC patients, constructed a novel IRG-based risk modelevaluating the prognostic value of IRGs in OC,and explored the relationships between IRGs and the tumor microenvironment. GEO gene expression dataset was further used toverify the prognostic risk model. Finally, four hub geneswere verified using qPCR and IHC, and their potential utility as therapeutic targets and diagnostic markers was discussed.

2.1 Data collection and processing

Public gene expression data and clinical annotations were respectively retrievedfrom TCGA (https://cancergenome.nih.gov/) and GEO (https:// www.ncbi.nlm.nih.gov/geo/) databases. The training cohort included 380 patients with OC from GEO (GSE140082) including 20 at stage I and 360 at stage II, III, IV; while 1 eligible OC cohorts fromGEO dataset (GSE26193) represented the test cohort of our study,which consisted of 20 Stage I patients and 87 Stage II, III, IV patients. Cancer geneexpression data were extracted from TCGA for further validation, including mRNA expression data fromovarian samples, Somatic mutations, CNV of 378OC patients.IRGs were from the genecards database. 9163 genes were obtained by using "tumor invasion" as the keyword and "Score" greater than 1 as the screening threshold.

Besides, 11stage I tumor specimens and 28stage II, III, IV specimens from Sichuan Provincial Hospital were used to perform differentialexpression analysis.

2.2 Consensus clustering analysis andaberrantly expressed genes identification

A consensus clustering algorithmwas used to classify OC cohorts into subtype identificationin the training set, and test the corresponding stability based on survival-relatedgenes. ConsensuClusterPlus package was used to performthe above steps and 1000 repetitions were conducted toguarantee the corresponding stability. K-M curves were plotted using the R package "survival" for patients with different subtypes. The limma package was used to analyze the differential expression of the subtypes with the best and worst prognosis against the count matrix after voom conversion.

Differentially expressed genes (DEGs) among stage I samples and stage II-III-IV samples were identified using “limma” package in R with an adjusted Pvalue< 0.001 and a |logFC|>1⁹.Weselected the top 351 upregulated and 300 downregulated genes for our study. Theheatmap and volcanoplot were performed with ‘gplots’ and ‘ggplot2’ packages of Rsoftware.

2.3 Enriching pathway of DEGs

In order to identify the pathways with the most significantly involved genes, upregulated and downregulated DEGs were conducted into DAVID¹⁰ for Gene ontology (GO)¹¹ analysis, in which the biological process (BP), molecular function (MF) enrichment and Kyoto Encyclopedia of Genes and Genomes(KEGG)¹² pathway analysis are all enriched. R software with ‘ggplot2’was used to visualize in bubble chart. ClusterProfiler package was used to performe the prognosis-related DE-IRGs. P < 0.05 was considered a significant difference.

2.4 Constructing protein-protein interaction(PPI) network and screening the hub genes

STRING(version 11.0) database online was used to identified to evaluate the interactive relationships among the DEGs with a map¹³. The combined score > 0.4 were considered as statistically significant. The intersection of prognosis-related genes, DEGs and IRGs was taken to obtain the prognosis-related genes. Cytoscape software (version 3.8.2)¹⁴ was used to visualize the high-expressed or low-expressed genes. According to the CytoHubba (version 0.1)plugging, the top 4 high-expressed genes were selected. In addition, the molecular complex detection (MCODE, version 2.0.0) plugin in Cytoscape software was used to identify of the important molecules of hub genes in PPI networks¹⁵.

2.5 Mutation of DE-IRGs in OC

The somatic mutations in OC was taken from the TCGA database and divided into low- and high-risk groups. The difference was detected as mutation annotation format (MAF), and visualized using “maftools” package of R software.

2.6 Construction of prognostic modeling and relationship between the hub genes and prognosis of OC

The survival package was used to screen DE-IRGs by univariate COX regression analysis (R package survival) and least absolute shrinkageand selection operator (LASSO, R package glmnet).A prognostic model was constructed usingleast absolute shrinkage and selection operator (LASSO)regression. Thereafter, the risk score of each patient was calculated.Using the median risk score as the cutoff, the training cohort was divided into low- and high-risk groups.Survival curves were generated by the Kaplan–Meiermethod and the two groups were compared using thelog-rank test. A time-dependent receiver operating characteristic (ROC) curve analysis was used to study themodel prediction accuracy. Cox regression was used toassess the independent prognostic value of the risk scoreand other clinicopathological features. To provide a reference for predicting the prognosis of OC patients, weused the "rms" R package to construct a nomogram basedon the risk score and clinicopathological features, and acalibration plot was used to assess the prognostic abilityof the nomogram.

2.7 Relationship between IRGs and TME

To quantifythe immune cell infiltration in each sample, single-sample Gene Set Enrichment Analysis (ssGSEA) was used toassess the enrichment of immune cells in the tumorsamples. The correlations between IRGs and immune cell content was calculated by spearman analysis.

2.8 Relationship between risk scores and clinicopathological features

Risk scores and clinical characteristics were included in univariate and multivariate Cox regression analyses to screen for independent prognostic factors, and columnar plots were constructed using the R-package rms.GO and KEGG enrichment analysis of differentially expressed genes between high and low risk groups. KEGG enrichment analysis.

2.9 Expression of prognostic genes

Demonstrate the expression levels of prognostic genes in various samples in the training and validation sets (OC early stage (stage I) and OC metastatic samples (stage II, stage III, stage IV)). The expression levels of prognostic genes in different ovarian cancer cell lines were obtained by entering prognostic genes into the CCLE database.

2.10 Protein expression level of the hub genes in OC

In this study, the protein expression ofhub genes was used IHC staining analysis. Sections from paraffin-embedded tissues of ovarian cancerwere collected, dewaxed and rehydrated. After blocked with 5% gout serum,the slides were incubatedwith primary antibodyat 4°Covernight. Then, the slides were incubated with secondaryantibody, and the images were captured using an OLYMPUS IX71 microscope. Anti-CSF1(colony-stimulating factor1) (1:200) andanti-NOS1 (Nitric oxide synthase isoform 1 ) (1:200) were purchased from Beyotime (No. AG1705; No. AF6792, China); Anti-CYP4B1( Cytochrome P450 4B1) (1:200) and anti-LPAR3(lysophosphatidic acid receptor 3) (1:200) were purchased from Proteintech (No. 11771-1-AP; No. 19509-1-AP, China);. The magnification of the IHCimages was 20×.

Statistical analysis

Data were processed, analyzed and presented using the R software (v.4.0.3) and relatedsoftware packages. Kaplan-Meier curves were used to describe the relationship between patient survival timeand survival probability. We visualized the risk-related information through charts, including risk score distribution, risk-related survival status, heat maps of prognostic genes, etc. The ROC analysis was used to analyzethe sensitivity and specificity of survival prediction using the gene signature risk score, and using AUC as anindicator of prognostic accuracy. Univariate and multivariate Cox regression analysis were used to verify theindependence of signatures. In addition, we assessed the correlation between the signature and clinical parameters. p < 0.05 was considered statistically significant.

Identification of two subtypes basedon TCGA data in OC

After preprocessing,the TCGA-OC dataset had 378 samples, theGSE140082 dataset had 380 samples (including 20 stage I OC tissues and 360 stage II-III-IV OC tissues), and theGSE26193 dataset had 107 samples (including 20 stage I OC tissues and 87 stagesII-III-IV OC tissues) (Table 1). We firstextracted the expression of 1038 invasion-related genesfrom the TCGA database and then used ConsensusClusterPlus to cluster these genes. At k=2, the samplescould be clustered together (Fig. 2A). The expression of 1038 invasion-related genes in the two subclasses is shownin Fig. 2B; 351 genes were highly expressed in the C1group and 687 genes expressed at low levels in the C2 group(Fig. 2C). Wefurther analyzed the prognostic relationship betweenthe two groups, and the results revealed significant differences between C1 and C2 (Fig. 2D; p < 0.05). The gene expression and prognosis of OC patients were analyzed by using KM for survival. It was suggested that the survival probability of C1 was higher than that of C2(P=0.015, Fig. 2D).DEGsbetween stages I and stages II-III-IV ovarian cancer were determined by the limma package in GSE140082 dataset. In total, 696 DEGs werefiltered according to the thresholds of false discoveryrate (FDR) < 0.05 and |log₂ fold change (FC)|> 1. Therewere 341 upregulated genes and 355 downregulated genes(Fig. 2E and 2F.). Then we listed the top351 upregulated and 300 downregulated genes according to the value of |log2FC|(Table S1), which demonstrated that these genes might play vital roles in theoccurrence of OC. The ‘gplots’ and ‘ggplot2’ packages of R software were used todraw heatmap and volcanoplot of thegenes (Fig. 2B,C,E,F).

PPI network construction and module analysis of DEGs

Interactions between the identified DEGs were revealed by constructing a PPInetwork. After combiningthe results of DEGs1, DEGs2 and IRGs, 65 hub genes were determined(Fig. 2A),including MMP16,NR2F2, IFIH1, FOSB, MRAS, IDO1, HTR3A, ALDH3B2, CSF1, etc (Table S2). After combiningthe results of MCODE and CytoHubba plugins, 8 hub genes were determined,including CALB2, PEG3, MCC, FBN1, MIR200C, BTC,IRF1, NOD2 (Fig. 2B). The 8 hub genes were loaded into the STRING database, to obtain the PPI data among them,and PPIs with highest interaction score(confidence > 0.3, P < 0.05) were selected. In total,there were 8 nodes and 18 edges in the network (Fig. 2B). The red circle represented positive correlation between the two genes, and the blue trianglerepresented negative correlation between the two genes.

GO function and KEGG pathway analyses for the DEGs

The results of the GO analysis indicated that in biologicalprocess terms, the upregulated DEGs were mainly enriched in regionalization, pattern specification process, forebrain development, downregulated DEGs were significantly enriched in gamma-catenin binding, co-receptor binding (Fig. 2C, Table S3). In molecular function terms,upregulated DEGs were mainly enriched in BMP receptor binding, transmembrane receptor protein serine/threonine kinase binding, and receptor serine/threonine kinase binding, whereasdownregulated DEGs were mainly enriched in gamma-catenin binding,and co-receptor binding (Table S3).

KEGG pathway analysis demonstrated that upregulated DEGs were enriched in Pathways in cancer, TGF-beta signaling pathway, whereas downregulated DEGs were significantlyenriched in Basal cell carcinoma, Hedgehog signaling pathway, Phenylalanine metabolism (Fig. 2D, Table2).

IPA classic analysis demonstrated thatupregulated DEGs were enriched in Pathways in Role of osteoblasts in Rheumatoid Arthritis Signaling pathway, Ga12/13 Signaling, whereas downregulated DEGs were significantlyenriched in Pulmaonary fibrosis idiopathic signaling pathway (Fig. 2E).

Mutations of DE-IRGs in OCs

By analyzing somatic mutation data fromTCGA-OC patients, the differences in genomic alterations were investigated. The oncoprint map showed the top 30 genes with the highest prevalence in DE-IRGs of OCs (Fig. 3A). Missense mutations were the most common mutation type. Among the DE-IRGs, LAMA1 has the highest mutation rate, while the mutation rates of other genes are relatively low. The chromosomal locations of somatic mutations in DE-IRGs are shown in Fig. 3B. According to Fig. 3C, EYA2 and MMP16 showed extensive CNV amplification, while CDH2 and CDH5 showed CNV deletion.

Construction and evaluation of a prognostic risk model

The training set data were used for further analyses.Univariate Cox analysis was performed on 287 DEGsamong molecular subtypes, and 4 genes associated with prognosis were identified as CSF1, LPAR3, CYP4B1 and NOS1. LASSO regression was used to reducethe gene numbers in the risk model. First, we analyzedthe trajectory of each independent variable, as shown inFig. 4A. Next, we used ten-fold cross validation to construct the model and confidence interval under eachlambda, as shown in Fig. 4B, C. 4 genes were chosen from the graph when lambda was 0.00204.The final 4-gene signature formula is as follows:

Riskscore=-0.2043863*EXP(CSF1)-0.0988571*EXP(LPAR3)-0.0705539*EXP(CYP4B1)-0.1775706*EXP(NOS1).

Patients were categorized into high and low risk groups based on the optimal threshold (-5.072868). PCA principal component analysis was performed on the samples, and Fig. 4D showed the ability of the risk score to discriminate between the samples. The prognostic Kaplan–Meier (KM) curveswere performed in the high- and low-risk groups. As shown in Fig. 4E, the survival rate of patients in the low-risk group was significantly higher than that of the high-risk group. The results suggested that the more patients in the training set survived and died, the higher their risk scores were,as shown in Fig. 4E and 4F.

The ROC curve (Receiver Operating Characteristic Curve) was also used toevaluate the prognostic model, and the area under the curve(AUC) was used to predict the accuracy.The AUC of the prognostic signature is 0.778, 0.709, and 0.692 for the 1-year, 2-year, and 3-year ROC curves, respectively (Fig. 4H), suggested that the prognostic model could predict the invasion of OC.

Validation of the GEO dataset

To validate the prognostic model with the screened genes as prognostic features, we used the dataset GSE26193 to evaluate the results of the training set. The risk score was calculated using the coefficients of the above genes, and the patients were categorized into high-risk and risk groups based on the optimal threshold (-7.079729).

For GSE26193, the KM curve showed that the survival rate of patients in the low-risk group was significantly higher than that of patients in the high-risk group (P < 0.05, Fig. 5A).As shown in Fig. 5B,C,D, the more patients in the validation set survived and died, the higher their risk scores were. The AUCof prognostic signature at 1, 2, and 3 years was 0.812, 0.737, and 0.628, respectively (Fig. 5E), which showedsome predictivity.

Relationship between risk score and tumor microenvironment

The tumor microenvironment significantlyaffected the therapeutic effect and prognosis of tumor.The ssGSEA was used to assess the immunestatus of the two groups. The difference in immune cells between high and low risk groups was compared through Wilcoxon analysisand seven significant immune cells were obtained, such as Type 2 T helper cell, Plasmacytoid dendritic cell, Monocyte, Macrophage, Immature dendritic cell, Gamma delta T cell, Central memory CD4 T cell(Fig. 6A).

Additionally, we further explored the correlationsbetween the risk score and immune cell contentby spearman analysis. Theresults showed that CSF1 was positively correlated with Plasmacytoid dendritic cell and Monocyte, and negativelycorrelated with Type 2 T helper cell; CYP4B1 was positively correlated with Plasmacytoiddendritic cell, and negativelycorrelated with Monocyte and Macrophage; NOS1 and LPAR3 wasnegativelycorrelated with Monocyte and Macrophage (Fig. 6B).

Univariate and multivariate analyses of the risk score for different clinical features

To assess the independence of the risk score signature modelfor clinical application, we used Univariate and multivariate Cox regression to analyze the clinical informationand the risk score. By comparing the distribution of the risk score among clinical feature groups,significant differences were found among age, stage,and grade; (Fig. 7A, F, G, H; P< 0.05) and no significant difference was found among the grade, treatment and whether histology serous (Fig. 7A; P> 0.05). The results also showed that the high-risk groups had poorer prognosis than low-risk groups. In the TCGA dataset, multivariate Coxregression analysis showed that the risk score (HR=2.641,95% CI 1.733–4.024, p <0.001) was significantly associated with prognosis ( Fig. 7B). Heatmap showed the relationship between clinical feature and risk groupings ( Fig. 7I). The findings show that therisk score signature model had good prediction in clinical.

Construction and evaluation of the prognostic risk genes

First, the Limma package was used to compare the differences in gene expression levels between the high- and low-risk groups for the training set, with the conditions of differential gene screening: |log₂FC|>0.5, P.adjust<0.05, and the volcano plots were plotted according to the results obtained, with the samples from the low-risk group as the control. Taking the samples of the low-risk group as the control, a total of 265 differentially expressed genes were obtained from the training set by differential expression analysis, of which 81 genes were up-regulated genesand 184 were down-regulated, and the heatmap and volcano map were generated in Fig. 8A and Fig. 8B.

Functional enrichment analysis of differentially expressed genes

To understand the differentially expressed genes between the high and low risk groups, GO and KEGG enrichment analysis of differentially expressed genes between the high and low risk groups was performed. As shown in Fig. 8C, 21pathways were gathered throughfunctional enrichmentanalysis with GO, such as embryonic organ morphogenesis, sensory organ morphogenesis and positive regulation of G protein-coupled receptor signaling pathway. For KEGG, six pathways were involved, such as Estrogen signaling pathway and Malignant pleural mesothelioma, et al,as shown in Fig. 8D.

Expression of the signature genes in OC tissues

To verify the accuracy of the 4-gene signature, we examined the expression of the signature genes(CSF1, CYP4B1, LPAR3 and NOS1)in clinical samples from OC patients by qPCR and IHCanalyses. The expression levels of signature genes in the training and validation sets were shown as Fig. 9Ain early OC (stage I) and Fig. 9Bin transfer OC (stage II~IV). In trainingset, theexpressions of CSF1, CYP4B1 and LPAR3were higher in stage II~IV than in stage I and NOS1 lower than in stage I. While in validation set, theexpressions of CSF1, CYP4B1LPAR3 and NOS1 were the same trends as the trainingset, but the difference was not statistically significant. The IHC showed the same trends of the four signature genes expression in both data sets (Fig. 9C).

OC invasion is the biggest obstacle to the effective cure of OC. The potential targets against OC invasion include FOXM1¹⁶, Dysbindin¹⁷, CD147 and HE4. Currently, there is no prognostic models based on invasion related genes (IRGs) have not yet effectively predicted the prognosis of OC patients.

In this study, weidentified 65 invasion-relatedgenes, which were significant differences in prognosis and the tumor microenvironment between stageⅠand stage Ⅱ~ Ⅳovarian cancers. Additionally, we developed a prognostic risk model based on the 4 hugIRGs, whichdemonstrated a powerful independent prognostic tool with better predictiveperformance. Furthermore, joint analysis of the array-based and sequence-based data of OC maybe a novel analytical strategy. Finally, we identified 4 hub genes among the model genes, and preliminarily verified them in our own IRG cohort using qPCR and IHC.

The GO analysis results indicated that the upregulated DEGs were primarilyassociated with regionalization, pattern specification process, forebrain development, BMP receptor binding, transmembrane receptor protein serine/threonine kinase binding, and receptor serine/threonine kinase binding, while the downregulated DEGs were primarily enriched in gamma-catenin binding, co-receptor binding, gamma-catenin binding,and co-receptor binding.As we all know, transmembrane receptor protein serine/threonine kinase binding can transmit extracellular signals into the cell, and then affects gene transcription to achieved a variety of biological functions. Its main ligands are transforming growth factor-βs (TGF-βs). Under pathological conditions, overexpressed TGF-β caused epithelial-mesenchymal transition (EMT), extracellular matrix (ECM) deposition, cancer-associated fibroblast (CAF) formation, whichled to fibrotic disease, and cancer. Given the critical role of TGF-β and its downstream molecules in the progressionof fbrosis and cancers, therapeutics targeting TGF-β signaling appears to be a promising strategy¹⁸. Furthermore, the enriched KEGG pathways of upregulated DEGsincluded pathways in cancer, TGF-β signaling pathway.The result of KEGG supportedthe GO analysis. Downregulated DEGswere related to Basal cell carcinoma, Hedgehog signaling pathway, Phenylalanine metabolism.The Hedgehog signaling pathway is essential in the intricate process of embryonic development. Dysregulation of this pathway has been linked to developmental abnormalities and various forms of cancer, such as Gorlin syndrome and sporadic malignancies¹⁹. The IPA classic analysis demonstrated thatupregulated DEGs were enriched in Pathways in Role of osteoblasts in Rheumatoid Arthritis Signaling pathway, Gα12/13 Signaling, whereas downregulated DEGs were significantlyenriched in Pulmaonary fibrosis idiopathic signaling pathway.Heightened levels of both Gα12 and Gα13, along with their augmented signaling, have been linked to the initiation and advancement of various malignancies. Moreover, the regulation of Gα12/13 protein expression is governed by a multitude of transcriptional and post-transcriptional processes, many of which are commonly disrupted in cancer.With the burgeoning comprehension of tumorigenic mechanisms orchestrated by Gα12/13 proteins, it is evident that directing Gα12/13 signaling in a tailored fashion could offer a novel approach to enhance therapeutic efficacy in various solid malignancies²⁰.

In this study, laminin chain genes LAMA1 had the highest mutation rate. EYA2 and MMP16 showed extensive CNV amplification, while CDH2 and CDH5 displayed CNV deletion. LAMA1 has mononucleotide repeats in the coding sequences that might be mutation targets in the cancers with microsatellite instability (MSI)²¹. Biallelic mutations in LAMA1 was the cause of Cerebellar dysplasia with cysts (CDC)²² and Poretti-Boltshauser Syndrome(PTBHS)²³. Previous studies indicated thatEYA2 suppresses the progression of hepatocellular carcinoma²⁴, lung cancer²⁵, and breast cancer²⁶.FBXO7 bound and stabilized SIX1 co-transcriptional regulator EYA2, stimulatingmesenchymal gene expression and antigen presentationmachinery.Targeting EYA2 Tyr phosphatase activity decreased mesenchymal phenotypes and enhanced cancer cell immunogenicity, resulting in attenuated tumor growth and metastasis,while increased infiltration of cytotoxic T and NK cells, and enhanced anti-PD-1 therapy response in mouse tumormodels²⁷. Nevertheless, the clinical significance of these genetic modifications must be validated through a more extensive sample size, and the precise mechanism of these genetic alterations necessitates both in vitro and in vivo verification. MMPs were involved in angiogenesis, which promoted cancer cell growth and migration. Targeting MMP16 could inhibits the Cancer Cell Proliferation, Migration and Invasion in Cutaneous Squamous Cell Carcinoma{Mu, 2020 #326}. At the primary tumor site, tumor cells have the capacity to acquire mesenchymal characteristics, a phenomenon known as the epithelial-to-mesenchymal transition (EMT). This intricate process is partly under the control of cadherin (CDH) genes, and it plays a crucial role in the development of metastatic lesions. High expressions of CDH1/2/4/11/12/13 were significantly correlated with cell adhesion, the extracellular matrix remodeling process, the EMT, WNT/beta-catenin, and interleukin-mediated immune responses, those were thought to be potential biomarkers for breast cancer progression and metastasis. {Ku, 2022 #328}.

Most of the IRGs in the risk model have been reported to be associated with cancer. As a result, we identified four hub genes: CSF1, CYP4B1, LPAR3 and NOS1. The potential mechanisms of these four hub genes in OC deserve further discussion.CSF1 can promotethe differentiation and survival of macrophages. The c-fms proto-oncogene encodes for thetyrosine kinase receptor for CSF1. CSF1 also plays an important role in lactation, ovulation, preimplantation,placental function with trophoblastic invasion and regulation of the estrous cycle²⁸. The CSF1 op/op mouse (null mutant for CSF1) displaysimpaired mammary development and impaired reproductive phenotype, with a near-absenceof tissue macrophages²⁹.CSF1 and/or its receptor have beenfound to be expressed by the tumor epithelium in several human epithelial cancers, includingthose of breast, lung, endometrial, prostatic, trophoblastic and ovarian³⁰. Elevatedlevels of CSF1 and c-fms are associated with poor prognosis in several epithelial cancers, such as Gastric Cancer and ovarian caner^31,32. It has reported that CSF-1-stimulated invasiveness , including ovarian cancers, is mediated through theactions of the urokinase plasminogen activator (uPA)³³. To date, blockade of colony stimulating factor-1 or its receptor represents the only truly selective approach to manipulate macrophages in cancer patients³⁰, such as PLX3397 (PLX108-01)³⁴.

CYP4B1,a mammalian cytochrome P450 monooxygenase, belongs to the CYP4 family, acts at the interface of endo- and xenobiotic metabolism³⁵. The contribution of CYP4B1 in cancer is of particular interest as CYP4B1 gene expression was found tobe altered and presumably associated with certaincancers³⁶.The involvement of CYP4B1 in cancer was considered via activation of procarcinogens and neovascularization.The silencing of CYP4B1 gene expression caused a decline in neovascularization. Upon injury of the corneal epithelium, the CYP4B1enzymatic pathway was activated leading to the production of 12-HETrE, which was a potent angiogeniceicosanoid³⁷.Higher levelof CYP4B1 mRNA was found in recurrent serousovarian cancer patients than that of cured ovariancancer patients³⁸. Suicide gene approach combining reengineered CYP4B1 andprodrug 4-ipomeanol (4-IPO) has shown a promising potential for targeted cancer therapy³⁹.

LPAR3 is capable of encoding the third Gprotein-coupled receptor for LPA⁴⁰.LPA is speciallycorrelated to ovarian malignancies, as presented by the highconcentration observed in the ascites of patients with ovarian cancer⁴¹. High circLPAR3 expression ispositively correlated with esophageal squamous cell carcinoma (ESCC)stage⁴². CircLPAR3 promoted ESCC cell migration, invasion, and metastasis. LPAR3 increased the phosphorylation of Akt and MAPK, affect the RAS/MAPK and PI3K/Akt pathways through weaking its suppression on c-MET expression. miR-15b affected the proliferation, apoptosis, and senescence of human ovarian cancer cells by binding to LPAR3 with the involvement of the PI3K/Akt pathway⁴³.LPAR3 expression increased along with poorer differentiation⁴¹.

NO (Nitric oxide), is a soluble endogenous gas working as a critical metabolic regulator in metabolism reprogramming of tumorgenesis⁴⁴. NOS1 was observed to be increasingly expressed in various types of cancer, and its expression has been associated with tumor progression. The GEO database revealed that NOS1 expression level was higher in ovarian cancer compared with paired normal ovarian tissues. NOS1 affected cell function in vitro, including proliferation, migration and invasion as well as chemoresistance to cispatin (DDP) treatment in OVCAR3 cells⁴⁵.

Theresults of ssGSEAshowed that CSF1 was positively correlated with Plasmacytoid dendritic cell and Monocyte, and negativelycorrelated with Type 2 T helper cell; CYP4B1 was positively correlated with Plasmacytoiddendritic cell, and negativelycorrelated with Monocyte and Macrophage; NOS1 and LPAR3 wasnegativelycorrelated with Monocyte and Macrophage. The tumor microenvironment, as the hotbed of thetumor, has significant immune cell infiltration. In accordance with the complexity of the tumor microenvironment, immune cells recruited to the tumor tissues havedual tumor-promoting and tumor-antagonizing characteristics.TME is composed of distinct cellpopulations in a complex matrix, which is characterized byinefficiencies of oxygen and nutrition delivery due to limitedor poorly differentiated vasculature. Inthis harsh TME, infiltrating immune cells are forced to undergorelevant metabolic adaptations, which in turn disrupt the antitumor effects of immune cells{DeBerardinis, 2016 #330}{Liu, 2022 #331}. The immune microenvironment plays a keyrole in tumor development and treatment⁴⁶.MacDonald KPA⁴⁷ et al have reported c-fms is up-regulated in all dendritic cell subsets during differentiation. Colony-stimulating factor 1 (CSF1) is a primary regulator of the survival, proliferation, and differentiation of monocyte/macrophage that sustains the protumorigenic functions of tumor-associated macrophages (TAMs)⁴⁸. Schneider KM⁴⁹ et al stated that chronically elevated levels of glucocorticoids drove the generation of an inflammatory subset of enteric glia that promoted monocyte- and TNF-mediated inflammation via CSF1. CSF1-dependent DC subsets link the innateand adaptive immune responses in T helper 2 (Th2)cell-mediated allergic lung inflammation⁵⁰. In line with this, CSF1 was positively correlated with Plasmacytoid dendritic cell and Monocyte, and negativelycorrelated with Type 2 helper cell. To ate, little is know about CYP4B1 correlated with Plasmacytoiddendritic cell, Monocyte and Macrophage. It has been reported that NOS1-derivedNO facilitates low-density lipoprotein (LDL) uptake by macrophages and has been identified as a regulator of macrophage and endothelial cell interaction⁵¹, which is opposite to our findings. Pei S⁵²et al detected thatin vitro experiments withco-culturedmonocytesandneutrophilsdemonstrated that monocytes from Lpar3(-/-) mice promoted the formation of Neutrophil Extracellular Traps(NETs), suggesting that LPA(3) acting on monocytes inhibits the formation of NETs. Chen L⁵³et al confirmed that the LPA-LPA 1, 3 receptor-HNF1α pathway was involved in the downregulation of Fut8, leading to diminished foam cell migration. Modulation of the Agpat4/LPA/p38/p65 axis regulated macrophage polarization, T-cell activity and colorectal cancer (CRC) progression. Notably, combined therapy with LPA and regular chemotherapy drugs synergistically suppressed CRC development⁵⁴.These conclusions further indicated that IRGs might play importantroles in the altered immune response in OC.

Regarding therelationship between the clinicalmanifestations of patients with prognosis of OC. The results showed that age, stage,and grade were found the independent influencing factors. Prognosis of ovarian cancer is typically determined by the cancer stageand grade. Theincidence of ovarian cancer and mortalityrates increase with age. Most cases of ovariancancer occur in women older than 50 years,but the diagnosis may be made at any age,including infancy. The incidence and mortality of OC were 6.6 and 2.0 per 100,000 women respectively when age at 20 to 49 years old, while 55.6 and 55.2per 100,000 women respectively when age over 75 years old {ROETT, 2009 #329}.

TheqPCR and IHC. Our verification results confirmed that the 4 hub genes were both overexpressed in OC tissues. Thereafter, using X-tileand the 3-yearOS of patients, we determined the optimal cutoff valuefor CSF1, CYP4B1, LPAR3 and NOS1expression. Based on each cutoffvalue, we subdivided the OC patients into high- andlow-expression groups, and further investigated the differences in the tumor microenvironment betweenthe pairs of groups. We found that the prognosis of thehigh-expression groups was poorer. Finally, univariate andmultivariate Cox regression analyses indicated that the genes were independent prognostic factorsin OC. These experimental results showed that highexpression of the hub genes was closely related to theprognosis and immunosuppressive microenvironmentof OC, which was consistent with our bioinformaticsanalyses, and further verified the stability and accuracy ofthe risk model.

In summary, we performed unsupervised clustering analysis and survival analysis between stage I and II-III-IV OC to obtain the differentially expressed genes in the training set. An invasion-related 4-genesignature was established to evaluate the prognosis of OC patients. Furthermore, to verify the signature, we used the GSE26193datasets for testingand finally obtained a good risk prediction effect in thosedatasets. Moreover, the results of the qPCR and immunohistochemistry staining assays revealed that CSF1, CYP4B1, LPAR3 and NOS1 were upregulatedin OC tissues. These results suggest that the signature couldpotentially be used to evaluate the prognostic risk of OCpatients. Nevertheless, experimental verification and direct clinical application tests are needed with larger sample sizes

All procedures followed the 2013 Declaration of Helsinki. The Hospital Ethics Committee has approved the study, and all participants gave their consent for participation, with the signed informed consent provided (No: 2022-445).

Acknowledgements

We would like to thank TCGA and GEO databases for the availability of thedata.We acknowledge support from the Department of Gynecology and Obstetrics, Sichuan Provincial People’s Hospital and Sichuan Provincial Key Laboratory for Human Disease Gene Study.

Author contribution

Dan Su designed and supervised the present study; Dan Su, Xiaoxia Guo and Xiaojing Wu wrote the main manuscript text. Xiaoxia Guo, Xiaojing Wu and Ping Yu, Yao Xie and Jie Mei performed the formal analysis and methodology. Shiqing Yu, Lijuan Kong, Congrui Liu and Xinru Feng collected and analyzed the data with the Softwares. Ping Yu, Yao Xie and Jie Mei prepared figures 1-3. Dan Su, Shiqing Yu, Lijuan Kong, Congrui Liu and Xinru Feng prepared figures 4-8. Xiaoxia Guo and Xiaojing Wu prepared figures 9. All authorsreviewed and approved the final manuscript.

Funding

This work was supported by the Science & Technology Department of Sichuan Province (grant number: 2022YFS0088), Sichuan Provincial Cadre Health Research Project (grant number: 2022-212) and the Science & Technology Department of Chengdu, Sichuan (grant number: 2019-YF05-00244-SN).

Competing interests

The authors declare no competing financial interests and competing non-financialinterest.

Data Availability Statement

High-throughput gene expression data of OC wereextracted from the TCGA Data Portal (https://tcga-data.nci.nih.gov/tcga). The mRNA expression data from ovary samples of OC patients (Ovary), somatic mutations, CNV and IRGs from genecards database. In addition, gene expression profiles of datasets GSE140082and GSE26193 were obtained from GEO. The data underlying this article will be shared on reasonable request to the corresponding author.

Siegel RL, Miller KD, Wagle NS, Jemal A. Cancer statistics, 2023. CA Cancer J Clin 2023; 73(1): 17-48.
Jazwinska DE, Kulawiec DG, Zervantonakis IK. Cancer-mesothelial and cancer-macrophage interactions in the ovarian cancer microenvironment. Am J Physiol Cell Physiol 2023; 325(3): C721-C30.
Javadi S, Ganeshan DM, Qayyum A, Iyer RB, Bhosale P. Ovarian Cancer, the Revised FIGO Staging System, and the Role of Imaging. AJR Am J Roentgenol 2016; 206(6): 1351-60.
Valastyan S, Weinberg RA. Tumor metastasis: molecular insights and evolving paradigms. Cell 2011; 147(2): 275-92.
Palma Flores C, Garcia-Vazquez R, Gallardo Rincon D, et al. MicroRNAs driving invasion and metastasis in ovarian cancer: Opportunities for translational medicine (Review). Int J Oncol 2017; 50(5): 1461-76.
Bilotta MT, Antignani A, Fitzgerald DJ. Managing the TME to improve the efficacy of cancer therapy. Front Immunol 2022; 13: 954992.
Levy SE, Boone BE. Next-Generation Sequencing Strategies. Cold Spring Harb Perspect Med 2019; 9(7): a025791-a802.
Tao Z, Shi A, Li R, Wang Y, Wang X, Zhao J. microarray bioinformatics in cancer- a review. J BUON 2017; 22(4): 838-43.
Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015; 43(7): e47.
Jiao X, Sherman BT, Huang DW, et al. DAVID-WS: A Stateful Web Service to Facilitate Gene/Protein List Analysis. Bioinformatics 2012; 28(13): 1805-6.
The Gene Ontology C. The Gene Ontology Resource: 20 years and still GOing strong. Nucleic Acids Res 2019; 47(D1): D330-D8.
Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res 2017; 45(D1): D353-D61.
M K, S G. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000; 28(1): 27-30.
Doncheva NT, Morris JH, Gorodkin J, Jensen LJ. Cytoscape StringApp: Network Analysis and Visualization of Proteomics Data. J Proteome Res 2019; 18(2): 623-32.
GD B, CW H. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics 2003; 4: 2.
Liu C, Barger CJ, Karpf AR. FOXM1: A Multifunctional Oncoprotein and Emerging Therapeutic Target in Ovarian Cancer. Cancers (Basel) 2021; 13(12).
Lv X, Guo X, Ru Y, et al. Dysbindin facilitates invasion and metastasis by promoting phosphorylation of ERK in epithelial ovarian cancer. J Cancer 2020; 11(10): 2821-9.
Peng D, Fu M, Wang M, Wei Y, Wei X. Targeting TGF-beta signal transduction for fibrosis and cancer therapy. Mol Cancer 2022; 21(1): 104.
Skoda AM, Simovic D, Karin V, Kardum V, Vranic S, Serman L. The role of the Hedgehog signaling pathway in cancer: A comprehensive review. Bosn J Basic Med Sci 2018; 18(1): 8-20.
Rasheed SAK, Subramanyan LV, Lim WK, Udayappan UK, Wang M, Casey PJ. The emerging roles of Galpha12/13 proteins on the hallmarks of cancer in solid tumors. Oncogene 2022; 41(2): 147-58.
Choi MR, An CH, Yoo NJ, Lee SH. Laminin gene LAMB4 is somatically mutated and expressionally altered in gastric and colorectal cancers. APMIS 2015; 123(1): 65-71.
Aldinger KA, Mosca SJ, Tetreault M, et al. Mutations in LAMA1 cause cerebellar dysplasia and cysts with and without retinal dystrophy. Am J Hum Genet 2014; 95(2): 227-34.
Schiff ER, Aychoua N, Nutan S, et al. Variability of retinopathy consequent upon novel mutations in LAMA1. Ophthalmic Genet 2022; 43(5): 671-8.
Liu ZK, Li C, Zhang RY, et al. EYA2 suppresses the progression of hepatocellular carcinoma via SOCS3-mediated blockade of JAK/STAT signaling. Mol Cancer 2021; 20(1): 79.
Li Z, Qiu R, Qiu X, Tian T. EYA2 promotes lung cancer cell proliferation by downregulatingthe expression of PTEN. Oncotarget 2017; 8(67): 110837-48.
Ren L, Guo D, Wan X, Qu R. EYA2 upregulates miR-93 to promote tumorigenesis of breast cancer by targeting and inhibiting the STING signaling pathway. Carcinogenesis 2021.
Shen JZ, Qiu Z, Wu Q, et al. A FBXO7/EYA2-SCF(FBXW7) axis promotes AXL-mediated maintenance of mesenchymal and immune evasion phenotypes of cancer cells. Mol Cell 2022; 82(6): 1123-39 e8.
Munoz-Garcia J, Cochonneau D, Teletchea S, et al. The twin cytokines interleukin-34 and CSF-1: masterful conductors of macrophage homeostasis. Theranostics 2021; 11(4): 1568-93.
Ryan GR, Dai XM, Dominguez MG, et al. Rescue of the colony-stimulating factor 1 (CSF-1)-nullizygous mouse (Csf1(op)/Csf1(op)) phenotype with a CSF-1 transgene and identification of sites of local CSF-1 synthesis. Blood 2001; 98(1): 74-84.
Ries CH, Hoves S, Cannarile MA, Ruttinger D. CSF-1/CSF-1R targeting agents in clinical development for cancer therapy. Curr Opin Pharmacol 2015; 23: 45-51.
Liu H, Zhang H, Shen Z, et al. Increased Expression of CSF-1 Associates With Poor Prognosis of Patients With Gastric Cancer Undergoing Gastrectomy. Medicine (Baltimore) 2016; 95(9): e2675.
Woo HH, Chambers SK. Human ALKBH3-induced m(1)A demethylation increases the CSF-1 mRNA stability in breast and ovarian cancer cells. Biochim Biophys Acta Gene Regul Mech 2019; 1862(1): 35-46.
Chambers SK, Wang Y, Gertz RE, Kacinski BM. Macrophage Colony-stimulating Factor Mediates Invasion of Ovarian Cancer Cells through Urokinase. Cancer Res 1995; 55: 1578-85.
Yao W, Cao Q, Luo S, et al. Microglial ERK-NRBP1-CREB-BDNF signaling in sustained antidepressant actions of (R)-ketamine. Mol Psychiatry 2022; 27(3): 1618-29.
Roder A, Husken S, Hutter MC, et al. Spotlight on CYP4B1. Int J Mol Sci 2023; 24(3).
Eun HS, Cho SY, Lee BS, Seong IO, Kim KH. Profiling cytochrome P450 family 4 gene expression in human hepatocellular carcinoma. Mol Med Rep 2018; 18(6): 4865-76.
Seta F, Patil K, Bellner L, et al. Inhibition of VEGF expression and corneal neovascularization by siRNA targeting cytochrome P450 4B1. Prostaglandins Other Lipid Mediat 2007; 84(3-4): 116-27.
Barlin JN, Jelinic P, Olvera N, et al. Validated gene targets associated with curatively treated advanced serous ovarian carcinoma. Gynecol Oncol 2013; 128(3): 512-7.
Wiek C, Schmidt EM, Roellecke K, et al. Identification of amino acid determinants in CYP4B1 for optimal catalytic processing of 4-ipomeanol. Biochem J 2015; 465(1): 103-14.
Diao H, Li R, El Zowalaty AE, et al. Deletion of Lysophosphatidic Acid Receptor 3 (Lpar3) Disrupts Fine Local Balance of Progesterone and Estrogen Signaling in Mouse Uterus During Implantation. Biol Reprod 2015; 93(5): 123.
Yu X, Zhang Y, Chen H. LPA receptor 1 mediates LPA-induced ovarian cancer metastasis: an in vitro and in vivo study. BMC Cancer 2016; 16(1): 846.
Shi Y, Fang N, Li Y, et al. Circular RNA LPAR3 sponges microRNA-198 to facilitate esophageal cancer migration, invasion, and metastasis. Cancer Sci 2020; 111(8): 2824-36.
Li GC, Qin XL, Song HH, et al. Upregulated microRNA-15b alleviates ovarian cancer through inhitbition of the PI3K/Akt pathway by targeting LPAR3. J Cell Physiol 2019; 234(12): 22331-42.
López-Sánchez LM, Aranda E, Rodríguez-Ariza A. Nitric oxide and tumor metabolic reprogramming. Biochemical Pharmacology 2020; 176.
Zou Z, Li X, Sun Y, et al. NOS1 expression promotes proliferation and invasion and enhances chemoresistance in ovarian cancer. Oncol Lett 2020; 19(4): 2989-95.
Anderson NM, Simon MC. The tumor microenvironment. Curr Biol 2020; 30(16): R921-R5.
MacDonald KPA, Rowe V, Bofinger HM, et al. The Colony-Stimulating Factor 1 Receptor Is Expressed on Dendritic Cells during Differentiation and Regulates Their Expansion. The Journal of Immunology 2005; 175(3): 1399-405.
Fujiwara T, Yakoub MA, Chandler A, et al. CSF1/CSF1R Signaling Inhibitor Pexidartinib (PLX3397) Reprograms Tumor-Associated Macrophages and Stimulates T-cell Infiltration in the Sarcoma Microenvironment. Mol Cancer Ther 2021; 20(8): 1388-99.
Schneider KM, Blank N, Alvarez Y, et al. The enteric nervous system relays psychological stress to intestinal inflammation. Cell 2023; 186(13): 2823-38 e20.
Moon HG, Kim SJ, Jeong JJ, et al. Airway Epithelial Cell-Derived Colony Stimulating Factor-1 Promotes Allergen Sensitization. Immunity 2018; 49(2): 275-87 e5.
Roy A, Banerjee S, Saqib U, Baig MS. NOS1-derived nitric oxide facilitates macrophage uptake of low-density lipoprotein. J Cell Biochem 2019; 120(7): 11593-603.
Pei S, Xu C, Pei J, et al. Lysophosphatidic Acid Receptor 3 Suppress Neutrophil Extracellular Traps Production and Thrombosis During Sepsis. Front Immunol 2022; 13: 844781.
Chen L, Zhang J, Yang X, Liu Y, Deng X, Yu C. Lysophosphatidic acid decreased macrophage foam cell migration correlated with downregulation of fucosyltransferase 8 via HNF1alpha. Atherosclerosis 2019; 290: 19-30.
Zhang D, Shi R, Xiang W, et al. The Agpat4/LPA axis in colorectal cancer cells regulates antitumor responses via p38/p65 signaling in macrophages. Signal Transduct Target Ther 2020; 5(1): 24.

Table 1. Details of OC studies and associated microarray datasets from GEO database

Compare groups

Publication

Total differentially expressed genes

Up-regulated

Down-regulated

Technology/

Platform

Age

GSE140082

cluster1 and cluster2

Clinical Cancer Research

1038

351

687

GPL14951 Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip

21~80

Table 2. KEGG pathway analysis of the 5 genes associated with OC

ID	Term	Count	P value
hsa04350	TGF-beta signaling pathway	4	0.000942
hsa05200	Pathways in cancer	7	0.001255
hsa05217	Basal cell carcinoma	3	0.002873
hsa04340	Hedgehog signaling pathway	3	0.003025
hsa00360	Phenylalanine metabolism	2	0.003467

No competing interests reported.

Supplementarymaterial.docx

Download PDF

Editorial decision: Revision requested
30 Sep, 2024
Reviews received at journal
23 Sep, 2024
Reviewers agreed at journal
23 Sep, 2024
Reviewers invited by journal
17 Sep, 2024
Editor assigned by journal
17 Sep, 2024
Editor invited by journal
10 Sep, 2024
Submission checks completed at journal
07 Sep, 2024
First submitted to journal
16 Aug, 2024

You are reading this latest preprint version

Integrated Bioinformatical Analysis for Identification the invasion-related Genes in ovarian cancer

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials and Methods

3. Results

4. Discussion

Declarations

References

Tables

Additional Declarations

Supplementary Files

Status:

Version 1