Glia Maturation Factor β as a Novel Independent Prognostic Biomarker and Potential Therapeutic Target of Kidney Renal Clear Cell Carcinoma

Abstract

Kidney renal clear cell carcinoma (KIRC) has the highest mortality rate and potential for invasion among renal cancers. The diagnosis and treatment of KIRC are becoming challenging because of its diverse pathogenic mechanisms. Glia (GMFB) is a highly conserved growth and differentiation factor for glia cells and neurons, and it is closely associated with neurodegenerative diseases. However, its role in KIRC remains unknown. The present study integrated bioinformatics approaches with suitable meta-analyses to determine the position of GMFB in KIRC. There was a significant decrease in Gmfb expression in KIRC kidneys compared with normal controls. Gmfb expression was negatively associated with pathologic stage, T and M stages, and histologic grade. Univariate and multivariate analyses showed that elevated Gmfb expression was an independent factor for a favorable prognosis. Furthermore, the nomogram verified that Gmfb is a low-risk factor for KIRC. Knockdown of Gmfb in Caki-2 cells increased viability and decreased p21 and p27 levels. Overexpression of Gmfb inhibited Caki-2 cell proliferation, migration, and invasion and decreased mitochondrial membrane potential. Gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses considering Gmfb co-expressed differentially expressed genes (DEGs) showed that collecting duct acid secretion and mineral absorption ranked were the most important upregulated and downregulated DEGs, respectively. The upregulated hub genes for DEGs were mainly involved in nucleosome assembly, nucleosome organization, and chromatin assembly, and the downregulated hub genes were primarily associated with keratinization. The ratio of tumor-infiltrating immune cells in KIRC tissues was evaluated using CIBERSORTx. The results showed that the Gmfb expression was significantly positively correlated with macrophage M2 cells and mast resting cell infiltration levels and negatively correlated with T follicular helper, T regulatory, and B plasma cell infiltration levels. The former cell types were associated with a beneficial outcome, while the latter had a worse outcome in patients with KIRC. In summary, this study identified GMFB as a novel independent biomarker and therapeutic target for KIRC, and it provides a helpful and distinct individualized treatment strategy for KIRC with a combination of molecular targets and tumor microenvironment.

Introduction

Kidney cancer is one of the ten most common types of cancer worldwide (1, 2). It is estimated that 300,000 people worldwide develop kidney cancer each year, and approximately half of them die from this disease (3). Kidney renal clear cell carcinoma (KIRC), which originates from the epithelial cells of the proximal convoluted tubule of the nephron, has become one of the most rapidly evolving areas of oncology, accounting for approximately 75% of kidney cancers (3–5). In addition to being the most common subtype, KIRC has the highest aggressiveness and metastasis and the highest mortality rate (4), as well as the highest immunity and vascular infiltration of cancer types (6, 7).

The molecular mechanisms involved in KIRC pathogenesis are diverse (8), including direct chromatin configuration alteration, vasculogenesis, and glucose metabolism (8, 9). Over the past few years, many studies have focused on the oncogenes of renal cell carcinoma, resulting in the discovery of a repetitive mutation in KIRC (10). However, oncogenes and suppressors such as RAS genes, TP53, Rb, CDKN2A, BRAF, PIK3CA, PTEN, EGFR, and ERBB2, which are highly mutated in other adult epithelial cancers, play only a minor role in KIRC (10). The diagnosis and treatment of KIRC are becoming a challenge for the oncology community. Currently, the most reliable prognostic predictors for KIRC are T, N, and M stages (11). However, in these indicators, survival rates may vary considerably, even among patients with the same set of illnesses (11, 12). Therefore, there is an urgent need to identify novel, reliable prognostic biomarkers and therapeutic targets for patients with KIRC.

Recently, bioinformatics techniques have become crucial for deciphering the molecular mechanisms underlying KIRC. A considerable number of KIRC-related high-throughput microarray datasets and sequencing datasets have been accumulated, resulting in groundbreaking changes in KIRC research. Multi-omics have identified some genes, such as CHAC1 (13), METTL14 (14), P4HB (15), and FOXM1 (16), as novel prognostic indicators for KIRC. However, the analysis of the tumor microenvironment (TME) linked to the specific targets of KIRC is lacking.

It is well known that heterogeneity, an essential feature of tumors, may contribute to the distinct responses of cancer patients to the same treatment. Single-cell sequencing technology allows for dissecting tumor heterogeneity at the molecular and cellular levels. The TME is increasingly considered a target for drug treatment strategies (17). Recently, Hu et al. revealed intratumoral heterogeneity in KIRC, resulting in different clinical outcomes in single-cell transcriptomes (18).

Glia maturation factor-beta (GMFB), a highly conserved brain-specific protein, is a 17-kDa brain protein recognized as a growth and differentiation factor acting on neurons and glia (19, 20). It is essential for brain growth and development, and it is widely expressed in mammalian neuronal, glia, and parenchymal cells (21). GMFB is phosphorylated by protein kinase A (PKA) (19). Furthermore, PKA-phosphorylated GMFB inhibits extracellular signal-regulated kinase (ERK) and enhances p38 (22). It has been shown that the upregulation of GMFB is closely related to neurodegenerative disorders such as Parkinson’s disease (23). Our previous work revealed that GMFB functions as a novel independent prognostic biomarker and therapeutic target for liver hepatocellular carcinoma (LIHC) (24). However, the role of GMFB in KIRC remains unclear.

