Expression and Gene Regulation Network of GRK6 in Clear Cell Renal Cell Carcinoma based on Comprehensive Bioinformatics Analysis

Zexun Deng Dalian Medical University Xiaomei Li Dalian Medical University Xinxin Peng Dalian Medical University Weijie Zhao Dalian Medical University Tiancheng Xia Dalian Medical University Taoze Ji Dalian Medical University Junjie Yang Dalian Medical University Yinxia Wu Northern Jiangsu People's Hospital Junjie Yu (  urologistyjj@163.com ) Peking University First Hospital Liqun Zhou Peking University First Hospital


Background
Renal cell carcinoma (RCC), one of the 10 most common cancers, occupies approximately 3.9% of new carcinomas, with an increasing incidence in last twenty years [1]. Because of its to resistance to radiotherapy [2] and chemotherapy [3], surgery remains the major treatment for patients with RCC. The most common subtype of renal cell carcinoma is clear cell RCC, and ccRCC accounts for approximately 75-80% of RCC [4], which is de ned by its lipid-rich and glycogen-rich cytoplasm [5]. Although the recurrence rate and mortality rate has decreased with the use of robotic and laparoscopic surgery, nephron-sparing surgery and minimally invasive procedures, and biologic response modi ers are applied to patients with metastatic RCC, the prognosis of terminal cancer is very poor, with the 5-year survival rate of 5-10% [6,7]. Therefore, it is necessary to discern novel biomarkers and targets for the diagnosis, treatment and prognosis of RCC.
G protein-coupled receptor kinases (GRKs) are a large protein kinase family, which regulates the activity of G protein-coupled receptors via phosphorylating their intracellular domains [8,9]. It has been reported that GRKs participate in a wide variety of human diseases, including heart failure, homologous desensitization, opiate addiction and malignant tumor development [8,10,11]. G protein-coupled receptor kinase 6 (GRK6) is one of the GRK family member, which plays a important role in many diseases process [12][13][14].Since then, a new role of GRK6 in a number of cancer metastases was gradually revealed, such as lung cancer [15], medulloblastoma [16], and multiple myeloma [14]. Our previous studies have implied that GRK6 was overexpressed in ccRCC tissue and that GRK6 expression levels predict overall survival in ccRCC patients. However, it is indispensable to examine the speci c regulatory mechanisms and functions of GRK6 in ccRCC.
In our previous studies, we indicated that high GRK6 expression was found in ccRCC tumor tissues compared to normal kidney tissues, and it was associated with patient's age, pathological grading and so on. Moreover, Kaplan-Meier survival analysis showed that the overexpression of GRK6 was associated with progression-free survival and poor overall survival of ccRCC patients. These results imply that GRK6 is a novel protooncogene. Thus, we studied GRK6 expression and mutations in data from patients with ccRCC in The Cancer Genome Atlas (TCGA) and other public databases. We analyzed genomic alterations and functional networks related to GRK6 in ccRCC by using multi-dimensional analysis methods. Therefore, new targets and strategies may be revealed by our results for ccRCC diagnosis and treatment.

