Identification of MEG8/miR‐378d/SOBP axis as a novel regulatory network and associated with immune infiltrates in ovarian carcinoma by integrated bioinformatics analysis

Abstract Background To investigate the potential molecular mechanism of ovarian cancer (OC) evolution and immunological correlation using the integrated bioinformatics analysis. Methods Data from the Gene Expression Omnibus was used to gain differentially expressed genes (DEGs). Gene Ontology and Kyoto Encyclopedia of Gene and Genome pathway analysis were completed by utilizing the Database for Annotation, Visualization, and Integrated Discovery. After multiple validations via The Cancer Genome Atlas (TCGA), Genotype‐Tissue Expression (GTEx) projects, the Human Protein Atlas, Kaplan–Meier (KM) plotter, and immune logical relationships of the key gene SOBP were evaluated based on Tumor Immune Estimation Resource, and Gene Set Enrichment Analysis (GSEA) software. Finally, the lncRNAs‐miRNAs‐mRNAs subnetwork was predicted by starBase, TargetScan, miRBD, and LncBase, individually. Correlation of expression and prognosis for mRNAs, miRNAs, and lncRNAs were confirmed by TCGA, Gene Expression Profiling Interactive Analysis 2 (GEPIA 2), starBase, and KM. Results A total of 192 shared DEGs were discovered from the four data sets, including 125 upregulated and 67 downregulated genes. Functional enrichment analysis presented that they were mainly enriched in cartilage development, pathway in PI3 K‐Akt signaling pathway. Lower expression of SOBP was the independent prognostic factor for inferior prognosis in OC patients. The downregulation of SOBP enhanced the infiltration levels of B cells, CD8+ T cells, Macrophage, Neutrophil, and Dendritic cells. GSEA also disclosed low SOBP showed a significantly associated with the activation of various immune‐related pathways. Finally, we first reported that the MEG8/miR‐378d/SOBP axis was linked to the development and prognosis of OC through regulating the cytokines pathway. Conclusions Our study establishes a novel MEG8/miR‐378d/SOBP axis in the development and prognosis of OC, and the triple subnetwork probably affects the progression of the OC by regulating the cytokines pathway.


| INTRODUCTION
Ovarian cancer (OC) is the eighth common cancer in females and the second cause of gynecologic tumor death. Globally, 295,414 (1.6% of all cancer cases) new OC patients and 184,799 (1.9% of all cancer cases) deaths occurred in 2018. 1 Most OC patients were diagnosed with advanced-stage tumor (stage III or IV), and surgery and platinum/taxane chemotherapy are still the main treatment for OC patients. 2 However, over 75% of women with advanced-stage OC will die from disease recurrence. 3 Therefore, it is meaningful to investigate early screening biomarkers and potential therapeutic targets to improve outcomes of OC patients.
Tumor-infiltrating lymphocytes (TILs) play a crucial role in the immune responses of human cancers. 4 The certain antigens in the cell membrane divided TILs into several subsets, and each of the subsets has different functions in the development of malignant tumors. Several studies have reported that TILs improved prognosis in various types of malignant tumors, such as breast, lung, and colon cancer. [5][6][7] In addition, several studies have also shown a distinct correlation among different TIL intensities regarding the survival outcomes. [8][9][10] However, the clinical effect of immune cells in OC patients remains unclear. Therefore, the role of TILs in treatment decision making and prognostic assessment deserves further study.
The competing endogenous RNA (ceRNA) hypothesis was introduced by Salmena et al. in 2011, which described a brand-new regulatory mechanism among ncRNAs and mRNAs. 11 CeRNA hypothesis sustained that lncRNAs and mRNAs could cross talk by competitively binding to common miRNAs to achieve their biological functions. 12 But the characters of the lncRNA-miRNA-mRNA network in OC were rarely explored. In addition, the potential mechanism of action and relationship with immunity in the ceRNA regulatory network is also worth exploring.
This study aimed to identify useful signature genes that were related to the development, metastasis, and prognostic values of OC, and investigated the underlying mechanism of hub genes through a comprehensive bioinformatics approach. First, we combined four gene expression sets (GSE36668, GSE12470, GSE14407, and GSE27651) from Gene Expression Omnibus (GEO) to obtain common differentially expressed genes (DEGs). Second, Gene Ontology (GO) and Kyoto Encyclopedia of Gene and Genome (KEGG) pathway analysis were utilized to search the function of the common DEGs in OC. Third, six hub genes (RIPK4, SCNN1A, SLC4A11, ELF3, CLDN4, and SOBP) were found after confirmation of the expression and survival of multiple databases. Especially, lower expression of the SOBP gene was independently associated with inferior survival outcomes in OC patients. In addition, high levels of immune infiltration and activation of immune pathways were associated with the downregulation of SOBP. Finally, we predicted the MEG8/ miR-378d/SOBP (lncRNA-miRNA-mRNA) axis using integrated bioinformatics analysis, which was associated with the development and prognosis of OC by regulating immune responses in the cytokines pathway.