The present study employed The Cancer Genome Atlas (TCGA) database and the CIBERSORTx tool (https://cibersortx.stanford.edu) to explore the potential effects of GMFB and GMFB-related TME, and the significant findings of molecular signatures were confirmed by quantitative real-time PCR (qRT-PCR), Western blotting, and immunostaining in human renal clear cell carcinoma cell line Caki-2. The obtained results allow further understanding of the interaction between GMFB and the TME in KIRC and thus provide a novel strategy for GMFB-based KIRC immunotherapies.

Methods

Data Collection and Differential Expression Analysis

RNA sequencing (RNA-seq), somatic mutations, and related clinical data were downloaded from TCGA (https://portal.gdc.cancer.gov/), containing 11,069 samples from 33 types of cancer (25). RNA sequence data from 539 patients with KIRC and 72 normal tissues were also downloaded from TCGA database. RNA-seq data and patient clinical information (Workflow Type: HTSeq-FPKM) were acquired using UCSC Xena (https://xena.ucsc.edu/), an online tool for the exploration of gene expression and clinical and phenotypic data. UALCAN (http://ualcan.path.uab.edu) (26) and the Human Protein Atlas (HPA) (https://www.Proteinatlas.org/) (27) were used to analyze GMFB protein expression patterns and immunohistochemical results in normal kidney tissue and KIRC samples.

Survival Analysis and the Relationship Between Gmfb Expression and Prognosis

Survival and clinical phenotype data were extracted for each sample downloaded from TCGA database. Overall survival (OS) and disease-free survival (DFS) were evaluated to determine the relationship between Gmfb expression and patient prognosis. The Kaplan–Meier curve and Cox regression analyses were performed to assess the association between Gmfb expression and OS in KIRC patients. Subsequently, multivariate analysis was employed to estimate whether Gmfb expression was an independent prognostic factor for survival in patients with KIRC.

Gene expression profiling interactive analysis (GEPIA) (http://gepia2.cancer-pku.cn/#index) was used to confirm the correlation between Gmfb expression and clinicopathological information in multiple cancer types (28).

Enrichment Analysis

The top 2,000 upregulated and downregulated differentially expressed genes (DEGs) between the high and low Gmfb expression groups, defined by Spearman’s rho above 0.30, were used for enrichment analyses. Pathways within the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) terms were assessed for enrichment using the KEGG Orthology-Based Annotation System (KOBAS) platform (http://kobas.cbi.pku.edu.cn/) (29–32) and Metascape (https://metascape.org/gp/index.html#/main/step1) (33, 34).

Assessment of Immune Infiltration

CIBERSORTx (http://cibersort.stanford.edu/) is a deconvolution algorithm that predicts the cell type components of intricate tissues based on normalized gene expression profiles; its results are usually in agreement with true estimates for many types of cancer (35). As so, CIBERSORTx was used to evaluate the relative proportions of 22 infiltrating immune cell subtypes in KIRC and LIHC (35). The percentages of the 22 immune cell subtypes derived from CIBERSORTx were regarded as a correction at p < 0.05 (36). The relationship between the expression of 22 subtypes of tumor-infiltrating immune cells (TIICs) and the survival of KIRC and LIHC patients was conducted using GEPIA2021 (http://gepia2021.cancer-pku.cns/) (37). The results of the survival analysis were considered statistically significant at p < 0.05.

Protein–Protein Interaction Network Construction and Hub Gene Screening

Protein–protein interaction (PPI) networks were constructed using the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) database (38) (version 11.5) and Metascape (http://metascape.org). Hub genes were identified using the CytoHubba plugin (39) in the Cytoscape software (40) (version 3.7.7). The maximal clique centrality (MCC) algorithm was employed to identify hub genes.

Cell Lines

The human embryonic kidney cell line 293T purchased from American Type Culture Collection (ATCC) and Caki-2 purchased from the Fenghui Institute of Biotechnology were used in this study. The HEK293T and Caki-2 cells were separately maintained in Dulbecco’s modified Eagle’s medium (DMEM)/high glucose and DMEM/nutrient mixture F12, respectively, supplemented with 10% fetal bovine serum (FBS) (Cat. No. 10091148, Gibco, Waltham, MA, USA) and 1% penicillin/streptomycin (Cat. No. 60162ES76, Yeasen, Shanghai, China) at 37°C and 5% CO₂.

Cell Transfection

The three pairs of small interfering RNA (siRNA) sequences of human Gmfb are shown in Supplementary Table 1A. A scrambled sequence was used as a siRNA control. The pcDNA-GMFB plasmid and pcDNA3.1 were kindly provided by Wan Sun. Caki-2 cells were transfected using Lipofectamine 2000 (Cat. No. T101, Vazyme, Nanjing, China), according to the manufacturer’s instructions. After transfection for 12, 24, and 48 h, cells were collected for subsequent experiments.

Western Blotting

Caki-2 cells were lysed in radioimmunoprecipitation assay buffer (Cat. No. P0013B, Beyotime, Jiangsu, China) with protease inhibitors (Thermo Fisher Scientific, Waltham, MA, USA) to extract whole-cell proteins. Protein concentration was quantified using a bicinchoninic acid (BCA) kit (Cat. No. 20201ES76, Yeasen, Shanghai, China). A sample (30 μg) of each protein was then separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis, transferred to a polyvinylidene fluoride membrane (Millipore, Burlington, MA, USA), and blocked with 5% bovine serum albumin (BSA) at room temperature for 1 h. After incubation overnight at 4°C with primary antibodies against GMFB (1:1,000; Cat. No. 10690-1-AP, ProteinTech, Rosemont, IL, USA), GMFG (1:1,000; Cat. No. D121824, Sangon Biotech, Shanghai, China), p21 (1:1,000; Cat. No. 10355-1-AP, ProteinTech, Rosemont, IL, USA), p27 (1:1,000; Cat. No. #3686, CST, Danvers, MA, USA), and actin (1:10,000; Cat. No. 20536-1-AP, ProteinTech, Rosemont, IL, USA), membranes were incubated with a secondary antibody for 2 h. The secondary antibody used was horseradish peroxidase (HRP)-conjugated goat anti-rabbit IgG (1:10,000; Cat. No. SA00001-2, ProteinTech, Rosemont, IL, USA). The proteins were visualized using enhanced chemiluminescence (ECL) (Cat. No. E412-02, Vazyme, Jiangsu, China). Tanon 5200S (Tanon Science & Technology, Shanghai, China) was used to evaluate the density of protein bands, and relative protein levels were quantified using ImageJ software (https://imagej.nih.gov/ij/).

Quantitative Real-Time PCR

Total RNA was extracted with TRIzol reagent (Cat. No. R401-01, Vazyme) and reverse transcribed into cDNA using the Revert Aid Reverse Transcription System (Takara, Dalian, China) according to the manufacturer’s protocol. PrimeScript RT polymerase (Cat. No. RR036A, Takara, Dalian, China) and SYBR Green Master Mix (Cat. No. FP205-02, Tiangen Biotech, Beijing, China) were used for conducting qRT-PCR on a LightCycler 96 Detection System (Bio-Rad Laboratories, Hercules, CA, USA). Data were analyzed using the 2^−ΔΔCt method (24). Oligos were synthesized by Sangon Biotech. The oligo sequences are shown in Supplementary Table 1B.

Immunofluorescence

Immunofluorescence analysis was performed to detect the colocalization of GMFB and Arp2/3 complex. 293T and Caki-2 cells were fixed with 4% paraformaldehyde for 15 min, permeabilized with 0.3% Triton X-100 in phosphate-buffered saline (PBS) for 30 min, blocked with 3% BSA for 1 h, and incubated with primary antibodies against GMFB (1:200 Cat. No. 10690-1-AP, 60062-1-Ig, ProteinTech, Rosemont, IL, USA), Arp2 (1:50; Cat. No. D260865, Sangon Biotech, Shanghai, China), and Arp3 (1:50; Cat. No. sc-48344, Santa Cruz, Dallas, TX, USA) overnight with immunofluorescence buffer at 4°C. Cells were then incubated with secondary antibodies (Invitrogen, Waltham, MA, USA), goat anti-mouse IgG(H+L) cross-adsorbed secondary antibody, Alexa Fluor 488 (Cat. No. A-11001); goat anti-mouse IgG(H+L) cross-adsorbed secondary antibody, Alexa Fluor 555 (Cat. No. A-21422); goat anti-rabbit IgG(H+L) cross-adsorbed secondary antibody, Alexa Fluor 488 (Cat. No. A-11008); and goat anti-rabbit IgG(H+L) cross-adsorbed secondary antibody, Alexa Fluor 555 (Cat. No. A-21428) for 1 h at room temperature. After three washes for 10 min each with 1× PBS, cell sections were counterstained and mounted with DAPI Fluoromount-G (Cat. No. D9542, Sigma-Aldrich, St. Louis, MO, USA) and examined under a confocal microscope (Leica, Wetzlar, Germany).

Determination of Cell Viability

Caki-2 cells were transfected with siRNA and pcDNA-GMFB using Lipo2000 (Cat. No. T101-01, Vazyme, Jiangsu, China). At 12, 24, and 48 h post-transfection, cell viability was tested using a Cell Counting Kit-8 (Cat. No. C0037, Beyotime, Jiangsu, China), according to the manufacturer’s instructions. This experiment was repeated four times.

Measurement of Mitochondrial Membrane Potential

Caki-2 cells were transfected with an empty vector or Gmfb overexpression (OE) plasmid. JC-1 staining (Cat. No. C2003S, Beyotime, Jiangsu, China) was performed to determine mitochondrial membrane potential (MMP), according to the manufacturer’s instructions. Cell visualization was performed using fluorescence microscopy (Leica, Wetzlar, Germany).

Wound Healing Assay

Caki-2 cells were seeded in six-well plates and grown to 80%–90% confluence. Confluent cell monolayers were scratched using a sterile 200-µl pipette tip and rinsed with PBS to remove scratched cells according to our previous work (24). Wound closure was observed after 24 h. The wound closure rate was calculated as follows: wound closure (%) = (the area of initial wound − the area of the final wound)/area of initial wound × 100.

Transwell Migration and Invasion Assays

The migration and invasion abilities of cells were assessed using transwell migration (without Matrigel) and invasion (with Matrigel) assays in 24 wells according to our previous work (24). Cell suspensions in an FBS-free medium were added to the upper transwell chambers (24-well, 8-μm pores; BD, Franklin Lakes, NJ, USA), whereas 600 µl of the medium was supplemented with 10% FBS and added to the lower transwell chambers. The cell count after positive staining with crystal violet indicated migration or invasion.

Results

Clinical Characteristics of Kidney Renal Clear Cell Carcinoma Patients

Detailed information on the clinical and gene expression data of the downloaded 539 samples is presented in Table 1. The median age of participants was 61 years. Among the 539 patients, 353 were male (65.5%) and 186 were female (34.5%). The 269 samples of participants under the age of 60 accounted for 49.9% of all samples. The topographic distribution included 51.6% T1 (n = 278), 13.2% T2 (n = 71), 33.2% T3 (n = 179), and 2% T4 (n = 11) cases. Only 16 patients (6.2%) showed lymph node involvement, while 78 cases (15.4%) had distant metastases. Most patients (467, 87.8%) were white.

Survival Atlas of Gmfb in Multiple Cancer Types and Differential Expression of GMFB Between Tumor and Normal Tissue Samples

The contribution of Gmfb to pan-cancer OS and DFS in the GEPIA2 database was compared using the Mantel–Cox test. The criterion for significant differences was a false discovery rate (FDR) <0.05, also known as p-value adjustment. OS analysis showed that Gmfb expression was negatively correlated with the hazard ratio (HR) of KIRC but positively correlated with the HR of LIHC in pan-cancers (Figure 1A). This was consistent with a previous report (24), in which Gmfb expression was negatively correlated with the HR of KIRC but positively correlated with the HR of adrenocortical carcinoma (ACC) (Figure 1A). Based on the correlation of Gmfb with OS and DFS in pan-cancer analysis, GMFB may be a promising novel prognostic factor for KIRC.

Figure 1

Gmfb expression levels were analyzed in pan-cancer tissues using TCGA dataset. Analysis of unpaired samples of 63 types of cancer tissues showed that Gmfb expression differed in 15 types of cancer tissues (Figure 1B), and it was downregulated in KIRC (p < 0.001) (Figure 1B). In the paired KIRC samples, the expression of Gmfb was also downregulated (p < 0.001) (Figure 1C). Protein expression analysis indicated that GMFB was weakly expressed in KIRC (p < 0.001) (Figure 1E). To verify GMFB expression at the protein level, we used the immunohistochemical results of KIRC provided by the HPA database (Figure 1F), which were consistent with those obtained from TCGA database. Gmfb was generally downregulated in renal cancer tissues.

In vitro, we confirmed that the Gmfb mRNA (Figure 1D) and protein (Figure 1G) levels were significantly decreased in Caki-2 cells when compared with HEK293T cells, which was consistent with the results of immunostaining of GMFB in HEK293T and Caki-2 cells (Figures 1H–I). We also detected the colocalization of GMFB and Arp2/3, which were known interaction partners of GMFB (Figures 1H,I), supporting that GMFB was an intracellular protein.

Correlation of Gmfb Expression With Clinical Characteristics of Kidney Renal Clear Cell Carcinoma

The association between Gmfb expression and clinical features of KIRC patients is summarized in Table 2. T stage (p < 0.001), M stage (p < 0.004), primary therapy outcome (p = 0.016), pathologic stage (p < 0.001), gender (p = 0.002), histological grade (p < 0.001), and OS events (p < 0.001) were all significantly associated with Gmfb expression.

The differential expression of Gmfb was examined according to age, gender, pathologic stage, T classification, N classification, M classification, and histological grade (G) of KIRC patients (Figure 2). Expression of Gmfb was significantly downregulated in men compared to women (p = 0.004) (Figure 2A); it was also significantly associated with pathologic stage (stage IV compared to stage I, p = 0.002; stage III compared to stage I, p < 0.001) (Figure 2B), T stage (T3 compared to T1, p < 0.001) (Figure 2C), M stage (M1 compared to M0, p = 0.007) (Figure 2D), and histological grade (G4 compared to G2, p < 0.001; G3 compared to G2, p = 0.003) (Figure 2E). There was no significant relationship between Gmfb expression and age (Figure 2F) or N classification (Figure 2G).

Figure 2

Prognostic Value of Gmfb Across Kidney Renal Clear Cell Carcinoma

OS analysis showed that Gmfb was a low-risk gene in kidney tissue (p = 0.00018) (Figure 3A). The Kaplan–Meier DFS analysis also demonstrated that patients with high Gmfb expression levels survived longer than patients with low Gmfb expression (p = 0.00013) (Figure 3B). Cox OS analysis demonstrated that a significant prognostic difference existed between the high and low Gmfb expression groups in both univariate (HR 0.541, 95% CI, 0.397–0.736, p < 0.001) and multivariate (HR 0.578, 95% CI, 0.374–0.892; p = 0.013) analyses (Table 3). A receiver operating characteristic (ROC) curve was used to evaluate the prediction accuracy of intensity. The Gmfb score, in terms of the ROC curve, was 0.905 (Figure 3C), indicating that Gmfb was a very sensitive biomarker for KIRC.

Figure 3

The Sankey diagram further demonstrated the association between age, gender, Gmfb expression, and clinical status (Figure 3D). Each column represents a characteristic (gender; T, N, and M stages; and expression of Gmfb; different colors represent the various states of these characteristics, and lines represent the distribution of the same sample in other characteristic variables). As shown in the diagram, most female and male patients progressed to stage I. Most patients with stage I disease had high Gmfb expression and a high survival rate. Conversely, patients with more severe diseases tended to have lower Gmfb expression and lower survival rates. A nomogram of OS was constructed to integrate Gmfb expression and some prognostic factors, including T, N, and M classifications and age (Figure 3E). A higher nomogram score represented worse predictive probability. The calibration curve assessed the performance of the nomogram for Gmfb, and the C-index for OS was 0.755 (Figure 3F). In conclusion, this nomogram is a good model for predicting the survival of patients with KIRC.

The Effects of Overexpression of Gmfb in Caki-2 Cells on Viability, Proliferation, Invasion, and Mitochondrial Membrane Potential

To clarify the association between Gmfb expression and cell viability in KIRC, the effects of Gmfb OE and knockdown (KD) on the viability of Caki-2 cells were examined. First, we designed three siRNA sequences for Gmfb, and KD efficiency was determined in HEK293T by Western blotting. As shown in Figure 4A, the siRNA1 with the highest KD efficiency in HEK 293T cells was used for the subsequent experiment (Figure 4A). Then, Gmfb KD and OE in Caki-2 cells were confirmed by Western blotting (Figures 4B, C) and qRT-PCR (Figures 4D, E). In parallel, the expression of GMFG protein, a homolog of GMFB, remained unchanged after transfection with siGmfb or Gmfb OE plasmid, indicating that Gmfb KD specificity was reasonable (Figures 4F, G). After transfection with siRNA and OE plasmids at 12, 24, and 48 h, a significant increase in viability after Gmfb KD was obtained at 24 h (p < 0.05) and 48 h (p < 0.01) (Figure 4H), while a significant decrease in viability after Gmfb OE was observed at 24 h (p < 0.01) and 48 h (p < 0.001) (Figure 4I). These results indicated that GMFB expression level was negatively correlated with tumor cell viability. Moreover, we investigate the effect of altered Gmfb expression on the Caki-2 cell cycle. At 48 h post-transfection of siGmfb in Caki-2 cells, the expression levels of p21 (CDKN1A) and p27 (CDKN1B) significantly decreased when compared with those in the scramble group (Figures 4J, K). The downregulation of p21 and p27 was required for the cellular transition from quiescence to the proliferative state (41). Gmfb OE in Caki-2 cells led to a lower number of migrated cells than in the empty vector group (Figures 4L, M). Invasion assays showed that Gmfb OE inhibited Caki-2 cell invasion (Figures 4L, N) and significantly reduced wound closure (Figures 4O, P). Detection of MMP using the JC-1 assay revealed a shift from red to green fluorescence in Gmfb OE Caki-2 cells (Figure 4Q), indicating that a high level of GMFB led to MMP loss. Taken together, we proposed that Gmfb was a suppressor gene in KIRC.

Figure 4

Functional Enrichment Analysis of Gmfb in Kidney Renal Clear Cell Carcinoma Patients

To elucidate the molecular mechanism of GMFB in KIRC, Gmfb co-expressed DEGs were used for functional enrichment analysis. The 3493 DEGs identified based on |log2 fold-change| > 1 (adjusted p < 0.05) corresponded to 912 upregulated genes and 2,581 downregulated genes. Detailed information on the DEGs is presented in Supplementary Table 2. DEG expression is shown in the volcano plot (Figure 5A), and the heatmap exhibits a corresponding hierarchical clustering analysis (Figure 5B). The details of the top 10 DEGs with upregulation and downregulation are listed in Table 4 and Supplementary Table 3, respectively. Most of these DEGs have not been reported in previous KIRC studies, except for ALDH6A1 and NCR3LG1. Six DEGs with unknown biological functions and three non-coding RNAs were not identified, including AC003072.1, AC073896.4, AC108673.3, AP001505.1, AC012510.1, AC084036.1, PWAR5, SCARNA10, and SNORA53. Among these, the mRNA levels of seven upregulated and four downregulated DEGs with known functions were further validated by qRT-PCR. Due to the very high cycle threshold values for ZNF196 and SELENOP (38 to 39), we did not redesign the primer for the validation experiment. Consistent with the results of the meta-analysis, INAFM1, HCFC1R1, and METTL26 were downregulated, and ATP6V0D2, TBC1D14, OXCT1, ALDH6A1, FAM160A1, and NCR3LG1 were upregulated in Gmfb OE Caki-2 cells (Figures 5C, D). The details of the nine validated DEGs are listed in Table 4, and the remaining 11 DEGs that were not validated are provided in Supplementary Table 3.

Figure 5

The KEGG pathway and GO enrichment analyses using KOBAS are shown in Figure 6. The neuroactive ligand–receptor interaction ranked 1 in the top 10 DEGs enriched KEGG pathways, followed by collecting duct acid secretion, linoleic acid metabolism, alpha-linolenic acid metabolism, and mineral absorption (Supplementary Figure 1A). GO enrichment analysis showed that cornification, extracellular region, and serine-type endopeptidase activity were the most enriched terms in biological process (BP), cellular compartment (CC), and molecular function (MF) categories, respectively (Supplementary Figure 1B and Supplementary Table 4). The top ten enriched KEGG pathways, GO terms, and the heatmap of associated DEGs are shown in Figures 6A–D (upregulated DEGs) and Figures 6E–H (downregulated DEGs). Upregulated DEGs were most enriched in collecting duct acid secretion (KEGG pathway), cellular protein metabolic process (BP category of GO), an integral component of the plasma membrane (CC of GO), and serine-type endopeptidase activity (MF of GO) (Figures 6A–D). Downregulated DEGs were most enriched in mineral absorption (KEGG pathway), cornification (BP of GO), extracellular region (CC of GO), and antigen-binding (MF of GO) (Figures 6E–H). There was more overlap in the KEGG pathway and GO analyses between total DEGs and downregulated DEGs than between upregulated DEGs, indicating that the downregulated DEGs may play a more critical and dominant role in KIRC. To this end, we selected five genes involved in the mineral absorption pathway (MT2A, 1G, 1F, TF, and S100G) for qRT-PCR validation of the RNA-seq results. The qRT-PCR revealed that MT2A, MT1G, and MT1F were downregulated, whereas S100G and TF did not show a significant change in Gmfb OE Caki-2 cells (Figures 6I, J).

Figure 6

To clarify the links between terms, PPI network and MCODE component analyses of the top 2,000 DEGs (Supplementary Figure 1C), upregulated (Figure 7A) and downregulated (Figure 7B), were performed using Metascape. Supplementary Table 5 provides detailed data. The PPI network showed that the predominant GO annotations of MCODE1 were cornified envelopes and keratinization in the total DEGs. Among the DEGs associated with upregulated Gmfb, GO annotations of MCODE1 were mainly linked with nucleosome assembly, nucleosome organization, and chromatin assembly. Among the DEGs related to downregulated Gmfb, GO annotations of MCODE1 were specifically associated with G-protein-coupled receptor (GPCR) downstream signaling, GPCR signaling, and GPCR ligand binding (Supplementary Table 5).

Figure 7

Identification of Hub Genes Using STRING and Cytoscape

Upregulated DEGs were projected separately into the human PPI network in the STRING database, and the hub genes were filtered using the MCC algorithm (Figure 8). The details of the top ten hub genes are shown in Supplementary Table 6. The top six upregulated hub genes (HIST1H2BB, HIST1H1A, HIST1H4F, HIST1H3J, HIST2H2AB, and HIST1H1B) (Figure 8B) belong to the histone cluster family and are mainly involved in nucleosome assembly, nucleosome organization, and chromatin assembly. The top ten downregulated hub genes included SPRR1B, SPRR1A, IVL, SPRR3, SPRR2A, TGM1, SPRR2D, SPRR2E, PI3, and CASP14 (Figure 8D). Except for CASP14 and TGM1, these genes were mainly associated with keratinization, and none has been reported in KIRC. These results suggest that GMFB may play a vital role in modulating nucleosome structure, cellular senescence, mitotic prophase, and keratinization in KIRC.

Figure 8

Correlation Between Gmfb Expression, Immune Cell Infiltration, and Prognostic Value of Tumor-Infiltrating Immune Cells in Kidney Renal Clear Cell Carcinoma

Based on the importance of TME, the relationship between immunity and Gmfb expression was further explored. Tumor RNA-seq data (retrieved from TCGA) of 33 different tumor patients and matched normal tissue samples from the Genomic Data Commons (GDC) website were downloaded to analyze the relationship between Gmfb expression and 22 types of TIICs. The heatmap of KIRC and LIHC analysis indicated Gmfb mRNA expression significantly correlated with TIICs (Figure 9A). The analysis by CIBERSORTx showed that the Gmfb mRNA expression positively correlated with the infiltration of T CD4+ memory resting cells (p = 0.00007, R = 0.172), monocytes (p = 5.68E−06, R = 0.196), mast resting cells (p = 0.047, R = 0.086), mast activated cells (p = 0.0017, R = 0.136), macrophage M2 cells (p = 1.81E−09, R = 0.257), macrophage M1 T cells (p = 0.0028, R = 0.13), eosinophil (p = 7.74E−11, R = 0.278), and B naïve cells (p = 9.95E−08, R = 0.229) but negatively correlated with the infiltration of regulatory T cells (Tregs) (p = 1.09E−31, R = −0.479), T follicular helper cells (TFH) (p = 1.43E−15, R = −0.337), T CD8+ cells (p = 4.97E−14, R = −0.319), natural killer activated cells (p = 7.42E−13, R = −0.305), B plasma cells (p = 0.0058, R = −0.12), and B memory cells (p = 0.00001, R = −0.185) (Figure 9A). However, our previous work showed that higher expression levels of Gmfb predicted poor prognostic outcomes in LIHC (26). To explore the possible mechanism underlying the differential effects of Gmfb on LIHC and KIRC, the correlation between Gmfb expression and TIIC in LIHC was analyzed. As shown in Supplementary Figure 2, the infiltration of B memory cells, mast activated cells, and neutrophils was significantly increased in LIHC (Supplementary Figures 2A–B) and was uniformly correlated with worse OS (Supplementary Figures 2C–E). None of these three immune cell types affected the OS of KIRC patients (Supplementary Figure 3). Moreover, the levels of B memory cells, mast activated cells, and neutrophils did not correlate significantly with Gmfb expression in the LIHC group (Figure 9A).

Figure 9

The association between TIICs and KIRC tumor tissues was examined. The relationship between the 22 TIICs and KIRC prognosis was investigated using GEPIA2021. The Kaplan–Meier survival curves (p < 0.05) are shown in Figures 9B–F, and further results are summarized in Supplementary Figure 3. The fractions of macrophages M2 cells (p = 0.05) and resting mast cells (p = 0.000403) were associated with improved survival (Figures 9B, C), but high levels of TFH cells (p = 0.05), Tregs (p = 0.00138), and B plasma cells (p = 0.00844) were associated with worse outcomes (Figures 9D–F).

Taken together, the higher expression of Gmfb in KIRC predicted a better prognosis, partially due to the optimal TME.

Discussion

KIRC is a cancer with high prevalence and mortality (3). However, it can be diagnosed early and successfully treated with surgical or ablative strategies; up to a third of cases present with or develop metastases (42). Thus, a deeper understanding of the molecular mechanisms underlying KIRC is needed. The present study integrated different KIRC databases with suitable meta-analysis approaches, identified GMFB as a novel independent prognostic biomarker, positively correlated with OS and DFS in KIRC patients, and further dissected the possible contribution of GMFB-relevant TME to good clinical outcomes in KIRC patients.

Two studies have shown that GMFB is a prognostic biomarker for astrocytomas (43) and serous ovarian cancers (44). However, using GEPIA2.0 and TCGA in pan-cancer analysis, we failed to obtain consistent results with previous reports (data not shown), possibly due to dynamic updates of the cancer database with increasing sample data deposition. We previously found that Gmfb KD significantly inhibited the proliferation, migration, and invasion of Hep3B cells, and higher Gmfb expression predicted a poor prognosis in LIHC (24). In the present study, a completely distinct relationship between Gmfb and OS/DFS in KIRC was obtained; namely, higher expression of GMFB was associated with a better prognosis for KIRC patients (Figure 3). To confirm the meta-analysis results, Gmfb KD and OE were performed in Caki-2 cells. We found that Gmfb KD promoted cell viability partly due to the downregulation of p21 and p27 protein levels in Caki-2 cells (Figures 4J, K). The downregulation of p21 and p27 protein levels indicated the transition of the cell cycle from quiescent to proliferation (41). Moreover, Gmfb OE inhibited the proliferation, migration, and decreased MMP in Caki-2 cells, and when MMP was disturbed, reactive oxygen species caused the oxidation of mitochondrial pores (45). Taken together, the in vitro results support that high expression of GMFB is beneficial for KIRC patients’ survival. The distinct molecular mechanisms underlying GMFB in LIHC and KIRC merit further investigation.

To reveal the molecular mechanism underlying Gmfb effects in KIRC, we enriched the DEGs with Gmfb co-expression for bioinformatics analysis. KEGG pathway and GO analyses were performed for the top 2,000 downregulated and upregulated DEGs. We validated nine significant DEGs. Consistent with the results from the meta-analysis, INAFM1, HCFC1R1, and METTL26 were downregulated, and ATP6V0D2, TBC1D14, OXCT1, ALDH6A1, FAM160A1, and NCR3LG1 were upregulated in GMFB overexpressed Caki-2 cells (Figures 5C, D). The validated DEGs with GMFB co-expression may participate in multiple pathways in KIRC biology, including enhancing response in KIRC cell to NK cell (NCR3LG1) (46), negatively regulating autophagy (TBC1D14) (47), and being regulated by lactate (ATP6V0D2) (48) in TMEs.

In the upregulated pathways, collecting duct acid secretion was the most enriched, containing DEGs ATP6V0A4, ATP6V1G3, ATP6V0D2, ATP6V1B1 (vacuolar proton pump subunit), SLC4A1 (chloride/bicarbonate exchanger), and CLCNKB (voltage-gated chloride channels). Lower CLCNKB expression has been reported in renal cancer (49). We hypothesized that the upregulated collecting duct acid secretion might play a role in maintaining acid–base homeostasis in the renal tubule.

Mineral absorption, neuroactive ligand–receptor interaction, and estrogen signaling in the downregulated panel were the three most enriched pathways (Figure 6E). Metallothionein (MT) proteins, such as MT1A, B, E, G, F, H, 349 HL1, X, M, MT2A, TF, and S100G, were enriched in the mineral absorption pathway. MTs have a high content of cysteine residues that bind various heavy metals. In this family, MT1A, B, E, G, F, H, HL1, X, M, MT2A, TF, and S100G (vitamin D-dependent calcium-binding protein) are enriched in the mineral absorption pathway. The dysregulation of MTs in renal cell cancer has also been reported. The toxic metals cadmium and lead are associated with renal cell cancer progression (50). A decrease in MT levels in renal cell cancer has also been reported (51). In the present study, we conducted qRT-PCR to validate some of the mRNA levels of several MTs in Caki-2 cells. The results showed that MT2A, MT1G, and MT1F were downregulated, but S100G and TF did not show a significant change in Gmfb OE Caki-2 cells (Figure 6I). This indicates that low expression of Gmfb may lead to increased circulating levels of toxic heavy metals, thus contributing to the progression of KIRC. Currently, there are no relevant reports on S100G and TF in KIRC, and their role, therefore, remains unclear. Cornification, keratinization, and the classical complement activation pathway were significantly enriched in the BP category of GO. There were 23 shared genes between the BP terms cornification and keratinization (Supplementary Table 7). Only CK7 (52, 53) and AE1AE3 (54) cytokeratins have been examined in KIRC. The fundamental functions of other cornification-related genes remain elusive, probably implicating an alteration of epithelial and mesenchymal immunophenotypes for KIRC, which requires further investigation. KOBAS-based GO enrichment analysis revealed downregulated DEGs in the extracellular region, extracellular space, and Ig complex circulating GO-CC category (Figure 6G). Most genes in the extracellular region (Supplementary Table 4) encode MMPs, cytokines, chemokines, growth factors, and others. These results, combined with the DEGs enriched in GO-MF (antigen binding, Ig receptor binding, and serine-type endopeptidase activity) (Figure 6H), suggest cell–cell interactions.

Analyzing PPI networks is increasingly recognized as meaningful to characterize the underlying biology of genes associated with complex diseases (55). As shown in Figure 8, histone cluster family members were enriched in CytoHubba with upregulated DEGs (Figure 8B), and the formation of the cornified envelope was enriched in CytoHubba with downregulated DEGs (Figure 8D). The former is possibly linked to epigenetic regulation in KIRC. In contrast, the latter may contribute to metastasis by disrupting the formation of the cornified envelope, which requires further validation.

The actual context of KIRC is quite diverse due to intratumoral or intertumoral TME heterogeneity. The enriched pathways and GO terms, including extracellular region and neuroactive ligand–receptor binding, encouraged us to analyze the Gmfb-relevant TME.

TME is a complex ecosystem, and it plays a crucial function in cancer progression and response to immunotherapy (7, 56). The TME is composed of adaptive immune cells (T and B lymphocytes) and innate immune response cells (including dendritic cells, mast cells, macrophages, neutrophils, myeloid-derived suppressor cells, and natural killer cells) as well as cancer cells and stroma (i.e., endothelial cells, fibroblasts, pericytes, and mesenchymal cells) (56). Immune and inflammatory responses are strongly associated with survival outcomes in patients with KIRC. Zhang et al. recently revealed insights into the TME in KIRC using single-cell sequencing technology (57) and elucidated an active role for tumor epithelia in promoting immune cell infiltration. In contrast, the immune cell composition in the TME may impact clinical outcomes. Infiltration by CD8+ T cells is associated with a worse prognosis in KIRC (58). CCR4 has been identified as a TME target for renal cancer (59). In the present study, we further integrated the OS data obtained from GEPIA2021 with the results of CIBERSORTx. We found a significant association between GMFB-associated immune cells and OS in patients with KIRC. A significant correlation between GMFB expression and KIRC is shown in Figure 9A. GMFB expression was positively correlated with macrophage M2 cells and resting mast cells but negatively associated with T follicular helper cells, Tregs, and plasma cells (Figure 9A). In addition, high levels of M2 macrophages (p = 0.05) and resting mast cells (p = 0.000403) favored a better prognosis (Figures 9B, C), while high levels of TFH cells (p = 0.05), Tregs (p = 0.00138), and plasma cells (p = 0.00844) caused deleterious clinical outcomes in patients with KIRC (Figures 9D–F). We summarized the role of Gmfb and Gmfb-related TME in KIRC in Figure 10.

Figure 10

The Gmfb-related TME in LICH was analyzed. Unexpectedly, the Gmfb-related TME in LIHC was quite different from that in KIRC (Supplementary Figure 2). Thus, we speculated that the distinct effects of Gmfb expression on the OS of LIHC and KIRC patients might be due to different immune cell infiltration.

This study had several limitations. The first one was mainly based on bioinformatics analysis, especially datasets from bulk RNA-seq, and lacked validation experiments with human KIRC samples. Secondly, the Gmfb-related TME also needed to be validated with human KIRC samples. The critical molecular signatures identified in KIRC merit further clinical investigation.

In conclusion, our study demonstrated that GMFB could be considered a novel independent prognostic biomarker for KIRC. Modulation of the TME in KIRC by GMFB intervention may represent a novel immunotherapeutic strategy for KIRC and probably improve clinical outcomes.

Funding

This paper was supported by the Ministry of Science and Technology of China (2020YFA0113101), the National Natural Science Foundation of China (32070719, 81770942), and the Fundamental Research Funds for the Central Universities (22120220009).

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

• 1

  SiegelRLMillerKDJemalA. Cancer Statistics, 2019. CA Cancer J Clin (2019) 69(1):7–34. doi: 10.3322/caac.21551

  • CrossRef
  • Google Scholar

• 2

  JiangHChenHWanPChenN. Decreased Expression of HADH is Related to Poor Prognosis and Immune Infiltration in Kidney Renal Clear Cell Carcinoma. Genomics (2021) 113(6):3556–64. doi: 10.1016/j.ygeno.2021.08.008

  • CrossRef
  • Google Scholar

• 3

  MitchellTJTurajlicSRowanANicolDFarmeryJHRO'BrienTet al. Timing the Landmark Events in the Evolution of Clear Cell Renal Cell Cancer: TRACERx Renal. Cell (2018) 173(3):611–23.e17. doi: 10.1016/j.cell.2018.02.020

  • CrossRef
  • Google Scholar

• 4

  ShuchBAminAArmstrongAJEbleJNFicarraVLopez-BeltranAet al. Understanding Pathologic Variants of Renal Cell Carcinoma: Distilling Therapeutic Opportunities From Biologic Complexity. Eur Urol (2015) 67(1):85–97. doi: 10.1016/j.eururo.2014.04.029

  • CrossRef
  • Google Scholar

• 5

  Pena-LlopisSVega-Rubín-de-CelisSLiaoALengNPavía-JiménezAWangSet al. BAP1 Loss Defines a New Class of Renal Cell Carcinoma. Nat Genet (2012) 44(7):751–9. doi: 10.1038/ng.2323

  • CrossRef
  • Google Scholar

• 6

  KohMYSayeghNAgarwalN. Seeing the Forest for the Trees-Single-Cell Atlases Link CD8(+) T Cells and Macrophages to Disease Progression and Treatment Response in Kidney Cancer. Cancer Cell (2021) 39(5):594–6. doi: 10.1016/j.ccell.2021.03.008

  • CrossRef
  • Google Scholar

• 7

  BraunDAStreetKBurkeKPCookmeyerDLDenizeTPedersenCBet al. Progressive Immune Dysfunction With Advancing Disease Stage in Renal Cell Carcinoma. Cancer Cell (2021) 39(5):632–48.e8. doi: 10.1016/j.ccell.2021.02.013

  • CrossRef
  • Google Scholar

• 8

  HakimiAAOstrovnayaIJacobsenASusztakKColemanJARussoPet al. Validation and Genomic Interrogation of the MET Variant Rs11762213 as a Predictor of Adverse Outcomes in Clear Cell Renal Cell Carcinoma. Cancer (2016) 122(3):402–10. doi: 10.1002/cncr.29765

  • CrossRef
  • Google Scholar

• 9

  SingerEABratslavskyGMiddeltonLSrinivasanRLinehanWM. Impact of Genetics on the Diagnosis and Treatment of Renal Cancer. Curr Urol Rep (2011) 12(1):47–55. doi: 10.1007/s11934-010-0156-y

  • CrossRef
  • Google Scholar

• 10

  DalglieshGLFurgeKGreenmanCChenLBignellGButlerAet al. Systematic Sequencing of Renal Carcinoma Reveals Inactivation of Histone Modifying Genes. Nature (2010) 463(7279):360–3. doi: 10.1038/nature08672

  • CrossRef
  • Google Scholar

• 11

  Martinez-SalamancaJIHuangWCMillánIBertiniRBiancoFJCarballidoJAet al. Prognostic Impact of the 2009 UICC/AJCC TNM Staging System for Renal Cell Carcinoma With Venous Extension. Eur Urol (2011) 59(1):120–7. doi: 10.1016/j.eururo.2010.10.001

  • CrossRef
  • Google Scholar

• 12

  DelahuntBDagherJEgevadLYaxleyJVarmaMSamaratungaHet al. Is the UICC/AJCC Pt2 Staging Category for Clear Cell Renal Cell Carcinoma Meaningful? Am J Surg Pathol (2019) 43(9):1249–52. doi: 10.1097/PAS.0000000000001308

  • CrossRef
  • Google Scholar

• 13

  LiDLiuSXuJChenLXuCChenFet al. Ferroptosis-Related Gene CHAC1 is a Valid Indicator for the Poor Prognosis of Kidney Renal Clear Cell Carcinoma. J Cell Mol Med (2021) 25(7):3610–21. doi: 10.1111/jcmm.16458

  • CrossRef
  • Google Scholar

• 14

  WangQZhangHChenQWanZGaoXQianWet al. Identification of METTL14 in Kidney Renal Clear Cell Carcinoma Using Bioinformatics Analysis. Dis Markers (2019) 2019:5648783. doi: 10.1155/2019/5648783

  • CrossRef
  • Google Scholar

• 15

  XieLLiHZhangLMaXDangYGuoJet al. Autophagy-Related Gene P4HB: A Novel Diagnosis and Prognosis Marker for Kidney Renal Clear Cell Carcinoma. Aging (Albany NY) (2020) 12(2):1828–42. doi: 10.18632/aging.102715

  • CrossRef
  • Google Scholar

• 16

  SongJSongFLiuKZhangWLuoRTangYet al. Multi-Omics Analysis Reveals Epithelial-Mesenchymal Transition-Related Gene FOXM1 as a Novel Prognostic Biomarker in Clear Cell Renal Carcinoma. Aging (Albany NY) (2019) 11(22):10316–37. doi: 10.18632/aging.102459

  • CrossRef
  • Google Scholar

• 17

  AlbiniASpornMB. The Tumor Microenvironment as a Target for Chemoprevention. Nat Rev Cancer (2007) 7(2):139–47. doi: 10.1038/nrc2067

  • CrossRef
  • Google Scholar

• 18

  HuJChenZBaoLZhouLHouYLiuLet al. Single-Cell Transcriptome Analysis Reveals Intratumoral Heterogeneity in ccRCC, Which Results in Different Clinical Outcomes. Mol Ther (2020) 28(7):1658–72. doi: 10.1016/j.ymthe.2020.04.023

  • CrossRef
  • Google Scholar

• 19

  InagakiMAoyamaMSobueKYamamotoNMorishimaTMoriyamaAet al. Sensitive Immunoassays for Human and Rat GMFB and GMFG, Tissue Distribution and Age-Related Changes. Biochim Biophys Acta (2004) 1670(3):208–16. doi: 10.1016/j.bbagen.2003.12.006

  • CrossRef
  • Google Scholar

• 20

  ZuoPFuZTaoTYeFChenLWangXet al. The Expression of Glia Maturation Factors and the Effect of Glia Maturation Factor-Gamma on Angiogenic Sprouting in Zebrafish. Exp Cell Res (2013) 319(5):707–17. doi: 10.1016/j.yexcr.2013.01.004

  • CrossRef
  • Google Scholar

• 21

  YinGDuMLiRLiKHuangXDuanDet al. Glia Maturation Factor Beta is Required for Reactive Gliosis After Traumatic Brain Injury in Zebrafish. Exp Neurol (2018) 305:129–38. doi: 10.1016/j.expneurol.2018.04.008

  • CrossRef
  • Google Scholar

• 22

  ZaheerALimR. Protein Kinase A (PKA)- and Protein Kinase C-Phosphorylated Glia Maturation Factor Promotes the Catalytic Activity of PKA. J Biol Chem (1997) 272(8):5183–6. doi: 10.1074/jbc.272.8.5183

  • CrossRef
  • Google Scholar

• 23

  WernerCJHeyny-von HaussenRMallGWolfS. Proteome Analysis of Human Substantia Nigra in Parkinson's Disease. Proteome Sci (2008) 6:8. doi: 10.1186/1477-5956-6-8

  • CrossRef
  • Google Scholar

• 24

  SunWHuCWangTWangJZhangJGaoFet al. Glia Maturation Factor Beta as a Novel Biomarker and Therapeutic Target for Hepatocellular Carcinoma. Front Oncol (2021) 11:744331. doi: 10.3389/fonc.2021.744331

  • CrossRef
  • Google Scholar

• 25

  HutterCZenklusenJC. The Cancer Genome Atlas: Creating Lasting Value Beyond Its Data. Cell (2018) 173(2):283–5. doi: 10.1016/j.cell.2018.03.042

  • CrossRef
  • Google Scholar

• 26

  ChandrashekarDSBashelBBalasubramanyaSAHCreightonCJPonce-RodriguezIChakravarthiBVSKet al. UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia (2017) 19(8):649–58. doi: 10.1016/j.neo.2017.05.002

  • CrossRef
  • Google Scholar

• 27

  ThulPJLindskogC. The Human Protein Atlas: A Spatial Map of the Human Proteome. Protein Sci (2018) 27(1):233–44. doi: 10.1002/pro.3307

  • CrossRef
  • Google Scholar

• 28

  TangZKangBLiCChenTZhangZ. GEPIA2: An Enhanced Web Server for Large-Scale Expression Profiling and Interactive Analysis. Nucleic Acids Res (2019) 47(W1):W556–60. doi: 10.1093/nar/gkz430

  • CrossRef
  • Google Scholar

• 29

  BuDLuoHHuoPWangZZhangSHeZet al. KOBAS-I: Intelligent Prioritization and Exploratory Visualization of Biological Functions for Gene Enrichment Analysis. Nucleic Acids Res (2021) 49(W1):W317–25. doi: 10.1093/nar/gkab447

  • CrossRef
  • Google Scholar

• 30

  AiCKongL. CGPS: A Machine Learning-Based Approach Integrating Multiple Gene Set Analysis Tools for Better Prioritization of Biologically Relevant Pathways. J Genet Genomics (2018) 45(9):489–504. doi: 10.1016/j.jgg.2018.08.002

  • CrossRef
  • Google Scholar

• 31

  WuJMaoXCaiTLuoJWeiL. KOBAS Server: A Web-Based Platform for Automated Annotation and Pathway Identification. Nucleic Acids Res (2006) 34(Web Server issue):W720–4. doi: 10.1093/nar/gkl167

  • CrossRef
  • Google Scholar

• 32

  XieCMaoXHuangJDingYWuJDongSet al. KOBAS 2.0: A Web Server for Annotation and Identification of Enriched Pathways and Diseases. Nucleic Acids Res (2011) 39(Web Server issue):W316–22. doi: 10.1093/nar/gkr483

  • CrossRef
  • Google Scholar

• 33

  YuanSXuYYiTWangH. The Anti-Tumor Effect of OP-B on Ovarian Cancer In Vitro and In Vivo, and its Mechanism: An Investigation Using Network Pharmacology-Based Analysis. J Ethnopharmacol (2021) p:114706. doi: 10.1016/j.jep.2021.114706

  • CrossRef
  • Google Scholar

• 34

  ZhouYZhouBPacheLChangMKhodabakhshiAHTanaseichukOet al. Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat Commun (2019) 10(1):1523. doi: 10.1038/s41467-019-09234-6

  • CrossRef
  • Google Scholar

• 35

  ChenBKhodadoustMSLiuCLNewmanAMAlizadehAA. Profiling Tumor Infiltrating Immune Cells With CIBERSORT. Methods Mol Biol (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12

  • CrossRef
  • Google Scholar

• 36

  NewmanAMLiuCLGreenMRGentlesAJFengWXuYet al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

  • CrossRef
  • Google Scholar

• 37

  LiCTangZZhangWYeZLiuF. GEPIA2021: Integrating Multiple Deconvolution-Based Analysis Into GEPIA. Nucleic Acids Res (2021) 49(W1):W242–6. doi: 10.1093/nar/gkab418

  • CrossRef
  • Google Scholar

• 38

  SzklarczykDGableALLyonDJungeAWyderSHuerta-CepasJet al. STRING V11: Protein-Protein Association Networks With Increased Coverage, Supporting Functional Discovery in Genome-Wide Experimental Datasets. Nucleic Acids Res (2019) 47(D1):D607–13. doi: 10.1093/nar/gky1131

  • CrossRef
  • Google Scholar

• 39

  ChinCHChenSHWuHHHoCWKoMTLinCY. Cytohubba: Identifying Hub Objects and Sub-Networks From Complex Interactome. BMC Syst Biol (2014) 8 Suppl 4:S11. doi: 10.1186/1752-0509-8-S4-S11

  • CrossRef
  • Google Scholar

• 40

  ShannonPMarkielAOzierOBaligaNSWangJTRamageDet al. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

  • CrossRef
  • Google Scholar

• 41

  KajiNOhashiKShuinMNiwaRUemuraTMizunoKet al. Cell Cycle-Associated Changes in Slingshot Phosphatase Activity and Roles in Cytokinesis in Animal Cells. J Biol Chem (2003) 278(35):33450–5. doi: 10.1074/jbc.M305802200

  • CrossRef
  • Google Scholar

• 42

  JonaschEGaoJRathmellWK. Renal Cell Carcinoma. BMJ (2014) 349:g4797. doi: 10.1136/bmj.g4797

  • CrossRef
  • Google Scholar

• 43

  KuangXYXuSLCuiWJiangXF. Association of GMF-Beta Expression With Ki-67 and its Significance in the Prognostic Evaluation of Astrocytoma. Zhonghua Bing Li Xue Za Zhi (2021) 50(11):1252–6. doi: 10.3760/cma.j.cn112151-20210627-00459

  • CrossRef
  • Google Scholar

• 44

  LiYLYeFChengXDHuYZhouCYLüWGet al. Identification of Glia Maturation Factor Beta as an Independent Prognostic Predictor for Serous Ovarian Cancer. Eur J Cancer (2010) 46(11):2104–18. doi: 10.1016/j.ejca.2010.04.015

  • CrossRef
  • Google Scholar

• 45

  Garcia-MartinezTVendrell-FlotatsMMartínez-RoderoIOrdóñez-LeónEAÁlvarez-RodríguezMLópez-BéjarMet al. Glutathione Ethyl Ester Protects In Vitro-Maturing Bovine Oocytes Against Oxidative Stress Induced by Subsequent Vitrification/Warming. Int J Mol Sci (2020) 21(20):7547. doi: 10.3390/ijms21207547

  • CrossRef
  • Google Scholar

• 46

  ShefferMLowryEBeelenNBorahMAmaraSNMaderCCet al. Genome-Scale Screens Identify Factors Regulating Tumor Cell Responses to Natural Killer Cells. Nat Genet (2021) 53(8):1196–206. doi: 10.1038/s41588-021-00889-w

  • CrossRef
  • Google Scholar

• 47

  LuTLiYPanMYuDWangZLiuCet al. TBC1D14 Inhibits Autophagy to Suppress Lymph Node Metastasis in Head and Neck Squamous Cell Carcinoma by Downregulating Macrophage Erythroblast Attacher. Int J Biol Sci (2022) 18(5):1795–812. doi: 10.7150/ijbs.68992

  • CrossRef
  • Google Scholar

• 48

  LiuNLuoJKuangDXuSDuanYXiaYet al. Lactate Inhibits ATP6V0d2 Expression in Tumor-Associated Macrophages to Promote HIF-2alpha-Mediated Tumor Progression. J Clin Invest (2019) 129(2):631–46. doi: 10.1172/JCI123027

  • CrossRef
  • Google Scholar

• 49

  MurakamiTSanoFHuangYKomiyaABabaMOsadaYet al. Identification and Characterization of Birt-Hogg-Dube Associated Renal Carcinoma. J Pathol (2007) 211(5):524–31. doi: 10.1002/path.2139

  • CrossRef
  • Google Scholar

• 50

  TarabayHHAbol-EneinHAwadallaAMortadaWIAbdel-AzizAF. Gene Expression and Oxidative Stress Markers Profile Associated With Toxic Metals in Patients With Renal Cell Carcinoma. Mol Biol Rep (2022) 49(2):1161–9. doi: 10.1007/s11033-021-06944-3

  • CrossRef
  • Google Scholar

• 51

  IshiiKUsuiSYamamotoHSugimuraYTatematsuMHiranoKet al. Decreases of Metallothionein and Aminopeptidase N in Renal Cancer Tissues. J Biochem (2001) 129(2):253–8. doi: 10.1093/oxfordjournals.jbchem.a002852

  • CrossRef
  • Google Scholar

• 52

  TruongLDShenSS. Immunohistochemical Diagnosis of Renal Neoplasms. Arch Pathol Lab Med (2011) 135(1):92–109. doi: 10.5858/2010-0478-RAR.1

  • CrossRef
  • Google Scholar

• 53

  WilliamsonSRGaddeRTrpkovKHirschMSSrigleyJRReuterVEet al. Diagnostic Criteria for Oncocytic Renal Neoplasms: A Survey of Urologic Pathologists. Hum Pathol (2017) 63:149–56. doi: 10.1016/j.humpath.2017.03.004

  • CrossRef
  • Google Scholar

• 54

  StepanAEMărgčritescuCStoicaLEStepanMDSimionescuCE. Clear Cell Renal Cell Carcinomas - Epithelial and Mesenchymal Immunophenotype. Rom J Morphol Embryol (2018) 59(4):1189–94.

  • Google Scholar

• 55

  GoniJEstebanFJde MendizábalNVSepulcreJArdanza-TrevijanoSAgirrezabalIet al. A Computational Analysis of Protein-Protein Interaction Networks in Neurodegenerative Diseases. BMC Syst Biol (2008) 2:52. doi: 10.1186/1752-0509-2-52

  • CrossRef
  • Google Scholar

• 56

  GrivennikovSIGretenFRKarinM. Immunity, Inflammation, and Cancer. Cell (2010) 140(6):883–99. doi: 10.1016/j.cell.2010.01.025

  • CrossRef
  • Google Scholar

• 57

  ZhangYNarayananSPMannanRRaskindGWangXVatsPet al. Single-Cell Analyses of Renal Cell Cancers Reveal Insights Into Tumor Microenvironment, Cell of Origin, and Therapy Response. Proc Natl Acad Sci U S A (2021) 118(24):e2103240118. doi: 10.1073/pnas.2103240118

  • CrossRef
  • Google Scholar

• 58

  GiraldoNABechtEVanoYPetitprezFLacroixLValidirePet al. Tumor-Infiltrating and Peripheral Blood T-Cell Immunophenotypes Predict Early Relapse in Localized Clear Cell Renal Cell Carcinoma. Clin Cancer Res (2017) 23(15):4416–28. doi: 10.1158/1078-0432.CCR-16-2848

  • CrossRef
  • Google Scholar

• 59

  BerlatoCKhanMNSchioppaTThompsonRManiatiEMontfortAet al. A CCR4 Antagonist Reverses the Tumor-Promoting Microenvironment of Renal Cancer. J Clin Invest (2017) 127(3):801–13. doi: 10.1172/JCI82976

  • CrossRef
  • Google Scholar