Identification of the hub genes RUNX2 and FN1 in gastric cancer

Abstract Background This study identified key genes in gastric cancer (GC) based on the mRNA microarray GSE19826 from the Gene Expression Omnibus (GEO) database and preliminarily explored the relationships among the key genes. Methods Differentially expressed genes (DEGs) were obtained using the GEO2R tool. The functions and pathway enrichment of the DEGs were analyzed using the Enrichr database. Protein–protein interactions (PPIs) were established by STRING. A lentiviral vector was constructed to silence RUNX2 expression in MGC-803 cells. The expression levels of RUNX2 and FN1 were measured. The influences of RUNX2 and FN1 on overall survival (OS) were determined using the Kaplan–Meier plotter online tool. Results In total, 69 upregulated and 65 downregulated genes were identified. Based on the PPI network of the DEGs, 20 genes were considered hub genes. RUNX2 silencing significantly downregulated the FN1 expression in MGC-803 cells. High expression of RUNX2 and low expression of FN1 were associated with long survival time in diffuse, poorly differentiated, and lymph node-positive GC. Conclusion High RUNX2 and FN1 expression were associated with poor OS in patients with GC. RUNX2 can negatively regulate the secretion of FN1, and both genes may serve as promising targets for GC treatment.


Introduction
Gastric cancer (GC) is a common malignant digestive tract tumor with a poor prognosis. GC incidence ranks fifth globally for mortality of all malignancies and third for cancer-related mortality of all malignancies, and the majority of GC cases occur in the developing countries [1][2][3]. In China, GC is one of the most common malignancies and was the second leading cause of cancer-related death in males and the fifth leading cause of cancer-related death in females in 2014 [4]. In the Western world, the majority of patients with GC are diagnosed at an advanced stage, so radical surgical treatment is no longer possible at the time of diagnosis, leading to poor prognosis. The expected 5-year survival for patients is approximately 26% in Western countries. Surgical resection is still the only curative therapy for nonmetastatic GC, while neoadjuvant chemotherapy, adjuvant chemotherapy chemoradiotherapy, and targeted therapy can reduce recurrence and prolong survival [5][6][7].
Complex signaling molecules lead to poor prognosis by affecting the occurrence and development of tumors [8]. Currently, the application of high-throughput platforms in gene expression analysis is becoming more valuable in many aspects of clinical research [9][10][11]. Through second-generation sequencing technology, hundreds of differentially expressed genes (DEGs) have been discovered and are involved in different tumor cell biological processes; the relationships between these genes must be very complicated. RUNX2 is an evolutionarily conserved regulator of cell fate, which contains a highly conserved runt domain. RUNX2 can recognize and directly bind to the consensus sequences (TGTGGT or ACCACA) and mediate transcriptional activation or repression of target genes [12]. Boregowda et al. showed that RUNX2 is overexpressed in melanoma cells and mediates their migration and invasion by affecting the expression of focal adhesion kinase [13]; Cao et al. reported that RUNX2 was associated with EMT by regulating galectin-3 expression in hepatocellular carcinoma [14]. RUNX2 has been shown to be involved in many molecular mechanisms of tumor progression [15][16][17]. Although Guo et al. found that the transcription factor RUNX2 upregulates the chemokine receptor CXCR4 to promote the invasive and metastatic potential of human GC [18], research on RUNX2 in GC has been rarely been conducted. The specific mechanism underlying the participation of RUNX2 in GC pathogenesis remains unclear. Furthermore, studies have shown that FN1 is involved in tumor invasion and metastasis [19,20]. However, no reports on the relationship between RUNX2 and FN1 in GC have been reported until now.
In this study, we aimed to screen the hub genes (genes with a high degree of connections in the gene expression network, which do not involve betweenness) in the development of GC. By constructing a PPI network and using cytoHubba, we found that RUNX2 and FN1 are hub genes in DEGs, that FN1 ranks number 1, and that these genes are connected. RUNX2 is a transcription factor, and we found that FN1 may be a target gene of the transcription factor RUNX2 through the online website Harmonizome (http://amp.pharm.mssm.edu/Harmonizome/). Therefore, we decided to inhibit the expression of RUNX2 in MGC803 cells to observe the changes in FN1 by transfection with lentiviral vectors. Finally, we obtained the prognostic value (the association of the protein level with the survival time of patients with GC) of RUNX2 and FN1 in GC from the KM plotter (http://kmplot.com/analysis/).

Microarray data and identification of DEGs
The dataset GSE19826 is based on the platform of GPL 570 (Affymetrix Human Genome U133 Plus 2.0 Array) and contains 12 GC tissue, 12 noncancer tissue, and three normal gastric tissue samples from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/ geo). GEO2R, supplied by NCBI (https://www.ncbi.nlm. nih.gov/geo/geo2r/), is an interactive web tool that was used to identify the DEGs between GC tissue samples and noncancerous tissue samples. GEO2R performs comparisons on original submitter-supplied processed data tables using the GEOquery and limma R packages from the Bioconductor project [21]. In the current study, genes with |log FC| ≥ 1.5 and p < 0.05 (adjusted by the false discovery rate) were regarded as DEGs [22]. Also, the GSE79973 dataset was used to validate the findings.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG)
To elucidate the potential biological processes, we performed GO enrichment analysis utilizing the clusterProfiler package in R software. KEGG pathway enrichment analysis was also carried out by the clusterProfiler package in R software to identify promising signaling pathways correlated with the DEGs. An adjusted p-value < 0.05 was considered significant.

Protein-protein interaction (PPI) network building
STRING database version 10.5 (http://string-db.org) was applied to evaluate the PPI information. DEGs were uploaded into STRING, a confidence score >0.9 defined them as significant [23], and then, the PPI network was constructed. We downloaded the resulting data in a table in a tab-separating values format. The data obtained were uploaded into Cytoscape software, which was used to construct PPI relationship subnetworks. The cytoHubba plug-in was used to identify the hub genes.
A hub gene is defined as a gene with the highest degree of connectivity in the hub module [23]. The hub genes showed a strong association with other node proteins (more than 10) [23].

Cell culture
The human GC cell line MGC-803 was obtained from the Department of General Surgery, Beijing Friendship Hospital (Beijing, China). The cells were maintained in the Dulbecco's modified Eagle's medium (HyClone, Logan, UT, USA) containing 10% FBS (HyClone, Logan, UT, USA), 100 mg/mL penicillin, and 100 mg/mL streptomycin. Cells were cultured at 37℃ in a humidified atmosphere with 5% CO 2 .

The prognostic value of RUNX2 and FN1
The prognostic value of RUNX2 and FN1 was obtained from the KM plotter (http://kmplot.com/analysis/).

Statistical analysis
All statistical analyses were performed using GraphPad Prism 6.0 software (GraphPad Software, Inc., USA). Data are presented as the mean ± standard deviation. Student's t-test was used to compare the differences between the two groups. Differences were considered statistically significant at a p value of < 0.05.

Identification of DEGs
A total of 134 DEGs were identified from the GSE19826 dataset and selected for further analysis. There are 69 upregulated genes and 65 downregulated genes in GC tissues when compared with those in noncancer tissues ( Table 1). Interestingly, the log FC value of RUNX2 was +1.87, and the log FC value of FN1 was −1.72, and a negative correlation might exist between them. The dataset GSE79973 further validated the findings (Supplementary Table 1).

PPI network construction
The
structure organization, collagen fibril organization, and connective tissue development (Figure 2a). KEGG pathway analysis showed that the DEGs were mainly enriched in protein digestion and absorption, ECMreceptor interaction, and focal adhesion (Figure 2b). FN1 is involved in an extracellular matrix organization and endodermal cell differentiation in BP term enrichment analysis and ECM-receptor interaction, focal adhesion, amoebiasis, AGE-RAGE signaling pathway in diabetic complications, human papillomavirus infection, and PI3K-Akt signaling pathway in the KEGG pathway analysis (Figure 2). RUNX2 was categorized in ossification, osteoblast differentiation, endochondral ossification, and osteoblast fate commitment in BP term enrichment analysis and parathyroid hormone synthesis, secretion, and action and transcriptional misregulation in cancer in the KEGG pathway analysis (not shown in Figure 2 due to a relatively low score).

Downregulation of RUNX2 increased FN1 expression in MGC803 cells
We found 20 hub genes in the GSE19826 dataset. Both RUNX2 and FN1 are hub genes, and FN1 ranks No. 1 with the highest score of the 20 hub genes; thus, we concluded that FN1 may play an important role in GC.
Moreover, studies about the hub gene RUNX2 in GC are lacking, so we are very interested in its role in GC. Taken together, we decided to conduct a preliminary study of the correlation between the two. To explore the changes in FN1 in GC cells with RUNX2, we transfected MGC803 cells with shRUNX2 and Scramble with an MOI of 10 in the presence of 6 µg/mL polybrene (Figure 3a). qRT-PCR was used to test the transfection efficacy of shRUNX2#1, shRUNX2#2, and shRUNX2#3; Lentivirus-shRUNX2#1 had the best transfection efficacy; thus, we chose MGC-803/shRUNX2#1 for further experiments (Figure 3b). qRT-PCR and Western blot analyses confirmed RUNX2 downregulation in MGC803 cells (p < 0.0001), and the expression of FN1 in MGC803 cells increased significantly   with the downregulation of RUNX2 (p < 0.0001) (Figure 3c-e).

The prognostic value of RUNX2 and FN1
The prognostic value of RUNX2 and FN1 was obtained from the KM plotter (http://kmplot.com/analysis/). According to the expression of a particular gene, the patients with GC were split into two groups with high and low expression. The overall survival (OS) of patients with GC was evaluated using a KM plot. The hazard ratio (HR) with 95% confidence intervals and log-rank P values are shown on the web page. As shown in Figure 4a and b, GC patients with higher expression of RUNX2 had a lower probability of survival than those with lower RUNX2 expression. However, those with higher expression of FN1 also had a lower probability of survival, which contradicted the results obtained in the cell experiments. Then, subgroup (those with diffuse, poorly differentiated, and lymph nodepositive GC) analysis was performed. As shown in Figure 4c and d, the median survival of patients with higher expression of RUNX2 was significantly longer than those with lower expression, whereas patients with lower expression of FN1 had a significantly longer median survival time than those with higher expression of FN1.

Discussion
GC is still a disease with high rates of morbidity and mortality due to its high heterogeneity. Although surgery is the main treatment method, the other treatments, including chemotherapy, radiotherapy, targeted therapy, and gene therapy, are available; however, the 5-year survival rate was still less than 30% [24]. GC is a main contributor to the global burden of disabilityadjusted life-years from cancer in men [25]. The burden of GC remains very high in Asia, Latin America, and Central and Eastern Europe, whereas it is no longer a common cancer in North America and most Western European countries [26]. It is essential to explore the mechanisms of GC progression to prevent its occurrence, improve treatment efficacy, and improve the survival rate of patients with GC.
In this study, 134 DEGs were screened, consisted of 69 upregulated genes and 65 downregulated genes. These genes were mainly enriched in extracellular matrix organization, collagen fibril organization, and protein complex subunit organization in biological processes. Among these DEGs, 20 genes had high degrees in the PPI network, including FN1, COL1A1, COL1A2, BGN, SPARC, and RUNX2. FN1 ranked first in the hub genes. The hub genes we found are partly the same as the hub genes found in other studies [22,27]. However, there were inconsistencies. The reason may be that multiple data sets intersected in some studies, which will inevitably lead to the loss of some hub genes [22,27].
As a member of the RUNX family, RUNX2 is identified by the runt-homology domain. In contrast to RUNX1 and RUNX3, whose mutations are closely linked to the promotion of leukemia and GC, respectively [28,29], the initial studies suggest that RUNX2 acts as a master regulator of osteoblast differentiation and bone development. During bone development, RUNX2 facilitates the differentiation of mesenchymal stem cells into osteoblast lineage cells [30]. Some studies suggest that RUNX2 has the ability to transactivate its downstream target genes, such as MMP9, MMP13, VEGF, survivin, and IL-8, which are involved in tumor progression, invasion, and metastasis [31][32][33][34][35][36]. FN1 is a member of the ligand glycoprotein family, which is widely expressed in various cell types and is involved in cell adhesion and migration [37]. As a transcription factor, the studies of RUNX2 in GC were relatively less, and we have a keen interest in it. FN1 ranks first in the DEGs of GSE19826, so we decided to conduct a preliminary exploration of the correlation between RUNX2 and FN1. We detected that the expression level of FN1 was upregulated in MGC803 cells by RT-qPCR and Western blot analyses. Therefore, we initially concluded that RUNX2 can act negatively on the FN1 gene in GC. This result seems inconsistent with the results of the prognostic value of the two genes obtained from the KM plotter. We performed a subgroup analysis of the two genes in the same datasets; to our surprise, we found that among the patients with diffuse and poorly differentiated GC, those with higher expression of RUNX2 had a significantly higher median survival time than those with lower RUNX2 expression, whereas patients with lower expression of FN1 had a significantly longer median survival time than those with higher FN1 expression (Figure 4c and d). Based on our findings, we speculate that the expression of RUNX2 in GC tissue changes with the degree of tumor malignancy, and its effect on FN1 also changes, and with a high degree of malignancy, the two genes are negatively correlated. To date, the role of RUNX2 in carcinogenesis and cancer progression remains unclear and may be tissue-and context-dependent [38]. RUNX2 may play a promoting role in invasive bone cancer [39], prostate cancer [31], pancreatic cancer [40], and GC [18]. Nevertheless, Chimge et al. showed that RUNX2 may possess tumor suppressor properties in breast cancer [41]. The overexpression of RUNX2 leads to the upregulation of Bax and increased sensitivity to apoptosis [42]. FN1 can downregulate P53 and inhibit apoptosis in colorectal cancer [43]. Considering that apoptosis is one of the mechanisms of drug resistance [44], we hypothesized that the mechanism by which RUNX2 negatively regulates FN1 may be involved in chemoresistance in GC.
This preliminary study was based on bioinformatics analyses with only single-cell experiments and with only qRT-PCR and Western blot analyses to verify the relationship between RUNX2 and FN1. The specific roles of the two hub genes in GC need to be verified by further research. In the future, more cell lines will be added, and clinical samples will be combined to further explore the mechanisms of two hub genes involved in the invasion, metastasis, and drug resistance of GC.
In summary, FN1, COL1A1, COL1A2, COL3A1, BGN, COL5A2, THBS2, SPARC, FBN1, COL5A1, SPP1, COL6A3, TIMP1, SERPINH1, COL12A1, RUNX2, BMP1, COL10A1, NID2, and COL8A1 might be hub genes of GC. The transcription factor RUNX2 may act negatively on the FN1 gene, and they might be correlated with the prognosis of GC. This study might shed new light on the correlation between RUNX2 and FN1, which may be promising therapeutic targets for GC treatment. Further investigation is needed to elucidate the specific biological mechanisms of RUNX2 and FN1 in GC in the future, and the relationships among the other hub genes need to be clarified with additional research.