Introduction
The runt-related transcription factor (RUNX) family includes RUNX1, RUNX2 and RUNX3, which are indispensable for normal lineage-specific development.1 Among which, RUNX2 is necessary during embryogenesis for skeletal development and is a key regulator for morphogenesis of organs like thyroid and breast.2–4 Moreover, RUNX2 exerts its biological functions via the activation of downstream signaling pathways such as epithelial mesenchymal transition (EMT), cell cycle, cell apoptosis, NOTCH, WNT, RAS and transforming growth factor-β (TGF-β)/Smad, etc.5,6
RUNX2 is increasingly implicated in tumor biology for potential tumorigenicity. It was reported that enhanced RUNX2 could transcriptionally induce gene expression involved in tumor initiation, development and metastasis.7 For example, uncontrolled RUNX2 expression was observed in the carcinogenic process.5 In melanoma, RUNX2 participated in the regulation of the EMT process.8 Upregulation of RUNX2 expression level could recover the cancer cell invasion in cervical cancer.9 Additionally, in the late-stage thyroid, breast and prostate cancer, upregulated expression of RUNX2 could accelerate bone metastasis.10,11 These revealed that RUNX2 could function as an oncogene and emerged as a promising therapeutic target to treat cancer. Therefore, extensive work has been done. It was confirmed that in cervical cancer, downregulated RUNX2 could slow the tumor growth and metastasis through inhibiting the cancer cell proliferation, invasion, migration, and reversing the EMT process.12 Suppressing the expression of RUNX2 in metastatic prostate and breast cancer cells could extremely decrease tumor progression and metastasis in vivo.11,13 Knockdown of RUNX2 in esophageal carcinoma cells could obviously inhibit cell growth, migration and invasion, and increased cell apoptosis and suppressed tumor growth.14
Nearly 25% of all cancer-related morbidity was contributed to lung cancer (LC) in 2020.15 Even more dire, greater than half of all LC patients have metastasized when diagnosed (57%), which may explain the unsatisfactory 5-year survival rate (5%).15 Therefore, effective and strong prognostic biomarkers and potential drug targets are urgently needed to improve diagnosis and prognosis, and personalized drug treatments of LC.
The dysregulated expression level of RUNX2 and its relationship with LC have been partly reported. It was demonstrated that RUNX2 expression level was upregulated in lung squamous cell carcinoma (LUSC) and lung adenocarcinoma (LUAD) tissues.16,17 Moreover, RUNX2 overexpression could promote EMT and increase the expression of EMT related genes including TWIST1, SNAIL1 and VIMENTIN, which directly reflect an increased migratory and invasion capacity of LUAD cells.18 These results suggested that RUNX2 could serve as a candidate target for LC therapy. Actually, previous research has reported that targeting RUNX2 could inhibit LC cell growth.19 Downregulation RUNX2 expression could decrease non-small cell lung cancer (NSCLC) tumor volume and weight.20 It was also reported that the proliferation, migration, and invasion of human LC cells could be repressed by downregulating the expression level of RUNX2.21,22
Numerous RNA and DNA research indicated their critical role in determining the biological and biomedical functions with the development of second-generation sequencing technology.23 Nowadays bioinformatic analyses has become one of the hottest research priorities in medicine. Specifically, bioinformatic analyses were carried out to predict gene targets and analyze their potential roles in various cancers.23 However, biological information analysis has not been be applied comprehensively to explore the potential role of RUNX2 in LC. Based on the publicly available gene information, containing gene expression, gene variation and epigenetic mutations, etc., we aimed to identify the expression pattern and variation of RUNX2 in LC and explore its prognostic value and potential functions in LC. Unlike previous studies that explored the function of RUNX2 through experimental methods, the bioinformatics analyses performed in our study will provide a comprehensive evaluation of RUNX2 in LC.
Materials and methods
ONCOMINE database analysis
ONCOMINE is a database of high-throughput sequencing data of cancers (www.oncomine.org ), with which we could assess and visualize the differential expression of specific genes from all availably public datasets and conduct differential expression analyses. ONCOMINE was utilized to compare the RUNX2 mRNA expression difference between clinical cancer samples and corresponding normal tissue samples. A student’s t-test was used to calculate the significant expression differences between tumor and normal tissues. The cut-off values were set at a P-value of 0.05 and fold change of 2.
Gene Expression Profiling Interactive Analysis (GEPIA)
The online database Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/ ) is a newly developed interactive website to analyze the RNA sequencing expression data of 9,736 tumors and 8,587 normal samples from the Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) projects.24 It was used to find the expression differences of RUNX2 between LUAD, LUSC and corresponding normal lung samples. The threshold was set as following: P-value of 0.01, fold change of 2, matched TCGA normal and GTEx data.
Survival prediction analysis of lung cancer patients
The Kaplan-Meier Plotter database (www.kmplot.com ) is one of the most effective analytical tools to evaluate survival and prognosis information.25 The prognosis value of RUNX2 was predicted in patients bearing LC by Kaplan-Meier Plots, which provides an assessment of gene expression data and survival information of 2,437 LC patients. The survival curve of progression-free survival (FP), overall survival (OS) and post progression survival (PPS) in patients with LUAD and LUSC was generated using the JetSet best probe set with the median as the cutoff. Log-rank test results with P-value < 0.05 were considered statistically significant.
Immunohistochemistry validation by Human Protein Atlas
The Human Protein Atlas (HPA) (https://www.proteinatlas.org/ ) is an extensive database for exploration of the expression pattern of mRNA and protein in cancer and normal tissue from a
large volume of proteomics and transcriptomics data.26 Therefore, the cancer and normal tissue atlas of the HPA were used to generate a list of cancer types detailing the gene and protein expression of RUNX2.
RNA extraction and real-time quantitative polymerase chain reaction (RT-qPCR)
The human LUSC cell line SK-MES-1, LUAD cell line A549 and normal lung bronchus epithelial cell line HBE were obtained from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). HBE (control group) was cultured with dulbecco's modified eagle's medium (DMEM) (Hyclone, Logan City, UT) and SK-MES-1, A549 (experimental group) were cultured with RPMI1640 containing 10% fetal bovine serum (FBS, Gibco, New Zealand), 1% antibiotic-antimycotic (Gibco, Grand Island, NY) at a density of 3×104/cm2 onto the 10-cm culture dish for 24 hours. Then total RNA was extracted from the collected cells using Trizol Reagent (Invitrogen, Carlsbad, CA), and cDNA was converted from 2 µg total RNA using 200 U of reverse transcriptase (MMLV RT) and 20 pM oligo dT (Promega Corporation, Madison, WI, USA) at 30 °C for 10 min, 42 °C for 60 min, 95 °C for 5 min, and 5 °C for 5 min. The forward primer (5-TCTGGCCTTCCACTCTCAGT-3) and reverse primer (5-TATGGAGTGCTGCTGGTCTG-3) were synthesized based on the Runx2 mRNA sequence. Quantitative PCR was performed with a protocol as follows: initial denaturation at 95 °C for 10 min, followed by 40 cycles at 95 °C for 30 s, 58 °C and 72 °C for 45 s, using the Strata Gene Mx3000p (Agilent Technologies, Inc., Santa Clara, CA, USA). Each gene of interest was normalized to GAPDH and the fold change was compared relative to the control sample. Each assay was performed in triplicate. The statistical analyses were processed using the statistical software Statistical Package for the Social Sciences (version 19.0, SPSS, Chicago). The data are represented as the mean ± standard derivation. A student's t-test was employed to analyze the differences between the control and experimental group, with p < 0.05 considered statistically significant.
cBioPortal analysis
The LUAD (TCGA, Firehose Legacy) and LUSC (TCGA, Firehose Legacy) datasets containing 517 and 501 pathologically reported cases, respectively, were selected for further analyses of RUNX2 using the cBioPortal (https://www.cbioportal.org/ ). The spectrum of genomic alterations included mutations, putative copy-number alterations (CNA) from GISTIC, mRNA expression z-scores (RNA Seq V2 RSEM) and protein expression Z-scores (RPPA).
Construction of a protein–protein interactions (PPI) network
The STRING database (https://string-db . org/cgi/input.pl) aims to collect, score, and integrate all publicly available sources of PPI information and to complement these with computational predictions. RUNX2 and the targets closely related to RUNX2 alterations (top 100) were put into STRING to build the PPI network to achieve a comprehensive and objective global network, including direct (physical) as well as indirect (functional) interactions. The cutoff of the PPI confidence score was 0.7 and the disconnected nodes were hidden.
Gene ontology analysis
DAVID (http://david.abcc.ncifcrf.gov/ ) is a database for annotation, visualization and integrated discovery.27 Gene Ontology (GO) is a widely confirmed gene functional enrichment database which is generally used to explore the enriched GO terms.28 GO functional annotation enrichment analysis could predict the target gene function through three aspects including biological process (BP), cellular component (CC) and molecular function (MF). Therefore, GO enrichment analysis was conducted to identify functional categories of RUNX2 and the genes closely associated with RUNX2 alterations. Enriched GO terms were identified according to the cut-off criterion of P-value < 0.05 and gene counts more than 2. The top 20 GO terms of BP, CC and MF were selected to display.
Results
The transcriptional level of RUNX2 in multiple types of tumors
The mRNA expression level of RUNX2 in multiple types of cancers was compared with that in the normal samples using the ONCOMINE database (Fig. 1). As shown in Figure 1a, the expression level of RUNX2 was obviously increased in LC patients of three datasets. In Garber’s dataset,29 RUNX2 was upregulated in all LC types compared to normal samples, including LUAD with a fold change of 2.829 (Fig. 1b), LUSC with a fold change of 3.703 (Fig. 1c), large-cell lung carcinoma with a fold change of 3.051 (Fig. 1d) and small cell lung carcinoma with a fold change of 1.768 (Fig. 1e).
Association between RUNX2 expression and the clinicopathological parameters of lung cancer patients
The mRNA expression level of RUNX2 between LC and normal lung tissues was evaluated by GEPIA. The results indicated that RUNX2 expression level was significantly increased in LUAD and LUSC tissues than that in the lung tissues (Fig. 2a–b). The relationship of RUNX2 expression with LUAD and LUSC tumor stage was also analyzed, however, no significant difference was mined (Fig. 2c–d).
Increased RUNX2 expression indicated a poor prognosis of lung cancer patients
The prognostic role of RUNX2 in patients with LC was investigated using Kaplan–Meier plots. The results indicated that the upregulated RUNX2 mRNA expression level was negatively correlated with the FP and OS in all LC patients (Fig. 2e–f, p < 0.001). Despite a decreasing trend of PPS in patients with increased expression of RUNX2, there was no statistically significant difference (Fig. 2g, p > 0.05).
Increased RUNX2 mRNA and protein levels validated by HPA and RT-qPCR
The HPA database was utilized to identify the expression level of RUNX2 in different normal and cancer tissues. The results demonstrated that NK-cells, dendritic cells, T-cells, cervix uterine, tonsil, appendix, bone marrow, granulocytes, lymph node, spleen contained the top ten normal tissue expression of RUNX2 mRNA (Fig. 3a). In addition, the seven tissues achieved by comparison of RUNX2 protein expression level, were salivary gland, spleen, lymph node, tonsil, bone marrow, nasopharynx, and bronchus (Fig. 3b). Among all types of cancer, thyroid cancer, pancreatic cancer, breast cancer and LC showed the highest expression level of RUNX2 mRNA and urothelial cancer, LC, skin cancer and lymphoma presented the highest expression value of RUNX2 protein (Fig. 3c–d). Our finding demonstrated that RUNX2 mRNA and protein expression levels were both obviously increased in LC. Moreover, the increased protein expression level of RUNX2 was also detected in LUAD and LUSC tissues compared to normal lung tissues (Fig. 3e–g). As shown in Figure 3h, the LUAD (A549) and LUSC (SK-MES-1) cells revealed significantly higher expression of RUNX2 with the level of 1.850±0.5075-fold and 1.773 ± 0.3810-fold more than that of the normal lung bronchus epithelial cells (HBE). The analysis above suggests that RUNX2 was a specificity-increased biomarker in LC.
Analysis of RUNX2 and the altered neighbor genes in lung cancer patients
The RUNX2 alterations and neighbor genes altered by RUNX2 alterations were analyzed in LUAD and LUSC using the cBioPortal database. As shown in Figure 4, RUNX2 was altered in 44 out of 517 patients with LUAD (8.51%) (Fig. 4a, c) and 42 out of 501 patients with LUSC (8.38%) (Fig. 4b, d). The top 10 neighbor genes which were most frequently altered by RUNX2 in LUAD and LUSC are listed in Tables 1 and 2, respectively, and all altered neighbor genes information is shown in Table S1–2. The results indicated that cell proliferation, invasion and metastasis-related genes DPYSL3, VCAN, ANTXR1, SULF1, COL1A2, COL3A1, COL6A3 were closely associated with RUNX2 alterations in LUAD (Table 1), and LOXL2, LTBP2, TPM1, TIMP3, COL5A1, COL6A2 were closely associated with RUNX2 alterations in LUSC (Table 2). Furthermore, the proteins in the PPI network are considered core targets in the pathogenesis of RUNX2 for LUAD (Fig. 5a) and LUSC (Fig. 5b).
Table 1The top 10 significant genes correlated with RUNX2 in lung adenocarcinoma (cBioPortal database)
| Correlated Gene | Cytoband | Spearman's Correlation | p-Value | q-Value |
|---|
| DPYSL3 | 5q32 | 0.572521 | 2.35E−46 | 4.72E−42 |
| COL1A2 | 7q21.3 | 0.558711 | 8.92E−44 | 8.97E−40 |
| ALPK2 | 18q21.31–q21.32 | 0.54582 | 1.78E−41 | 9.93E−38 |
| VCAN | 5q14.2–q14.3 | 0.545571 | 1.97E−41 | 9.93E−38 |
| SERPINF1 | 17p13.3 | 0.544201 | 3.42E−41 | 1.37E−37 |
| COL3A1 | 2q32.2 | 0.542773 | 6.04E−41 | 1.88E−37 |
| ALDH1L2 | 12q23.3 | 0.542577 | 6.53E−41 | 1.88E−37 |
| ANTXR1 | 2p13.3 | 0.537946 | 4.06E−40 | 9.32E−37 |
| SULF1 | 8q13.2–q13.3 | 0.537882 | 4.16E−40 | 9.32E−37 |
| COL6A3 | 2q37.3 | 0.535679 | 9.84E−40 | 1.98E−36 |
Table 2The top 10 significant genes correlated with RUNX2 in squamous cell lung carcinoma (cBioPortal database)
| Correlated Gene | Cytoband | Spearman's Correlation | p-Value | q-Value |
|---|
| PMEPA1 | 20q13.31 | 0.666687883 | 1.08E−65 | 2.17E−61 |
| APCDD1L | 20q13.32 | 0.647424619 | 7.46E−61 | 7.53E−57 |
| LOXL2 | 8p21.3 | 0.6412603 | 2.23E−59 | 1.50E−55 |
| COL5A1 | 9q34.3 | 0.637548012 | 1.66E−58 | 8.40E−55 |
| COL6A2 | 21q22.3 | 0.624767517 | 1.37E−55 | 5.52E−52 |
| LTBP2 | 14q24.3 | 0.622393683 | 4.60E−55 | 1.55E−51 |
| TPM1 | 15q22.2 | 0.621275378 | 8.11E−55 | 2.34E−51 |
| ADAMTS14 | 10q22.1 | 0.619570308 | 1.92E−54 | 4.84E−51 |
| TIMP3 | 22q12.3 | 0.619182766 | 2.33E−54 | 5.23E−51 |
| INHBA | 7p14.1 | 0.616864852 | 7.44E−54 | 1.50E−50 |
GO enrichment analysis
The biological function of RUNX2 and the genes closed related to RUNX2 alterations were evaluated by performing GO enrichment analysis. As shown in Figure 6a, clinically important module genes of LUAD in the BP group were strongly enriched in cell adhesion, oxidation reduction, cell motion, extracellular matrix organization, blood vessel development, cell migration, regulation of cell proliferation and collagen fibril organization, etc. The genes in the CC terms were strongly enriched in membrane, extracellular matrix, membrane fraction and cell fraction, etc. (Fig. 6b). The genes in the MF group were significantly enriched in polysaccharide binding, oxidoreductase activity, acting on peroxide as acceptor and growth factor binding, etc. (Fig. 6c). For LUSC, the clinically significant module genes in the BP group were significantly enriched in cell adhesion, extracellular matrix organization, response to wounding, vasculature development, blood vessel development, collagen metabolic process, collagen fibril organization, cell-substrate adhesion, oxidation reduction, regulation of cell motion, blood vessel morphogenesis, regulation of cell migration, etc. (Fig. 6d). The genes in the CC group were mainly enriched in collagen, extracellular matrix, cell surface, etc. (Fig. 6e). And the genes in the MF group were significantly enriched in extracellular matrix structural constituent, calcium ion binding, growth factor binding, metalloendopeptidase activity, cytokine binding, integrin binding, etc. (Fig. 6f).
Discussion
Functional genome analysis revealed that gene expression is far more complicated than expected and required an ongoing and widespread regulatory landscape.30 Tumor initiation, progression and metastasis are an outcome of deregulation of multiple signaling pathways via the abnormal function of transcription factors.
RUNX2 is a lineage-specific transcription factor which is essential in bone development.31 However, emerging evidence has indicated that RUNX2 may be a pivotal oncogene in tumorigenesis in recent years. The oncogenic potential of RUNX2 was first demonstrated in T-cell lymphoma.32 After that, a large number of experiments have since proven the involvement of RUNX2 in tumorigenesis.11,33,34 In the initial step of tumor formation, RUNX2 was overexpressed and caused abnormal activation of pathways like PI3K/Akt, Wnt, TGF-β, BMP, MAPK/ERK and Notch, MET, which were a hallmark of various types of cancers.35 RUNX2 was documented to be increased in tumor cells and had a negative prognostic value in cancers of the thyroid, breast, colon and prostate.13,36,37 For example, it was reported that upregulated levels of RUNX2 accelerated prostate cancer aggressiveness by potentiating tumor cell progression and metastasis through promoting filopodia and adhesion in cancer cells.13,33 In T cell lymphomas, highly expressed RUNX2 contributed to the increased transcriptional activity of genes that accelerate survival, growth, and metastasis in the cancer cells.38 In pituitary tumors, RUNX2 upregulated the tumor-promoter and inhibited tumor-suppressors, which facilitate bone metastasis.39 In advanced thyroid, breast and prostate cancer, increased RUNX2 could dramatically accelerate tumor metastasis.10,11 These above studies revealed that RUNX2 can act as a potential oncogene and therapeutic target in various types of cancers.
Actually, it was demonstrated that downregulated RUNX2 could inhibit tumor development via suppressing the EMT process of tumor cells and blocking cell proliferation and invasion in cervical cancer.12 Downregulated RUNX2 in metastatic prostate and breast cancer cells could significantly reduce tumor metastasis.11,13 RUNX2 knockdown could effectively increase esophageal carcinoma cell apoptosis, inhibit cell invasion and eventually suppress tumor enlargement and metastasis.14
LC is one of the most common malignancies leading to substantial mortality worldwide, because of a high proportion of LC patients diagnosed in an advanced stage.40 Therefore, new diagnostic markers and therapeutic targets for LC are urgently required. The relationship between the dysregulated expression of RUNX2 and LC has been partly reported. It was yielded that RUNX2 overexpression could promote EMT and increase the migratory and invasion capacity of LUAD cells.18 Targeting RUNX2 could inhibit LC cell growth and greatly reduce tumor volume and weight.19,20
In recent decades, bioinformatic analyses are increasingly applied to explore the relationships between genes and cancers. Although the function of RUNX2 in the tumorigenesis and the prognostic values in many cancers have been partially validated, further bioinformatics analysis of LC and RUNX2 has not been performed. The current study is a first attempt to explore the mRNA and protein expression, potential functions and distinct prognostic values of RUNX2 in LC.
RUNX2 was reported to overexpress in primary or metastatic LC samples compared to that of normal tissue samples.18 In the present study, ONCOMINE and GEPIA databases manifested the increased expression of RUNX2 in human LC compared with lung tissues (Fig. 1, Fig. 2a, b). The prognostic value of RUNX2 in LC patients was documented using the Kaplan–Meier Plotter analysis. As shown in Figure 2e–f, the increased RUNX2 expression was obviously related to worse FP and OS in LC patients throughout a 200-month follow-up period. Furthermore, RUNX2 mRNA and protein were found to have a suitable high expression level in LC tissue when analyzing the expression profile using the HPA database (Fig. 3). Nevertheless, RUNX2 alteration information was detected by cBioPortal, and joint analysis was conducted using GO enrichment. It was implicated that the cell proliferation, invasion and metastasis-related genes, as well as the functions of cell proliferation, adhesion, motion, migration which contributed to tumor formation, development and metastasis were significantly associated with RUNX2 alterations (Fig. 6).
Future directions
Large scale experiments and multicenter clinical trials are still needed to validate this supposition. Besides, further analysis of online data concerning RUNX2 in lung cancer, ruling out confounders such as age, tumor stage, and recurrence status should be further refined.
Conclusions
The current study recognized the increased expression level of RUNX2 in LC. Moreover, the upregulated RUNX2 expression level manifested as the poor prognosis of FP and OS. RUNX2 alterations were closely related to cell proliferation, invasion and metastasis-related genes as well as the functions of cell proliferation, adhesion, motion and migration. We provided a novel method that utilizes available public data to mine the relationship between oncogene and cancer. Our findings showed that the elevated RUNX2 expression level in LC patients and its potential biological functions lead to its possible value as a new predictive factor and effective drug target. Certainly, additional large-scale studies and multicenter clinical trials are needed to corroborate this finding.
Supporting information
Supplementary material for this article is available at https://doi.org/10.14218/ERHM.2021.00009 .
Table S1
The genes correlated with RUNX2 in lung adenocarcinoma (cBioPortal database).
(PDF)
Table S2
The genes correlated with RUNX2 in squamous cell lung carcinoma (cBioPortal database).
(PDF)
Abbreviations
- RUNX2:
Runt-related transcription factor 2
- EMT:
epithelial mesenchymal transition
- TGF-β:
transforming growth factor-β
- LC:
lung cancer
- LUAD:
lung squamous cell carcinoma
- LUAD:
lung adenocarcinoma
- NSCLC:
non-small cell lung cancer
- GEPIA:
gene expression profiling interactive analysis
- FP:
progression-free survival
- OS:
overall survival
- PPS:
post progression survival
- HPA:
human protein atlas
- PPI:
protein–protein interaction
- GO:
gene ontology
- BP:
biological process
- CC:
cellular component
- MF:
molecular function
Declarations
Data sharing statement
All data are available from the corresponding authors.
Funding
The study is funded by the Interdisciplinary research of 9th People's Hospital affiliated to Shanghai Jiao Tong University School of Medicine (YG2019QNA13).
Conflict of interest
The authors have no conflicts of interest related to this publication.
Authors’ contributions
Study design and administration (JH, XBZ); performance of experiments, analysis and interpretation of data, manuscript writing (DX, KL); critical revision (JH); technical and material support (JC, YYG). All the authors have made a significant contribution to this study and have approved the final manuscript.