| Differential expression analysis
We found four ovarian tumor data sets (GSE36668, GSE12470, GSE14407, and GSE27651) from the GEO database (http://www.ncbi.nlm.nih.gov/geo). 13 The detailed data set information of the four OC data sets are shown in Table 1. The DEGs were obtained and visualized by conducting affy, limma, and ggplot2 packages in R software (version 3.5.3). |log2FC| >1.5 and adjusted p-value <0.05 means statistically significant. VENNY (http://bioin fogp.cnb.csic.es/tools/ venny/ index.html) 2.1.0 was applied to get common DEGs and draw the Venn diagrams.

| Functional analysis of common DEGs
For a further understanding of the common DEGs, the Database for Annotation, Visualization, and Integrated Discovery (DAVID) was utilized to perform GO and KEGG enrichment analyses, 14 and the results were visualized with ggplot2 and GO plot R packages. To identify the potential mechanism of SOBP on prognosis of OC, gene set enrichment analysis was conducted by Gene Set Enrichment Analysis (GSEA) software version 3.0. 15 We Conclusions: Our study establishes a novel MEG8/miR-378d/SOBP axis in the development and prognosis of OC, and the triple subnetwork probably affects the progression of the OC by regulating the cytokines pathway.

K E Y W O R D S
bioinformatics analysis, competing endogenous RNA, differently expressed genes, ovarian cancer, tumor-infiltrating lymphocytes divided 379 OC samples from The Cancer Genome Atlas (TCGA) into high SOBP group and low SOBP group by median value, other parameters criteria were set by default. Nominal p < 0.01, false discovery rate (FDR)<0.25, and |Normalized Enrichment Score (NES)| >1 were considered as statistically significant.

| Gene expression profiling interactive analysis 2 (GEPIA2)
GEPIA2 (http://gepia2.cance r-pku.cn) is an online tool for analyzing the transcriptional profiles of human cancers and normal tissues, by using the TCGA database and the Genotype-Tissue Expression (GTEx) projects. 16 We performed the analysis regarding the differential expression and prognostic values for common DEGs, and the relationship between MEG8 and SOBP was analyzed in the correlation analysis module. Genes with |log2FC| >1 and p < 0.05 were set as the cutoff criterion.

| Kaplan-Meier (KM) plotter
The KM Plotter (http://kmplot.com/analy sis/) was applied to verify the prognosis of the six key genes in OC. 17 The KM Plotter could assess the influence of 54,675 genes on survival, including 6234 breast cancers, 3452 lung cancers, 2190 OC, and 1440 gastric cancers. The prognostic value of the mRNAs, miRNAs, and lncRNAs were showed with a hazard ratio (HR), 95% confidence interval (CI), and p-value. The p-value <0.05 was reflected statistically significant.

| Immunohistochemistry
We used the Human Protein Atlas (https://www.prote inatl as.org/), which had mRNA and protein expression data on 44 different human normal tissues, as well as 17 common types of cancer. 18 The antibody-based protein data show the corresponding protein expression. The intensity of the staining was utilized to evaluate the six key genes in OC tissues at protein expression levels.

| Oncomine
Oncomine (http://www.oncom ine.org), 19 an online database, was used to analyze the difference in mRNA expression of SOBP among normal tissues and cancer samples. The thresholds of p-value and fold change (FC) were set as 0.001 and 1.5, respectively.

| Tumor immune estimation resource (TIMER) analysis
To date, surgery and platinum/taxane chemotherapy are still the main treatment for OC patients. 2 However, there is still no effective target index and drug application in the field of immunotherapy. Increasing evidence has indicated a close association between immune infiltration in cancer and clinical outcomes. 20,21 TIMER database (https://cistr ome.shiny apps.io/timer/) was used to systematically analyze the tumorinfiltrating immune cells in human cancers using more than 10,000 samples from the TCGA database. 22 We explored the association between SOBP expression and the abundance of infiltrating immune cells. The marker genes used for the analysis were based on data from previous studies. 23,24 TIMER also output multivariate Cox analysis results, including hazard ratios and statistical significance. p < 0.05 was considered statistically significant.

| Prediction and evaluation of upstream miRNAs and lncRNAs
Upstream miRNAs of SOBP were predicted by starBase, TargetScan, and miRBD, then, LncBase was applied to predict the upstream lncRNAs of miRNAs. [25][26][27][28] The common miR-NAs were achieved by VENNY 2.1.0. To verify the relationship between miRNA and SOBP, we conducted correlation and survival analysis of common miRNAs in TCGA, starBase, and KM, individually. Similar methods were used to confirm the relationship between lncRNAs and hmiR-378d. ImmLnc (http://bio-bigda ta.hrbmu.edu.cn), an online tool for identifying lncRNA regulators of immune-related pathways, was using to explore the immunological relationship of MEG8. p-value <0.05 was considered as statistically significant.

| Identification of common DEGs
After the analyses of GSE36668, GSE12470, GSE14407, and GSE27651 data sets, a total of 126 OC samples (94 tumor tissues and 32 normal tissues) were explored. There were 1825 DEGs screened from the GSE36668 data set, including 1066 upregulated and 759 downregulated genes; 3032 DEGs separated from the GSE12470 data set, including 1810 upregulated and 1222 downregulated genes; 3620 DEGs selected from the GSE14407 data set, including 2324 upregulated and 1296 downregulated genes. A total of 3051 DEGs were nominated from the GSE27651 data set, including 1490 upregulated and 1561 downregulated genes. The volcano plots of DEGs among each data set are shown in Figure 1A-D. By Venn diagram analysis, 192 common DEGs in the intersection of the four data sets were identified and selected for further analysis, containing 125 upregulation and 67 downregulation genes ( Figure 1E and F).

| Functional analysis of common DEGs
For a deeper understanding of these common DEGs, DAVID was applied to run GO and KEGG pathway enrichment analysis. GO analysis showed common DEGs were enriched in biological processes and cellular components (Figure 2A   and ptr01230 (Biosynthesis of amino acids). The detailed statistics are presented in Table 2.

| Survival analysis
To identify the key genes more accurately in OC, we explored the prognosis of the common DEGs by using GEPIA and KM, respectively. After the combined analyses of the expression and survival outcomes, we presented five upregulated key genes (RIPK4, SCNN1A, SLC4A11, ELF3, and CLDN4) were not only increased in OC samples, but also indicated poor overall survival (OS) in OC patients ( Figure 3A-E and Figure 3G), and only one downregulated hub gene (SOBP) was associated with a lower OS ( Figure 3F and G). As a validation, we performed survival analysis of the six candidate genes by using KM, only SOBP was associated with favorable relapse-free survival (RFS) in both TCGA and KM data sets ( Figure 3G-I). Moreover, SOBP was also an independent survival factor for OS, through a multivariable Cox proportional hazard model including designated clinical features and gene expression from TIMER (Table 3).

| Expression validation of candidate genes
Subsequently, TCGA OC samples and GTEx normal tissues were used to approve mRNA expression of six key genes. As shown in Figure 4, all the results of six genes (RIPK4, SCNN1A, SLC4A11, ELF3, CLDN4, and SOBP) were in accordance with previous data. Meanwhile, RIPK4 and SOBP mRNA expression was associated with the tumor pathological stage, which was consistent with their survival outcomes. The upregulated five genes showed a positive correlation with each other, while the expression of SOBP was negatively correlated with ELF3, CLDN4, and SLC4A11, significantly ( Figure 4I). Representative immunohistochemistry images of six gene protein expressions in OC tissues are shown in Figure 5.

| Correlation between SOBP expression and infiltration level of immune cells
Our findings showed that SOBP was a hub gene in OC in combination with the results of mRNA, protein expression levels, and prognostic assessment. Furthermore, we analyzed the relationship between SOBP expression and immune cell infiltration via the TIMER database. We performed a pancancerous analysis of mRNA expression by using Oncomine and TIMER databases, the results showed SOBP was a lower expression in various human cancers, including OC ( Figure 6A and B). We then noticed that low SOBP expression correlated with elevated infiltration of generally immune cell types in OC and showed a significant correlation with tumor purity. Specifically, the lower expression level of SOBP was correlated with the high infiltration of CD8+ T cells (r = −0.155, p = 6.56e-04), Macrophages (r = −0.155, p = 6.50e-04), Neutrophils (r = −0.233, p = 2.38e-07), and Dendritic cells (r = −0.204, p = 6.43e-06) in OC tissues ( Figure 6C). Next, we investigated the correlation among the transcription of SOBP and the rank of tumor-infiltrating immune cells in more detail, based on the immune marker genes expression in OC tissues from the TIMER databases. The results showed that the mRNA of SOBP was inversely related to the expression of marker genes of designated immune cells, including CD8+ T cell, T cell (general), Monocyte, TAM (tumor-associated macrophage), M1 Macrophage, M2 Macrophage, Neutrophils, Natural killer cells, Dendritic cells, Th 1 (T helper cell 1), Treg (regulatory T cell), and T cell exhaustion (Table 4).

| Construction of MEG8/miR-378d/ SOBP triple subnetwork
By a series of upstream prediction and correlation evaluations in expression, a brand-new lncRNA-miRNA-mRNA triple regulatory axis in OC was built. First, we predicted upstream miRNAs of SOBP in starBase, TargetScan, and miRBD separately and got 61 common miRNAs ( Figure 8A). After approval of expression and survival analysis, 10 key miRNAs (hsa-miR-10b-5p, hsa-miR-15a-5p, hsa-miR-16-5p, hsa-miR-17-5p, hsa-miR-200a-3p, hsa-miR-141-3p, hsa-miR-378a-3p, hsa-miR-378c, hsa-miR-378d, and hsa-miR-664b-3p) were negatively correlated with the expression of SOBP in TCGA samples ( Figure 8B, 8D-M), and only miR-378d indicated poor prognosis ( Figure 8C). Subsequently, we found the 115 upstream lncRNAs of miR-378d from the LncBase database, and draw the Venn diagram with 332 Negative Expressed LncRNAs (NELs) for miR-378d from TCGA analysis, only two lncRNAs (MEG8 and CECR7) showed predictability and negative correlation ( Figure 9A). Equally, we verified the expression levels and prognostic values of two key lncRNAs. The results showed both of them were a lower expression in OC and negatively correlated with miR-378d from GEPIA2 and starBase database, only the expression of lncRNA-MEG8 was related to better OS in KM and higher SOBP in TCGA ( Figure 9B-G). Altogether, the MEG8/miR-378d/SOBP competing endogenous RNA (ceRNA) axis could be the potential mechanism in affecting tumor progression and prognosis of OC. Finally, several immune pathways have been found to be associated with lncRNA-MEG8 by using ImmLnc, including antigen processing and presentation, cytokine receptors, TGFb family member, cytokines, and natural killer cell cytotoxicity (Table 5). In particular, the pathway of cytokines also showed significant correlations with SOBP in the previous GESA report. Therefore, the cytokines pathway must be the most probable immune response mechanism of the MEG8/ miR-378d/SOBP axis in OC.

| DISCUSSION
Despite continued advances in treatment for OC, only 19% of patients were diagnosed at an early stage due to the lack and treatment of OC is a matter of urgency. In this study, we integrated the four expression data sets (GSE36668, GSE12470, GSE14407, and GSE27651) and obtained 192 common DEGs. GO analysis showed that common DEGs were enriched in cartilage development and negative regulation of epithelial cell differentiation, palate development, negative regulation of keratinocyte proliferation, extracellular exosome, and nucleus. KEGG pathways showed that the pathways regarding the abovementioned DEGs were significantly associated with the biological behavior of OC, including PI3 K-Akt signaling pathway, Pathways in cancer, p53 signaling pathway, Rap1 signaling pathway, Cell cycle, Alanine, aspartate and glutamate metabolism, Melanoma, and Biosynthesis of amino acids. After systematic expression and survival analysis, six genes (RIPK4, SCNN1A, SLC4A11, ELF3, CLDN4, and SOBP) were recognized as key genes that may be associated with the development of OC. Expression levels of RIPK4,  SCNN1A, SLC4A11, ELF3, and CLDN4 were amplified and their upregulation was linked to poor prognosis of patients with OC. The five genes have been described to act as oncogenes in various human tumors and link to tumor progression. Moreover, several studies have also confirmed that the five genes may provide capable biomarkers for cancer. For example, RIPK4 expression is upregulated in nasopharyngeal, pancreatic, bladder urothelial, and cervical squamous cell carcinoma, and is associated with a poor outcome. [32][33][34][35] Increased SCNN1A contributed to unfavorite prognosis in pulmonary adenocarcinoma and OC. 36,37 Qin et al. also showed a high expression of SLC4A11 was an independent prognostic factor for poor OS in grade 3 or 4 serous OC through bioinformatics analysis. 38 ELF3 overexpression is significantly linked to poor outcomes in hepatocellular, colorectal cancer, and lung adenocarcinoma patients, and that ELF3 enhance cell growth, migration in these cancers. [39][40][41][42] Based on the previous studies, CLDN4 showed high expression in many epithelial malignant tumors, such as ovarian and pancreatic cancer. [43][44][45][46] All these studies including ours showed that RIPK4, SCNN1A, SLC4A11, ELF3, and CLDN4 genes may be the five key oncogenes in the progression of OC. Regarding SOBP (sine oculis-binding protein) gene (also named JXC1), it has rarely been reported in human diseases and has never been studied in tumors. Chen's results showed SOBP is crucial for cochlear growth, cell fate, and patterning of the organ of Corti. 47 Birk et al. indicated that SOBP is altered in intellectual disability and is overexpression in the brain limbic system. 48 In this study, we creatively explored the role of SOBP in OC patients. The results confirmed that SOBP was a low expression in mRNA and protein levels in OC tissues, and the downregulated of SOBP was related to poor OS and RFS in OC patients. Hence, SOBP could be considered a prognostic marker of OC, and potential mechanisms continue to be revealed. Recently, the microenvironment including several kinds of cell bunches (such as tumor cells, immune cells, and fibroblasts) in the tumor has been a hot topic. According to the mRNA expression level of SOBP, we obtained nine types of tumor-infiltrating lymphocytes in OC tissues, containing B cells, CD8+ T cells, Th 1, Treg, T cell exhaustion, macrophages (TAM, M1, and M2), neutrophils, dendritic, and natural killer cells from TIMER. The infiltration of these TILs was discovered to be inversely correlated with SOBP. Meanwhile, we identified 10 key upregulated pathways about the immune response in the low SOBP group, such as IL 5 signaling pathway, CTL (Cytotoxic T lymphocytes)mediated immune response against target cells, IL 17 signaling pathway, NO2-dependent IL 12 Pathway in NK cells, Dendritic cells in regulating TH1 and TH2 development, T helper cell surface molecules, B lymphocyte cell surface molecules, T cytotoxic cell surface molecules, Th1/Th2 differentiation, and Cytokine network. These findings suggested that SOBP expression was strictly related to the level of TILs in OC tissues, and further analyses will focus on detailed mechanisms of SOBP modifying these infiltrating immune cells in OC. The intersection between SOBP expression and  immune functions in OC implies a potential target of SOBP for the next treatment of OC patients. Meanwhile, the roles of lncRNAs have drawn more and more attention. 49 Several studies indicated that lncRNAs may act as competitive endogenous RNA (ceRNAs) in a new triple regulatory network by competitively binding their shared miRNAs. [50][51][52][53] However, neither lncRNAs nor miRNAs have been described to engage in the direct modulation of SOBP in OC. In this work, we noticed the existence of various lncRNAs and miRNAs which may lead to aberrant expression of SOBP. Furthermore, by correlation of expression and prognosis, we assumed a regulatory axis of lncRNA-miRNA-mRNA must affect the development and progression of OC. In our study, three databases (starBase, TargetScan, and miRBD) were used to predict miR-378d as the upstream miRNA of SOBP simultaneously. The negative correlation occurred between the expression of miR-378d and SOBP, and the high miR-378d was associated with poor prognostic effect in OC. Similarly, lncRNA-MEG8 could be regarded as the upstream of miR-378d with a negative correlation in expression, and the lower MEG8 was bound up with an unfavorable prognosis in OC. Even more, the same related immune pathway, the cytokines pathway, was found among MEG8 and SOBP. In previous studies, miR-378d was upregulated in colorectal cancer and might be a potential biomarker in diagnosis. 54,55 MEG8 suppresses activation of the epithelial-mesenchymal transition of hepatocytes by the Notch pathway. 56 Therefore, under a comprehensive analysis and previous reports, we establish the MEG8/miR-378d/SOBP axis and speculate that it may act a crucial role in the tumor progression and immune function of OC by regulating the cytokines pathway.
In general, we studied the differential expression of mRNA and protein levels in normal ovarian and tumor tissue samples. In addition, the key gene SOBP was creatively screened out for verification by combining the results of multiple databases. In addition, we reported the novel MEG8/miR-378d/SOBP ceRNA subnetwork first, and speculate that it may act a crucial part in the tumor progression and immune function of OC by regulating the cytokines pathway. Thus, SOBP could offer great clinical value in the diagnosis, prognosis evaluation, and immunotherapy of OC patients.
There were several limitations in this study. First, we explored the expression of SOBP in OC using multiple public databases, and we did not prove these findings in our clinical samples. Second, further studies should be conducted to support the consistency with protein level and causality relation between different variables using western blot. Third, because this was a retrospective study, the level of evidence is imperfect. Although we creatively launched a novel lncRNA-miRNA-mRNA subnetwork in OC and predicted potential immune-related pathways implicated in its function, more studies are warranted to elucidate the specific roles of MEG8/miR-378d/SOBP and underlying mechanisms that affect tumor development and progression.

| CONCLUSIONS
In summary, we detected six DEGs (RIPK4, SCNN1A, SLC4A11, ELF3, CLDN4, and SOBP) which can serve as tumor markers for the diagnosis and prognosis of OC. In particular, the gene SOBP was significantly low expressed in OC and negatively correlated with TILs. Moreover, a novel MEG8/miR-378d/SOBP axis probably affects the progression of OC by regulating the cytokines pathway and may be utilized as a potential target of immunotherapy in the future. Even so, further biological function validation is required to reinforce this conclusion.