Expression patterns and prognostic value of RUNX genes in kidney cancer

Kidney cancer is the third most common malignancy of the urinary system, of which, kidney renal clear cell carcinoma (KIRC) accounts for the vast majority. Runt-related transcription factors (RUNX) are involved in multiple cellular functions. However, the diverse expression patterns and prognostic values of RUNX genes in kidney cancer remained to be elucidated. In our study, we mined the DNA methylation, transcriptional and survival data of RUNX genes in patients with different kinds of kidney cancer through Oncomine, Gene Expression Profiling Interactive Analysis, UALCAN, Kaplan–Meier Plotter, cBioPortal and LinkedOmics. We found that RUNX1 and RUNX3 were upregulated in KIRC tissues compared with those in normal tissues. The survival analysis results indicated a high transcription level of RUNX1 was associated with poor overall survival (OS) in KIRC patients. Furthermore, KIRC tumor tissues had significantly lower levels of RUNX1 promoter methylation than that in paracancerous tissues, with decreased DNA methylation of RUNX1 notably associated with poor OS in KIRC. In conclusion, our results revealed that RUNX1 may be a potential therapeutic target for treating KIRC, and RUNX1 promoter methylation level shows promise as a novel diagnostic and prognostic biomarker, which laid a foundation for further study.

www.nature.com/scientificreports/ displayed various post-translational modifications such as acetylation 14 , methylation 15 , phosphorylation 16 and sumoylation 17 . All of these contribute to the functional complexity of RUNX genes. So far, a few studies have investigated the expression profiles of RUNX genes in RCC 18,19 . However, the relationship between RUNX genes and the onset, progression and prognosis of RCC remains controversial. Therefore, our study aimed to analyze the connection between the expression levels of RUNX genes and clinicopathological parameters of RCC especially KIRC patients by multi-dimensional analysis methods, and conducted a preliminary analysis of their regulation and potential functions.