Oncomine analysis
The mRNA expression and DNA copy number of GRK6 in ccRCC were analyzed within the Oncomine 4.5 database. Oncomine (www.oncomine.org), currently the world's largest oncogene chip database and integrated data mining platform, contains 715 gene expression data sets and data from 86,733 cancer tissues and normal tissues [17]. This analysis drew on a series of ccRCC studies, including Yusenko Renal, Gumz Renal, Jones Renal, TCGA Renal 2 and Lenburg Renal studies [18][19][20][21]. GRK6 expression was assessed in ccRCC tissue relative to its expression in normal tissue, and differences associated with p 0.01 were considered signi cant. UALCAN analysis UALCAN uses TCGA level 3 RNA-seq and clinical data from 31 cancer types, which is an interactive webportal to perform in-depth analyses of TCGA gene expression data [22]. One of the portal's user-friendly features is that it allows analysis of relative expression of a query gene(s) across tumor and normal samples, as well as in various tumor sub-groups based on individual cancer stages, tumor grade or other clinicopathological features. UALCAN is publicly available at http://ualcan.path.uab.edu. c-BioPortal analysis Page 4/16 The cBio Cancer Genomics Portal (http://cbioportal.org) is an open-access resource for interactive exploration of multidimensional cancer genomics data sets, currently containing 225 cancer studies [23]. We used c-BioPortal to analyze GRK6 alterations in the TCGA KIRC sample. The search parameters included mutation, CNVs, and mRNA expression. The tab OncoPrint displays an overview of genetic alterations per sample in GRK6. The tab Network visualizes the biological interaction network of GRK6 derived from public pathway databases, with color-coding and lter options based on the frequency of genomic alterations in each gene. Neighboring genes with alteration frequencies greater than 10% were included. We then performed GO and KEGG pathway enrichment analyses of the most frequently altered neighbor genes using the "clusterPro ler" package in R [24]. The GO annotation had three parts: cellular component (CC), biological process (BP), and molecular function (MF).

LinkedOmics analysis
The LinkedOmics database (http://www.linkedomics.org/login.php) is a Web-based platform for analyzing 32 TCGA cancer-associated multi-dimensional datasets [25]. The LinkFinder module of LinkedOmics was used to study genes differentially expressed in correlation with GRK6 in the TCGA KIRC cohort (n=533). Results were analyzed statistically using Pearson's correlation coe cient. The LinkFinder also created statistical plots for individual genes. All results were graphically presented in volcano plots, heat maps or scatter plots. The LinkInterpreter module of LinkedOmics performs pathway and network analyses of differentially expressed genes. The comprehensive functional category database in the Webbased Gene SeT AnaLysis Toolkit (WebGestalt) [26] was applied. Data from the LinkFinder results were signed and ranked, and GSEA was used to perform analyses of GO (CC, BP and MF), KEGG pathways, kinase-target enrichment, miRNA-target enrichment and transcription factor-target enrichment. The latter two network analyses were based on the Molecular Signatures Database (MSigDB) [27]. The rank criterion was an FDR < 0.05, and 500 simulations were performed.

GeneMANIA analysis
GeneMANIA (http://www.genemania.org) is a exible, user-friendly web interface for constructing proteinprotein interaction (PPI) network, generating hypotheses about gene function, analyzing gene lists and prioritizing genes for functional assays [28]. The website can set the source of the edge of the network, and it features several bioinformatics methods: physical interaction, gene co-expression, gene colocation, gene enrichment analysis and website prediction. We used GeneMANIA to visualize the gene networks and predict function of genes that GSEA identi ed as being enriched in KIRC: kinase IKBKB, miRNA 192, miRNA 215 and transcription factor NFKB.

GRK6 expression in ccRCC
We have preliminary analyzed GRK6 transcription levels in lots of ccRCC studies from the Cancer Genome Atlas database (TCGA) and Gene Expression Omnibus (GEO). The data in the Oncomine 4.5 database showed that DNA copy number variation (CNV) and mRNA expression of GRK6 were apparently higher in ccRCC tissues than those in normal tissues (p 0.01). Although the fold differences were within 2.5, GRK6 ranked within the top 5% based on DNA CNVs and within the top 30% based on mRNA expressiontop ( Figure 1). Further sub-group analysis of various clinic pathological characteristics of 533 Kidney renal clear cell carcinoma (KIRC) samples in the TCGA consistently showed high GRK6 transcription. According to a subgroup analysis based on age, race, gender, disease stage, and tumor grade, the GRK6 transcription level in ccRCC patients was signi cantly higher than that in healthy people ( Figure 2). Therefore, GRK6 expression can be used as a potential marker for diagnosis in ccRCC.

Genomic alterations of GRK6 in ccRCC
Frequency and type of GRK6 alterations in ccRCC Then, based on the sequencing data from KIRC patients in the TCGA database, we used cBioPortal to determine the type and frequency of GRK6 changes in ccRCC. GRK6 was altered in 97 of 446 (22%) KIRC patients ( Figure 3A). These alterations were mRNA upregulation in 25 cases (5.6%), downregulation in 3 cases (0.7%), ampli cation in 69 cases (15.5%), mutation in 4 case (1.1%), and multiple alterations in 4 cases (1.2%) ( Table 1). Thus, ampli cation is the most common type of GRK6 CNV in ccRCC.

Biological interaction network of GRK6 alterations in ccRCC
Next, we completed the biological interaction network of GRK6 in ccRCC. To do this, we used the tab in cbioportal to display GRK6 neighboring genes, which changed at > 10% frequency ( Figure 3B and Table   1). Use networks to identify frequently changing neighbor genesThe neighbor genes of GRK6 with the most frequent alterations were CCR8 (26.9%), HRH2 (24.0%) and GRM6 (19.1%). Use Network to identify frequently altered neighbor genes. Analysis of signi cantly enriched gene ontology (GO) terms revealed that these genes encode proteins localized mainly to the plasma membrane, and integral component of plasma membrane. These proteins are primarily involved in G-protein coupled receptor signaling pathway, in ammatory response and chemokine-mediated signaling pathway, while they also serve as chemokine receptor activity and G-protein coupled receptor activity ( Figure 4A-4C). Similarly, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed enrichment in chemokine signaling pathway, cytokine-cytokine receptor interaction and neuroactive ligand-receptor interaction ( Figure 4D, 4E). Therefore, GRK6 altered biological interaction network was involved in chemokine signaling pathway. Enrichment analysis of GRK6 functional networks in ccRCC Analysis of GO and KEGG pathways of co-expressed genes related to GRK6 in ccRCC The mRNA sequencing data of 533 KIRC patients in TCGA can be analyzed by the functional module of LinkedOmics. As shown in the volcano plot ( Figure 5A), there are 7,971 genes (dark red dots) that have a strong positive correlation with GRK6, while 4,622 genes (dark green dots) have a strong negative correlation (false discovery rate [FDR] < 0.01). As shown in the heat map, these 50 important gene sets are positively and negatively correlated with GRK6 ( Figure 5B, 5C). This result indicates that GRK6 has a wide range of effects on the transcriptome. As shown in Supplementary Figure 1A GRK6 networks of kinase, miRNA or transcription factor targets in ccRCC In order to nd the targets of GRK6 in ccRCC, we analyzed the kinase, miRNA and transcription factor target networks of positively related gene sets generated by GSEA. The top 5 most signi cant target networks were the kinase-target networks related primarily to the kinases inhibitor of nuclear factor kappa B kinase subunit beta (IKBKB), p21 (RAC1) activated kinase 1 (PAK1), LCK proto-oncogene, Src family tyrosine kinase (LCK), conserved helix-loop-helix ubiquitous kinase (CHUK) and serum/glucocorticoid regulated kinase 1 (SGK1) ( Table 2 and Supplementary Tables 5-7

Discussion
Many different functions of G protein-coupled receptor kinase 6 have been reported in various cancers, such as tumor progression, angiogenesis, metastasis and cell proliferation etc. Chen et al suggested that GRK6 plays an important role in cell adhesion and migration of prostate cancer cell lines, breast cancer cell lines and cervical cancer cell lines [29]. In addition, GRK6 provides a signal transduction scaffold that regulates cell adhesion and cytoskeletal organization through ArfGAP 1 interaction with G protein coupled receptor kinase and indirect transactivation of epidermal growth factor receptors. This in turn affects the migration and invasion of cancer cells [30]. Due to secondary messengers (including cAMP and calmodulin), GRK6 can affect the migration and invasion of cancer cells [31]. We conducted a bioinformatics analysis of public sequencing data to gain a more thorough understanding of GRK6's functions in ccRCC and its regulatory network, and it will better guide the future research of ccRCC.
RCC lacks particular symptoms at early stage and effective noninvasive methods for screening, leading to the present situation that the treatment of early RCC is hampered. Therefore, new RCC markers are needed to improve early diagnosis. Analysis of transcript sequencing data from more than 1,000 clinical samples in the GEO and TCGA databases including ve different ccRCC study data sets [18][19][20][21] con rmed that GRK6 mRNA levels and CNV in ccRCC were signi cantly higher than normal kidney tissues. We nd that GRK6 overexpression occurs in multiple cases of ccRCC, which is worthy of further clinical veri cation and has the potential to be a valuable marker for diagnosis and prognosis.
CNVs can have signi cant genomic implications, disrupting genes and altering genetic content, resulting in phenotypic differences [32]. The increase in the copy number of GRK6 in ccRCC was found by our research, and that the major type of GRK6 alteration was ampli cation, which was associated with shorter survival period. We speculate that altered GRK6 expression and GRK6 dysfunction in ccRCC may be caused by alterations in chromosomal structure. Due to GRK6 plays many important physiological functions, its alteration may lead to changes in a wide variety of downstream signaling pathways. In fact, neighboring gene networks close to GRK6 usually show varying degrees of ampli cation in ccRCC. Related functional networks are involved in chemokine signaling pathway, cytokine-cytokine receptor interaction and neuroactive ligand-receptor interaction. Thus, the network of GRK6 alterations is involved in the core node of post-transcriptional regulation, which is closely related to regulating various enzymes stability and functions. Accumulating results show the nonnegligible roles of modi cation enzymes in tumor progression, such as the ubiquitinases [33], deubiquitinases [34], kinases [35], and phosphatases [36]. GRKs were rst reported to catalyze the phosphorylation and desensitization of G protein-coupled receptors (GPCRs). Nevertheless, recent studies identi ed its GPCR-independent roles in regulating cellular functions [37]. For example, GRK2, GRK5, and GRK6 can mediate TNFa-induced NF-kB signaling via direct phosphorylation of IkBa [38]. Moreover, GRKs were reported to be involved in tumor development and progression. GRK5 can functionally regulate well-known cancer-related proteins such as Wnt and tumor suppressor p53 [39]. Another example is the upregulation of GRK5 in glioblastoma stem cells, indicating its signi cance in promoting tumor proliferation [40]. The above research is consistent with the results of our functional network analysis.
The important network of target kinases, miRNAs and transcription factors can be revealed by the enrichment analysis of target gene sets by GSEA. Our results suggest that the functional network of GRK6 participates primarily in cytokine metabolic process, interferon-gamma production and adaptive immune response. These ndings are consistent with the fact that GRK6 can mediate TNFa-induced NF-kB signaling via direct phosphorylation of IkBa [38]. It is critical to understand how alteration in a protein important for ensuring normal transcription can lead to major dysfunction and even to cancer such as ccRCC.
In recent years, many studies have found that the tumor microenvironment is closely related to tumor invasion, growth, and metastasis. The interaction between chemokines and chemokine receptors recruits different immune cell subgroups into the tumor microenvironment, and these groups have an important impact on the occurrence and development of tumors. We found that GRK6 in ccRCC is associated with a network of kinases including IKBKB, PAK1 and LCK. These kinases regulate chemokine, cytokines and the adaptive immune response [41,42]. In fact, IKBKB can enhance expression of the cytokines activated, and alter most of the NF-κB activity discovered in solid tumors. As a fundamental regulator of NF-κB activity, not surprisingly, the activity of IKBKB is closely related to tumor development and progression, such as breast cancer [43], skin cancer [44], and gastric carcer [45]. In ccRCC, GRK6 may regulate chemokine, cytokines and tumor microenvironment via IKBKB kinase.
In the promoters and enhancers of many genes, there are functional NF-κB binding sites, and activated NF-κB can bind to speci c sequences on the DNA chain to initiate and regulate immune responses, mediate cell adhesion, differentiation, proliferation, apoptosis and in ammation, and play an important role in the occurrence and development of various immune in ammatory diseases and tumors [46,47]. When the body is stimulated by stimulating factors such as viruses, radiation, oxygen free radicals, bacterial lipopolysaccharides, TNFa etc, the NFkB-IkBs complex is activated, and NF-κB is dissociated and transferred to the nucleus, and speci cally binds to the corresponding site to promote the corresponding gene transcription, so that the biological reaction occurs. At the same time, most of the activators in NF-κB and its signaling pathways are considered to play an important role in the occurrence and development of tumors [48]. NF-κB activation induces DNA damage, oncogenic mutations and genomic instability leading to tumor initiation. Chronic in ammation and NF-κB can also cause chromosomal instability and aneuploidy. NF-κB enhances the proliferation of initiated tumor cells by promoting the production of various cytokines, growth factors and cell cycle proteins. Further studies should test this hypothesis.
Our study identi ed some miRNAs that were associated with GRK6. Actually these short noncoding RNAs, normally involved in post-transcriptional regulation of gene expression, can contribute to human carcinogenesis [49]. The particular miRNAs in our study have been linked to tumor proliferation, apoptosis, invasion, metastasis, cytokine metabolic. In fact, miR-190-5p is regulated by NFκB1/p50, thereby restoring the apoptotic response of cells.   The dysregulation of these miRNAs is closely related to the occurrence and development of tumors. We analyze the relationship between GRK6 and these miRNAs through bioinformatics to guide further research.
Our study strongly demonstrates the importance of GRK6 in carcinogenesis and its potential as a ccRCC marker, and the results show that GRK6 overexpression in ccRCC has a profound impact on genome stability, multiple steps of gene expression (DNA replication, RNA splicing and protein translation) and cytokine metabolism. GRK6 is related to several tumor-associated kinases (such as IKBKB), miRNAs (such as miRNA-190), and transcription factors (such as NF-κB). Our research is based on the most prevalent bioinformatics theory and uses online tools to analyze target genes from tumor data from public databases. Meanwhile, this method has the advantages of low cost, large sample size and simplicity compared with traditional chip screening [54]. This enables large-scale ccRCC genomics research and subsequent functional studies.
However, the TCGA database has its own limitations. One is that the TCGA KIRC samples contain three ethnic groups. The genetic background and etiology of KIRC can differ signi cantly across ethnic groups. Another limitation is that ccRCC patients are not easy to be found in the early stage of the disease. Under normal circumstances, the disease has progressed to the middle and late stages when obvious symptoms appear. But the KIRC samples contain relatively few patients in stage 4. Therefore, we should increase the number of patients of different races and KIRC stages to enrich our clinical samples and make the results more accurate. The third limitation is that transcriptome sequencing can detect only static mutations; it cannot directly provide information on protein activity or expression level. These questions should be addressed in follow-up experiments.

Conclusions
GRK6 is overexpressed and ampli cation is the most common type of GRK6 CNV in ccRcc. The functional networks of GRK6 is in connection with cytokine metabolic process, chemokine regulation,interferon-gamma production and adaptive immune response. Functional network analysis suggested that GRK6 regulates proteasome, ribosome, base excision repair, DNA replication, RNA splicing and protein translation via pathways involving multiple tumor-associated kinases, miRNAs and transcription factors. This study potent demonstrates the importance of GRK6 in carcinogenesis and its potential as a ccRCC marker.

Tables
Due to technical limitations, table1 and Table 2 are only available as a download in the Supplemental Files section.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.