Materials and methods
Pan-cancer analysis of expression levels of RUNX genes by Oncomine and gene expression profiling interactive analysis (GEPIA). Initially, pan-cancer analysis was performed to assess the transcriptional levels of RUNX genes in different types of cancer and corresponding normal tissues using Oncomine, including 460, 450 and 447 total unique analyses for RUNX1, RUNX2 and RUNX3 respectively. Oncomine (https:// www. oncom ine. org/) is a platform that provide solutions for individual researchers and multinational companies with robust peer-reviewed analysis methods as well as a powerful set of analysis functions that compute gene expression signatures, clusters and gene-set modules and that allow extracting biological insights from the data automatically. This database contains 715 datasets, 86,733 normal and tumor samples 20 . Then, we utilized GEPIA in this study to verify the relationship between the expression of RUNX genes and kidney cancer. GEPIA (http:// gepia. cancer-pku. cn/) is an interactive web server for analyzing the RNA sequencing expression data of 9736 tumors and 8587 normal samples from the TCGA and the GTEx projects by using a standard processing pipeline 21 . In prognostic analysis, we grouped the patients by the median as a cut-off value.

UALCAN, cBio cancer genomics portal (cBioPortal) and CPTAC analyses. Then we obtained
and downloaded the clinical and FPKM-standardized RNA-seq data and the clinical information of 531 KIRC patients from the TCGA database. The clinical data were preprocessed by removal of samples without survival status and patients with survival time less than 30 days were also excluded because they might die of non-cancerrelated diseases. Furthermore, we performed univariate and multivariate Cox regression analysis to demonstrate the correlations between OS and clinical variables.
After the pan-cancer analysis, we explored the relationship between promoter methylation status and its mRNA expression of RUNX1 gene, using UALCAN and cBioPortal database. UALCAN (http:// ualcan. path. uab. edu/ index. html) is an interactive web resource for analyzing cancer data and providing gene expression analysis by using TCGA datasets (https:// www. cancer. gov) 22 . The Beta value indicates level of DNA methylation ranging from 0 (unmethylated) to 1 (fully methylated). Different beta value cut-off has been considered to indicate hypermethylation (Beta value: 0.7-0.5) or hypo-methylation (Beta-value: 0.3-0.25) 23,24 . The cBioPortal (http:// www. cbiop ortal. org/) is an open access for interactive exploration of multiple cancer genomics data sets, currently covering 282 cancer researches 25 . The National Cancer Institute's Clinical Proteomic Tumor Analysis Consortium (CPTAC, https:// prote omics. cancer. gov/) is a public data portal of proteomic sequence and corresponding genomic sequence datasets. Using CPTAC, we mined the protein expression of RUNX genes.

LinkedOmics analysis and gene set enrichment analysis (GSEA). After confirming that RUNX1
was associated with prognosis of KIRC, we preliminarily analyzed its specific cytological function and molecular mechanism. LinkedOmics (http:// linke domics. org/ login. php) is a publicly available portal that includes multiomics data from all 32 TCGA cancer types 26 . LinkedOmics was used to investigate the transcription networks of RUNX1 in KIRC. Data from the LinkFinder results were signed and ranked, meanwhile GSEA was used to perform Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways 27 , transcription factor-target enrichment and kinase-target enrichment. The statistical analyses were based on the Molecular Signatures Database (MSigDB) 28 ; 1000 simulations were performed and False discovery rate (FDR < 0.05) was used to select significantly enriched gene sets.

Ethical approval.
No ethical approval was obtained because this study did not involve a clinical evaluation, did not involve laboratory animals and invasive procedures.

Results
Analysis of expression level of RUNX genes in pan-cancer. As shown in Fig. 1, pan-cancer analysis results showed that RUNX1 increased significantly in 58 datasets, especially those of leukemia, head and neck cancer, colorectal cancer, kidney cancer, and breast cancer. As for RUNX2, 24 datasets displayed increased expression, whereas 14 datasets showed the opposite results. Similarly, higher expression of RUNX3 was found in 19 datasets, while lower expression was found in 10 datasets. Overall, the obtainedresults indicated significantly elevated expression of RUNX1 (with fold changes of 6.58 to 6.51) and RUNX3 (with fold changes of 2.7 to 8.32) in kidney cancer ( Table 1).

Verification of expression of RUNX genes in kidney cancer.
To further verify differential expression of RUNX genes in kidney cancer, we compared expression profiles of each gene between different kinds of kidney cancer samples and paired normal tissues by GEPIA. The results indicated that RUNX1 and RUNX3 were significantly overexpressed in KIRC ( Fig. 2A, C). For RUNX2, no significant expression difference was found (Fig. 2B). Furthermore, the relationship between expressions of RUNX1, RUNX2and RUNX3 and prognosis of KIRC were performed by GEPIA and the results revealed that higher expression of RUNX1 had poorer overall www.nature.com/scientificreports/   Fig. 2D). Especially, higher expression of RUNX1 was also found correlated to poor prognosis in KIRC (p < 0.001, Fig. 2E) and KIRP patients (p = 0.024, Fig. 2G), whereas there was no significant relation between RUNX1 and OS in KICH patients (p = 0.27, Fig. 2F). Besides, RUNX2 expression showed significant correlation to the prognosis of KIRC patients (p = 0.046, Fig. 2H), while RUNX3 did not (p = 0.77, Fig. 2I).  www.nature.com/scientificreports/ Correlation analysis between RUNX1 and clinicopathological characteristics in KIRC. Table 2 showed the basic clinical characteristics. In total, we identified 531 KIRC patients with RUNX1 expression data and clinical information. The results showed that RUNX1 was highly expressed in female patients and white race. Higher RUNX1 expression was associated with advanced TNM stage and poor histological grade stage.
The results were consistent with that RUNX1 might be an unfavorable factor for KIRC patients. Furthermore, we performed Cox regression analysis, and results showed stage and age were significantly associated with OS in KIRC patients (Fig. 3A). Then Multivariate Cox regression analysis showed that stage was an independent factor influencing KIRC prognosis (Fig. 3B).

Analysis of promoter methylation status and protein expression of RUNX1 gene in KIRC.
To explore the hinge of RUNX1 expression, we investigated the promoter methylation level of RUNX1 in KIRC by UALCAN. Twelve probes in RUNX1 promoter were used for detecting DNA methylation level of RUNX1 (Fig. 4A). Notably, primary tumor tissues had obviously lower promoter methylation levels than normal tissues (p < 0.001, Fig. 4B). Meanwhile, it was a significant inverse correlation between DNA methylation level of RUNX1 gene and its mRNA expression in KIRC samples (Spearman = − 0.69, p = 1.33e−46; Pearson = − 0.60, p = 1.19e−32, Fig. 4C) based on cBioPortal analysis. The results indicated that upregulated expression of RUNX1 was associated with DNA hypomethylation, and it could be considered as a risk factor for KIRC. Surprisingly, the prognosis analysis showed that patients with promoter hypomethylation of RUNX1 had a worse OS (p < 0.001, Fig. 4D). Furthermore, we compared the protein expression of RUNX1 in normal and KIRC tissues, verifying that the expression level of RUNX1 protein was indeed significantly elevated in KIRC (p < 0.001, Fig. 4E).

Enrichment analysis of RUNX1 gene functional networks in KIRC.
As shown in the volcano plot using LinkedOmics analysis (Fig. 5A), red spots represented genes positively correlated with RUNX1 in KIRC, while green spots represented genes with a negative correlation. Furthermore, the top 50 significant gene correlated positively or negatively with RUNX1 were shown in the heat maps (Fig. 5B, C), which suggested a widespread impact of RUNX1 in the transcription. In addition, GSEA showed varying expression of RUNX1 gene mainly in the receptor complex, cell-substrate junction and apical part of cell, which primarily participated in adaptive immune response, leukocyte migration and protein targeting. They played a key role in receptor-ligand activity, protein tyrosine kinase activity, transferase activity, and transferring acyl groups ( Fig. 6A-C). KEGG pathway analysis showed that the differentially expressed genes were mainly enriched in proteoglycans in cancer, microRNAs in cancer, focal adhesion and peroxisome (Fig. 6D-F). To further explore the targets of RUNX1 gene in KIRC, we analyzed the transcription factors and kinase targets of positively correlated gene sets generated by GSEA ( Table 3). The top 3 most significant kinase target networks related to LYN proto-oncogene, p21 activated kinase 1(PAK1) and HCK proto-oncogene3. The transcription factor target network was mainly related with ETS1, NFKAPPAB, MZF1, IRF, SRF.

Discussions
The RUNX family has been noticed to play an important role in leukemia and solid tumors. Known as one of the most frequently mutated genes in human leukemia, RUNX1 is originally identified to have a role in hematopoiesis 29 . Increasingly, it has been implicated in cancers of ovary 30 18 . According to a recent research from Kamikudo et al., as the RUNX family has a mechanism to compensate for loss among the family members, it is difficult to individually inhibit RUNX family proteins, RUNX family cluster regulation might be a cancer treatment strategy 34 . However, As Lie-a-ling mentioned, although it has become evident that expression level of RUNX1 can be used as a marker of tumor progression 35 , it is not yet fully uncovered how the alteration contributes to tumorigenesis, since both the amount and activation status of proteins can have effects. Therefore, we performed comprehensive analyses of the expression levels of RUNX genes in kidney cancer. www.nature.com/scientificreports/ Our study showed that compared to normal kidney tissues, the expression levels of RUNX1 and RUNX3 were increased in KIRC cases. Prognostic studies further suggested that higher expression of RUNX1 gene was significantly associated with poorer OS in KIRC. Interestingly, the lower promoter methylation level of the RUNX1 gene was found to be significantly related to its higher mRNA and protein expression and, consequently, to the poorer OS of KIRC. In a research reported by Matsumura T et al., defection of RUNX1 methylation in hematopoietic stem cells (HSCs) was found to inhibit apoptosis and provide cells with a growth advantage, which is one   Recent studies have found overexpression of RUNX1 in mouse models of kidney fibrosis, which is related to KIRC, indicating RUNX1 has a regulation of TGFβ-driven epithelial-to-mesenchymal transition (EMT) 40 . Recently, Young et al. indicated that RUNX1, which had a relationship with multiple signaling pathways including JAK/STAT , MAPK, p53 and VEGF, could be recognized as a novel therapeutic target and prognostic factor 41 . Kamikudo et al. found that moderate inhibition of RUNX1 most significantly increased the total level of RUNX family through "genetic compensation of RUNX family transcription factors", emphasized the role of RUNX1 in tumorigenesis 42 . Furthermore, Zhao et al. manifested that PRMT1-dependent methylation of RUNX1 likely contributed to its inhibitory activity 15 . Considering our findings, patients with higher transcription levels of RUNX1 had worse prognosis, which corresponded to lower levels of promoter methylation. The results suggested that the promoter methylation of RUNX1 occurred in KIRC cases and deserved to be seen as a potential diagnostic and prognostic marker.
Previous studies suggested that RUNX1 was critical in a variety of genes transcription and played a role in cellular regulation. It has been recognized that genomic instability and mutagenesis are essential features of cancer cells, and kinases and their associated signaling pathways help stabilize and repair genomic DNA 43,44 . Ballissimo et al. noticed that cells with inefficient RUNX1 showed defects in DNA repair, including base excision, homologous recombination and DNA interstrand crosslink repair 45 . Sanoji et al. also proved that RUNX1 was involved in cell cycle arrest and apoptosis, as a potential factor in cancer formation 46 . In order to clarify the role of RUNX1 in KIRC, we tried to locate its target kinases and transcription factors by LinkedOmics, and found an extensive connection with them, indicating that RUNX1 was rather involved in cellular regulation. Furthermore, GSEA was performed to identify significantly enriched or depleted groups of genes. The results showed that RUNX1 was mainly responsible for adaptive immune response, leukocyte migration and protein targeting. In addition, RUNX1 mainly participated in peroxisome, microRNA and proteoglycans in cancers. We assume that these participations are likely to make RUNX1 an initiator of KIRC.
In conclusion, we performed an integrated analysis about the expression and prognostic value of RUNX genes in kidney cancer. Our results showed that RUNX1 and RUNX3 were upregulated in the tissues of KIRC compared to the normal ones. Furthermore, the results revealed that RUNX1 was a potential therapeutic target for KIRC, and that lower promoter methylation level of RUNX1 indicated poorer survival. Generally, we assumed RUNX1 may be a potential diagnostic and prognostic marker for KIRC, and its abnormal promoter methylation may participate in tumorigenesis, which lay a foundation for further study.