Skip to main content

ORIGINAL RESEARCH article

Front. Oncol., 30 September 2021
Sec. Genitourinary Oncology
This article is part of the Research Topic Biomarkers in Genitourinary Cancers, Volume I View all 21 articles

Identification of Hub Genes Associated With Clear Cell Renal Cell Carcinoma by Integrated Bioinformatics Analysis

  • 1Department of Nephrology, Xiangya Hospital Central South University, Changsha, China
  • 2Department of Cell Biology, School of Life Sciences, Central South University, Changsha, China
  • 3Hunan Key Laboratory of Organ Fibrosis, Central South University, Changsha, China
  • 4Department of Otolaryngology-Head and Neck Surgery, Second Xiangya Hospital Central South University, Changsha, China

Background: Clear cell renal cell carcinoma (ccRCC) is a common genitourinary cancer type with a high mortality rate. Due to a diverse range of biochemical alterations and a high level of tumor heterogeneity, it is crucial to select highly validated prognostic biomarkers to be able to identify subtypes of ccRCC early and apply precision medicine approaches.

Methods: Transcriptome data of ccRCC and clinical traits of patients were obtained from the GSE126964 dataset of Gene Expression Omnibus and The Cancer Genome Atlas Kidney Renal Clear Cell Carcinoma (TCGA-KIRC) database. Weighted gene co-expression network analysis (WGCNA) and differentially expressed gene (DEG) screening were applied to detect common differentially co-expressed genes. Gene Ontology, Kyoto Encyclopedia of Genes and Genomes analysis, survival analysis, prognostic model establishment, and gene set enrichment analysis were also performed. Immunohistochemical analysis results of the expression levels of prognostic genes were obtained from The Human Protein Atlas. Single-gene RNA sequencing data were obtained from the GSE131685 and GSE171306 datasets.

Results: In the present study, a total of 2,492 DEGs identified between ccRCC and healthy controls were filtered, revealing 1,300 upregulated genes and 1,192 downregulated genes. Using WGCNA, the turquoise module was identified to be closely associated with ccRCC. Hub genes were identified using the maximal clique centrality algorithm. After having intersected the hub genes and the DEGs in GSE126964 and TCGA-KIRC dataset, and after performing univariate, least absolute shrinkage and selection operator, and multivariate Cox regression analyses, ALDOB, EFHD1, and ESRRG were identified as significant prognostic factors in patients diagnosed with ccRCC. Single-gene RNA sequencing analysis revealed the expression profile of ALDOB, EFHD1, and ESRRG in different cell types of ccRCC.

Conclusions: The present results demonstrated that ALDOB, EFHD1, and ESRRG may act as potential targets for medical therapy and could serve as diagnostic biomarkers for ccRCC.

Introduction

Renal cell carcinoma (RCC) is one of the most common genitourinary cancer types worldwide, and it has a number of heterogeneous histological subtypes, with clear cell RCC (ccRCC) accounting for ~85% of all cases (1). In total, 431,288 new patients were diagnosed with renal cancer, and 179,368 of these patients succumbed to the disease worldwide in 2020 (2). ccRCC is not susceptible to chemoradiotherapy (3). Although ccRCC is curable at an early localized stage by partial or total surgical nephrectomy, advanced or metastatic ccRCC remains a clinical challenge (4). Over the past years, antiangiogenic treatment, inhibitors of the mammalian target of rapamycin (mTOR) pathway, or immune checkpoint inhibition therapy have considerably evolved (5). However, due to diverse biochemical alterations and a high level of tumor heterogeneity, it is important to select highly validated prognostic biomarkers to identify subtypes of ccRCC early and apply precision medicine approaches (6).

The molecular mechanism of ccRCC is characterized by genetic diversity and chromosomal complexity. Loss of the heterozygosity of chromosome 3p, where the von Hippel–Lindau (VHL) gene is located, is found in over 90% of ccRCC cases, and it is considered the critical genetic event (79). A loss-of-function mutation in the VHL gene induces the aberrant regulation of a number of VHL-mediated targets, pathways, and processes, which is a significant step in the development of ccRCC (10, 11). The VHL protein, as an E3 ubiquitin ligase, is notably involved in the ubiquitylation of the prolyl hydroxylated transcription factors, hypoxia-inducible factor 1 α (HIF1 α) and HIF2 α, under normoxic conditions. HIF1 α and HIF2 α have an important role in the regulation of angiogenesis, erythropoiesis, glycolysis, and apoptosis (1214). Moreover, next-generation sequencing technologies have provided evidence that PBRM1, SETD2, or BAP1 mutations are the drivers of tumor evolution (15, 16). Although the molecular features of ccRCC have been increasingly defined by previous studies (1719), there remain numerous subtypes of ccRCC the pathogenic mechanisms of which have yet to be clearly determined at the genetic and molecular levels. Thus, it is important to identify more additional disease-related genes.

Benefiting from the rapid development of genome sequencing technology, bioinformatics can be used to study gene expression profiles in order to examine the molecular mechanism of tumors and identify tumor-specific indicators. Weighted gene co-expression network analysis (WGCNA) was developed by Horvath and Zhang in 2005 (20). At present, WGCNA is becoming a powerful approach to detecting gene modules, exploring the correlation of the modules and phenotypes, and discovering hub genes that regulate critical biological processes (21, 22).

In the present study, a gene expression profile of ccRCC from the Gene Expression Omnibus (GEO) was downloaded. WGCNA and differentially expressed gene (DEG) screening were applied to detect common differentially co-expressed genes. Then, The Cancer Genome Atlas Kidney Renal Clear Cell Carcinoma (TCGA-KIRC) data were used to establish the prognostic model of ccRCC. Single-cell RNA sequencing (RNA-seq) data from GEO were used to verify the expression profile of the prognostic genes in different cell types. This study aimed not only to understand ccRCC pathogenesis but also to determine its molecular mechanisms and provide insights into novel therapeutic targets for drugs.

Methods

Data Collection and Single-Cell RNA Sequencing Data Processing

The workflow for the current study is presented in Figure 1. Original data were collected from the GSE126964 dataset, which contained 55 ccRCC tumor tissues and 11 matched normal tissues (23). The GEO expression matrix was annotated with gene symbols using the information from the GPL20795 HiSeq X Ten platform file, as well as log2 transformed in R (version 4.0.4) and RStudio (version 1.2.5033) if necessary. Principal component analysis (PCA) was performed, and the outliers of GSM3619137 and GSM3619152 were excluded (Figure S1). In total, only 53 ccRCC sample and 11 normal sample data were used for subsequent analysis.

FIGURE 1
www.frontiersin.org

Figure 1 Workflow of the present study. (A) Identification workflow. (B) Verification workflow. DEG, differentially expressed gene; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; TCGA, The Cancer Genome Atlas; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; HPA, Human Protein Atlas; GSEA, gene set enrichment analysis.

RNA-seq data of TCGA-KIRC and corresponding clinical information were obtained from TCGA (https://portal.gdc.cancer.gov/).

Single-cell RNA-seq data from GSE131685 and GSE171306 were downloaded through GEO website. R package “Seurat” (version 4.0.2) was used to process the data (24). Three healthy kidney samples from GSE131685 (25) and two ccRCC samples from GSE171306 (26) were merged for further analysis. The single-cell RNA-seq data processing was described previously (27). The cell clusters were annotated manually based on previous knowledge and information from literatures (28, 29). Expression profiling of the genes were depicted by heatmap and violin plot using the function “FeaturePlot” and “VlnPlot.”

Differentially Expressed Gene Identification

The “limma” software package (version 3.48.0) (30) was used to conduct the DEG analysis between ccRCC and normal sample data from the GSE126964 dataset (30). An adjusted p-value <0.05 and a |log2fold change| of ≥2.0 were selected as the cutoff criteria. The volcano and heatmap plots were generated using ggplot2 (version 3.3.3) and pheatmap (version 1.0.12) packages, respectively. The DEGs of TCGA-KIRC dataset (https://portal.gdc.cancer.gov/) were obtained via Gene Expression Profiling Interactive Analysis (GEPIA2; http://gepia2.cancer-pku.cn/) (31) using the same aforementioned threshold.

Gene Ontology Enrichment and Kyoto Encyclopedia of Genes and Genomes Analysis of Differentially Expressed Genes

The “clusterProfiler” (version 3.18.1) R package was used for GO and KEGG enrichment analyses (http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (32). The three main processes in GO analysis are as follows: biological process (BP), molecular function (MF), and cellular component (CC). The p-value was conventionally set at 0.05. A circle plot was generated by “Goplot” R package (version 1.0.2).

Weighted Gene Co-Expression Network Construction

The “WGCNA” package (version 1.70-3) of R (20) was used to construct the co-expression networks. Genes with mean counts of over 5 were selected. A total of 64 samples were used to calculate the Pearson’s correlation matrices. The matrices of adjacency were created based on the Pearson’s correlation matrices. Then, the clinical trait data were uploaded, and the scale independence and mean connectivity were estimated. Subsequently, the topological overlap measure (TOM) matrix, which was created from the adjacency matrix, was used to estimate the network’s connectivity property. A hierarchical clustering dendrogram of the TOM matrix was constructed using the average distance with a minimum size threshold of 50 to classify the similar gene expression profiles into different gene modules. Finally, similar gene modules were merged, with a threshold of 0.20.

Co-Expression Network Construction and Hub Gene Identification

The Cytoscape software v3.7.2 was used to visualize the co-expression network in the turquoise module (33). The data were imported into Cytohubba, a Cytoscape plug-in for hub gene identification, and the maximal clique centrality (MCC) algorithm was used to calculate the scores of all nodes of the network. The top 30 nodes with the highest MCC scores were selected as the hub genes associated with ccRCC. The “real” key genes were identified as those intersecting between the top 30 nodes in turquoise module, DEGs from GSE126964 and DEGs from TCGA-KIRC.

Identification and Verification of Prognostic Gene Signatures

Univariate Cox regression analysis was performed to screen the genes significantly associated with overall survival (OS) in the TCGA-KIRC dataset. The OS-related genes with p < 0.1 were included in the least absolute shrinkage and selection operator (LASSO) regression analysis by using the R package “glmnet” (version 4.1-2). Then, a multivariate Cox regression model analysis was performed to establish a Cox proportional hazards regression prognostic model. We used the following formula to calculate the risk score of each patient:

Risk Score=i=1nβi×xi

In this formula, β is coefficient and x is the expression level of each prognostic gene i. The samples were divided into a high-risk group and a low-risk group according to the median risk score of the training cohort from TCGA-KIRC. Receiver operating characteristic (ROC) analysis and Kaplan–Meier analysis were conducted between the high-risk group and the low-risk group.

Validation of the Protein Expression Levels of Prognostic Genes in the Human Protein Atlas Database

The Human Protein Atlas (HPA) is a database that aims to map all the human proteins in cells, tissues, and organs using an integration of various omics technologies (https://www.proteinatlas.org/). We also verified the protein expression levels of the survival-related hub genes based on immunohistochemistry using the HPA database.

Gene Set Enrichment Analysis of Prognostic Genes

Gene set enrichment analysis (GSEA) was also used to detect the potential molecular mechanisms of the prognostic genes. Enriched terms predicted to be associated with the KEGG pathway in c2.cp.v7.2.symbols.gmt were screened by GSEA. Images were generated by “ggplot2” (version 3.3.3) package. The p-value of <0.05 was considered statistically significant.

Prognostic Gene Expression Profiles

The prognostic gene expression profiles were obtained from the GTEx Portal (https://gtexportal.org/home/).

Results

Differentially Expressed Gene Screening

The “limma” package was utilized to analyze DEGs in the GSE126964 dataset, with the threshold of |log2(fold-change)|>2.0 and adjusted p < 0.05. A total of 2,492 DEGs between ccRCC and normal control samples were filtered, revealing 1,300 upregulated genes and 1,192 downregulated genes (Figures 2A, B).

FIGURE 2
www.frontiersin.org

Figure 2 Screening for DEGs. (A) Volcano map of DEGs between ccRCC and normal samples in the GSE126964 dataset. The red plots in the volcano represent upregulated genes, and the blue points represent downregulated genes. (B) Heatmap of the 200 selected DEGs. The color in heatmaps from blue to yellow shows the progression from low expression to high expression, respectively. (C) GO and KEGG analyses of the DEGs. The outer circle shows the scatter plot of the assigned gene log2fold change for all terms: red points show genes that exhibited increased expression, whereas the blue points represent genes that exhibited decreased expression. The inner circle indicates the Z-score value and the number of genes. Red represents a higher z-score value, and purple represents a lower Z-score value. DEG, differentially expressed gene; BP, biological process; CC, cell component; MF, molecular function; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ccRCC, clear cell renal cell carcinoma.

The DEGs were mostly enriched in “T cell activation,” “leukocyte cell–cell adhesion,” “apical part of cell,” “external side of plasma membrane,” “collagen-containing extracellular matrix (ECM),” and “ion transmembrane transporter activity” in the GO analysis (Figure 2C). In the KEGG analysis, DEGs were enriched in “cytokine–cytokine receptor interaction,” “hematopoietic cell lineage,” “viral protein interaction with cytokine and cytokine receptor,” “cell adhesion molecules,” and “protein digestion and absorption” (Figure 2C).

We also evaluated the metabolic shift between ccRCC tissues and normal control tissues in the GSE126964 dataset. Similar to the finding of Clark et al. (19), glycolysis-associated genes were found to be significantly upregulated, and most oxidative phosphorylation (OXPHOS) and tricarboxylic acid (TCA) cycle-associated genes were significantly downregulated in the GSE126964 dataset (Figure S2).

Weighted Co-Expression Network Construction and Analysis

The sample clustering dendrograms of the ccRCC and normal samples are shown in Figure S3A. The soft-power threshold β was selected as 5 to ensure that both the scale-free topology model fit index (R2) and mean connectivity reached steady status (Figure 3A). Then, gene modules were detected based on the TOM matrix. A total of 25 modules were identified via average linkage hierarchical clustering, and each module was represented by a different color (Figure 3B). Among the modules, the turquoise module had the highest correlation with ccRCC traits (r = -0.97, p = 1e-39) (Figure 3C). A set of 400 selected genes were identified for the network heatmap construction (Figure S3B).

FIGURE 3
www.frontiersin.org

Figure 3 WGCNA of the ccRCC samples. (A) Analysis of the network topology for various soft-thresholding powers. The left plot shows the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). The horizontal red line shows x-axis = 0.85. The right plot displays the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis). The power was set as 5 for further analysis. (B) Hierarchical cluster analysis was conducted to detect co-expression clusters with corresponding color assignments. Each color represents a module in the constructed gene co-expression network, as assessed via WGCNA. (C) Module–trait relationships. Each row represents a color module, and every column represents a clinical trait. Each cell contains the corresponding correlation and p-value. (D) A scatter plot of GS for ccRCC vs. the MM in the turquoise module. (E) Identification of hub genes using the MCC method. Genes with the top 30 MCC values were colored red to yellow. Red refers to a relatively large MCC value, and yellow refers to relatively smaller MCC values. WGCNA, weighted gene co-expression network analysis; ccRCC, clear cell renal cell carcinoma; GS, gene significance; MM, module membership; MCC, maximal clique centrality.

Identification of Key Genes

An intramodular analysis of gene significance (GS) and module membership (MM) of the genes in the module turquoise was subsequently conducted. A high correlation coefficient of GS and MM was found in the turquoise module (cor = 0.97, p < 1e-200) (Figure 3D). The co-expression network of the turquoise module was constructed using Cytoscape software. Then, the module net was analyzed with the “Cytohubba” plug-in, and a network of the top 30 hub genes was constructed using the MCC algorithm (Figure 3E).

In order to identify the “real” key genes, we then obtained 796 DEGs, using a cohort of KIRC, from TCGA via GEPIA2, with the same threshold values. After comparing the DEGs in the GSE126964 dataset, TCGA-KIRC data and the top 30 hub genes from the turquoise module, a set of 13 key genes was identified (Figure 4A).

FIGURE 4
www.frontiersin.org

Figure 4 Prognostic analysis of the key genes. (A) Key genes belonging to both the hub genes and DEGs of the GSE126964 and TCGA-KIRC datasets. (B, C) LASSO regression complexity was controlled by lambda using the glmnet R package. (D) The multivariate analysis of risk factors in ccRCC. (E) Overall survival and ROC analysis between high-risk score and low-risk score groups in the training cohort and testing cohort. (F) The overall survival stratified by the high- and low-risk score groups was plotted for the training cohort and testing cohort. Detailed risk scores, survival information, and heat maps of gene expression are also included for each dataset. *p < 0.05; **p < 0.01; DEG, differentially expressed gene; TCGA-KIRC, The Cancer Genome Atlas Kidney Renal Clear Cell Carcinoma; ROC, receiver operating characteristic; AUC, area under the curve; ccRCC, clear cell renal cell carcinoma; LASSO, least absolute shrinkage and selection operator.

Validation of Key Genes via Survival Analysis

We randomly divided the patients in TCGA-KIRC into two cohorts, a training cohort (N = 266) and a testing cohort (N = 266). The univariate Cox regression analyses of 13 key genes with regard to OS of samples from the training cohort were performed (Table 1). Eight genes with p < 0.1 (GGT6, SLC22A8, FAM3B, PTH1R, ALDOB, ESRRG, SLC34A1, and EFHD1) were included in LASSO analysis (Figures 4B, C). Following the cross validation, seven genes achieved the minimum partial likelihood deviance. Then, we performed a multivariate Cox regression with these seven genes (GGT6, FAM3B, PTH1R, ALDOB, ESRRG, SLC34A1, and EFHD1) as covariants. We finally got three genes, including ALDOB, ESRRG, and EFHD1 without collinearity, and each of them could be an independent prognostic marker for ccRCC (Figure 4D). A prognostic model based on the three genes was established. The risk score for each individual patient was calculated with the following formula: risk score = (-0.105197) * ALDOB + (-0.275676) * ESRRG + (-0.269554) * EFHD1.

TABLE 1
www.frontiersin.org

Table 1 Univariate Cox regression analysis in the train cohort.

Then, the Kaplan–Meier analysis was performed. As shown in Figure 4E, the survival rate of patients in the high-risk group was significantly lower than that in the low-risk group in either training cohort (p < 5.81e-4) or testing cohort (p < 5.48e-20). The ROC curve was then used to evaluate the accuracy of the survival analysis. The areas under the curves (AUCs) were 0.717 and 0.699 in the training cohort and testing cohort, respectively (Figure 4E), which indicate that the prediction effect was good. We also plotted the distribution of risk scores in patients with ccRCC and the correlation between survival time and risk scores in the training cohort and testing cohort (Figure 4F). In addition, all of the three genes (ALDOB, ESRRG, and EFHD1) were significantly downregulated and associated with poor pathologic stages in the training cohort, testing cohort, or GSE126964 dataset (Figure 4F and Figure S4). Moreover, ALDOB, ESRRG, and EFHD1 were highly expressed in renal tissues among the different normal tissues, which indicated a critical regulatory function of these genes in the normal kidney (Figure S5).

Validation of Protein Expressions of Prognostic Genes

Immunohistochemistry staining results obtained from the HPA database revealed the protein expression levels of the key survival-related genes (Figure 5). The results showed the downregulation of ALDOB, EFHD1, and ESRRG protein in ccRCC samples compared with normal controls.

FIGURE 5
www.frontiersin.org

Figure 5 Immunohistochemistry staining of prognostic proteins based on the HPA. Protein expression levels of (A) ALDOB, (B) EFHD1, and (C) ESRRG in tumor tissues and normal tissues. HPA, Human Protein Atlas; Tub., tubules; Glo., glomeruli.

Gene Set Enrichment Analysis of Prognostic Genes

GSEA was conducted to search the KEGG pathways in which the prognostic genes and risk scores were enriched in the samples with high expression or high-risk levels from TCGA-KIRC. “Oxidative phosphorylation” and “Fatty acid metabolism” pathways were enriched with low-risk score and high expression of ALDOB, ESRRG, and EFHD1, while immune-related pathways, including “Cytokine–cytokine receptor interaction,” “Chemokine signaling pathway,” and “Primary immunodeficiency” were significantly enriched with high-risk score and low expression of the prognostic genes (Figures 6A–D).

FIGURE 6
www.frontiersin.org

Figure 6 Single-gene GSEA of (A) ALDOB, (B) EFHD1, (C) ESRRG and (D) risk score based on TCGA-KIRC. GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; TCGA-KIRC, The Cancer Genome Atlas Kidney Renal Clear Cell Carcinoma.

Single-Cell Transcriptomic Context of the Prognostic Genes

To further verify the relationship among ALDOB, ESRRG, and EFHD1 in ccRCC, single-cell RNA-seq data from GSE131685 and GSE171306 were employed (25, 26). After quality control, a total of 34,371 cells from two ccRCC samples and three normal kidney samples were profiled (Figure 7A). We identified 27 different cell clusters and five cell groups, including immune cells, epithelial cells, endothelial cells, mesenchymal cells, and tumor cells (Figures 7B, C). Consistent with previous research (28), proximal tubular epithelial cells account for over 90% of a normal renal cortical sample, while in ccRCC, over 50% was accounted for immune cells and approximate 20% for tumor cells (Figure 7D). Except the clusters of macrophage 1 (MC1) and T cell 2 (T2), most kinds of immune cells were identified from ccRCC patients, which depicted a tumor immune microenvironment of ccRCC (Figure 7E). We also identified four tumor cell clusters. Analysis of KEGG pathway in tumor cells suggested the increased glycolysis gluconeogenesis, cancer, and focal adhesion-associated metabolism in ccRCC, while oxidative phosphorylation-associated pathways were negatively enriched with tumor cells (Figures S6, S7).

FIGURE 7
www.frontiersin.org

Figure 7 Prognostic expression profile based on single-cell sequencing analysis. (A) Composition and distribution of single cells from GSE131685 and GSE171306. (B) UMAP embedding of 34,371 single cells from three human normal kidneys and two ccRCC samples. Labels refer to 27 clusters identified. (C) Scaled gene expression of the top 10 specific genes in each cluster. Each column is the average expression of all cells in a cluster. (D) Composition of different cell types in the five single-cell RNA-seq samples. (E) Number of cells per immune cell type and clinical parameter. The expression profile of (F) ALDOB, (G) EFHD1, and (H) ESRRG for each cell and the (I) violin plot. UMAP, Uniform Manifold Approximation and Projection; ccRCC, clear cell renal cell carcinoma; Im., immune; Epi., epithelial; Endo., endothelial; Mes, mesenchymal; CD, collecting duct; CT, connecting tubule; iEn, injured endothelial cells; Fib, fibroblast; Mast, mast cell; MC, macrophage; Mono, monocyte; PT, proximal tubule; iPT, injured proximal tubule; VR, vasa recta.

We then explored the expression profile of ALDOB, ESRRG, and EFHD1 in different types of cells. Similar with results from TCGA-KIRC and GSE126964, the expression of ALDOB, ESRRG, and EFHD1 were much lower in tumor cells than that in other intrinsic renal cells (Figures 7F–I).

Discussion

ccRCC is a common genitourinary cancer with a high mortality rate (3). There is an urgent requirement to identify additional potential targets for drugs and biomarkers for early diagnosis of ccRCC. In the present study, a novel prognostic model based on three genes (ALDOB, EFHD1, and ESRRG) for ccRCC was established. ALDOB, EFHD1, and ESRRG were also identified as novel independent prognostic markers for ccRCC in different datasets via integrated bioinformatics analysis, including DEG analysis, WGCNA, and single-cell analysis.

Metabolic Shift in Clear Cell Renal Cell Carcinoma

Metabolic disorder is a hallmark in different types of cancer, since sufficient energy and metabolite production are required for the malignant proliferation of cancer cells (34). Gebhard et al. (35) reported that ccRCC tissues were overloaded with glycogen and lipid compared with normal tissues, which suggested that the metabolism of lipids and glucose may be altered in ccRCC (36). In particular, a mutation in VHL is considered to be closely associated with metabolic reprogramming in ccRCC (37). Subsequent accumulation of HIF1 α leads to the expression of glucose transporter-1, thereby promoting cellular glucose uptake. In addition, it can induce lactate dehydrogenase, which promotes the conversion of pyruvate to lactate and switches energy production from the TCA to lactate fermentation (38). This phenomenon is widely known as the Warburg effect. Despite the well-known VHL–HIF axis, there are a number of altered levels of the biochemical enzymes, substrates, and metabolic intermediates or products that are involved in the metabolic reprogramming waiting to be discovered. Our analysis of the GSE126964 data and single-cell RNA-seq data from GSE131685 and GSE171306 supported the observations of the metabolic shift in ccRCC and demonstrated the upregulation of glycolysis-associated genes and downregulation of OXPHOS and TCA-associated genes at the transcriptome level (Figure 6 and Figures S2, S6, S7). Although a metabolic shift is advantageous for tumor progression, altered metabolic pathways in ccRCC may also be exploited as therapeutic targets and may be an important future research direction.

ALDOB, EFHD1, and ESRRG Can Be Novel Independent Prognostic Markers for Clear Cell Renal Cell Carcinoma

ALDOB encodes aldolase B, an enzyme that is expressed in the liver and kidneys and is involved in glycolysis process and fructolysis process. The function of which can cleave fructose-1,6-bisphosphate to yield glyceraldehyde and dihydroxyacetone phosphate (39). A research found that declined ALDOB expression was associated with multiple malignant characteristics of HCC and indicate a poor prognosis (40). Moreover, Bu et al. (41) shows that ALDOB upregulation is commonly found in the metastatic cell in liver during primary colon cancer proliferation by enhancing fructose metabolism and central carbon metabolism. A study by Wang et al. (42) found that the low expression of ALDOB is also important in ccRCC and predicts poor prognosis, which is consistent with our research results. It leads to a high level of fructose 1,6-bisphosphate (FBP) and protects ccRCC from oxidative stress (42). However, the mechanism and prognostic value of accumulated FBP remain unknown.

EFHD1 encodes a mitochondrial inner membrane protein, acted as a calcium sensor for mitochondrial flash activation (43), and induced metabolic changes during the development of pro-/pre-B cells (44). A recent report suggested that EFHD1 may interact with β-actin for its involvement in the Ca2+-dependent regulation of mitochondrial morphology (45). EFHD1 was significantly downregulated in both the GSE126964 and TCGA-KIRC cohorts (Figure 4) and may also have an impact on mitochondrial energy metabolism in ccRCC. However, the detailed mechanism of how EFHD1 regulates ccRCC pathogenesis is currently unknown.

ESRRG encodes a member of nuclear receptor superfamily of transcription factors and has been shown to be a tumor suppressor in different types of cancer (4648). A study by Huang et al. (49) also identified ESRRG as a co-expressed DEG in different datasets of hypertension-related RCC. Moreover, an experimental study on the mechanism of ESRRG conducted by Nam et al. (50) demonstrated that ESRRG suppressed the migratory and invasive abilities of behaviors in RCC cells. Our analysis and previous studies have all suggested that lower ESRRG expression may be a reliable predictor of a poor clinical outcome.

Hub Genes Including ARMH4, PTH1R, SLC22A8, and SLC34A1 Were Correlated With Cancer

ARMH4, also named C14ORF37, encodes a protein that contains an armadillo-like helical domain. It has been shown that ARMH4 can interact with and inhibit the function of mTOR complex 2 kinase activity and function as a tumor suppressor in hematological malignancies, which is driven by interleukin 6 (IL-6)–signal transducer and activator of transcription 3 (STAT3) signaling pathways (51, 52). Wang et al. (53) also predicted that ARMH4 may act as a modulator for QKI, KH domain containing RNA binding, one of the key RNA-binding proteins shown in TCGA-KIRC dataset, and may change its splicing regulation in kidney cancer (53). Our analysis further confirmed the importance of ARMH4 in ccRCC. However, the specific mechanism remains to be explored.

PTH1R encodes a G protein-coupled receptor of parathyroid hormone (PTH) and PTH-related protein and plays a central role in calcium homeostasis (54). Recently, the structure and dynamics of the active PTH1R have been shown by cryo-electron microscopy (55). Studies from different research groups have reported that the decreased expression of PTH1R was a poor prognosis factor in multiple types of cancer (5659). Although PTH1R was found to be highly expressed in normal kidney samples (Figure 5 and Figure S4), the detailed mechanism of PTH1R in renal function and ccRCC has yet to be fully elucidated and requires further study.

SLCs are a superfamily of membrane proteins responsible for the cellular uptake of a diverse range of substances. Among the SLCs, SLC22A8 and SLC34A1 both show kidney-specific expression. SLC22A8 is involved in the sodium-independent transport and excretion of organic anions, while SLC34A1 is a sodium-phosphate co-transporter that controls proximal tubule phosphate reabsorption (60, 61). Although defects in SLCs can lead to serious diseases (6264), there is lack of research in cancer, especially ccRCC. Kang et al. (65) examined the expression patterns and prognostic values of SLCs in the development of ccRCC using different bioinformatics methods. These authors demonstrated that the low expression levels of a cluster of SLCs, including SLC34A1, were correlated with ccRCC progression and poor prognosis (65). Since ccRCC shows a prominent metabolic shift effect, the production and accumulation of metabolites are also different from those in normal tissues. Therefore, SLCs may be critical in ccRCC.

Conclusion

Through a series of comprehensive bioinformatics analyses, including DEG screening, WGCNA, and single-cell analysis, a prognostic model based on ALDOB, EFHD1, and ESRRG was established, and these three genes were also identified as independent prognostic markers for ccRCC. The aforementioned prognostic genes have the potential to become therapeutic targets and biomarkers for ccRCC. However, these key survival-related genes should be tested in a large cohort of ccRCC cases and should be analyzed and validated in additional in vivo and in vitro experiments.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics Statement

The review board of the Xiangya Hospital of Central South University approved the present study.

Author Contributions

ZP, RX, and LT conceived and directed the project. HH and LZ collected the data and information. HH, CH, YD, and LF analyzed and interpreted the data. HH wrote the manuscript with the help of all the other authors. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the project funded by China Postdoctoral Science Foundation (2020TQ0363 and 2020M682598), the National Natural Science Foundation of China (82073918, 82090020, and 82090024), the Natural Science Foundation of Hunan, China (2021JJ40992), and the Fundamental Research Funds for the Central Universities of Central South University (2021zzts0078).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

Acknowledgments

We would like to acknowledge the reviewers for their helpful comments on this paper.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.726655/full#supplementary-material

Supplementary Figure 1 | PCA based on the whole gene list. PCA, principal component analysis.

Supplementary Figure 2 | Schematic of metabolic pathways and their selected genes for glycolysis, TCA and electron transport chain, with the log2FC level between ccRCC and control samples in the GSE126964 dataset. *p<0.05; **p<0.01; ***p<0.001. FC, fold change; ccRCC, clear cell renal cell carcinoma; TCA, tricarboxylic acid cycle.

Supplementary Figure 3 | (A) Sample clustering was conducted to detect outliers. All samples are located in the clusters and pass the cutoff thresholds after removing outliers. (B) Heatmap depicting the TOM of genes selected for weighted co-expression network analysis. The light color represents lower overlap, and red represents higher overlap. TOM, Topological Overlap Matrix.

Supplementary Figure 4 | The overall survival stratified by the high and low-risk score groups was plotted for the GSE126964 dataset. Detailed risk scores, survival information and heat maps of gene expression are also included.

Supplementary Figure 5 | Prognostic genes expression profiles in different normal tissues.

Supplementary Figure 6 | Single-gene GSEA of four tumor clusters for KEGG pathway enrichment. GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; FC, fold change.

Supplementary Figure 7 | The expression profile of (A) OXPHOS-associated and (B) glycolysis-associated genes for each cell. OXPHOS, oxidative phosphorylation.

Supplementary Figure 8 | The expression profile of other key genes for each cell.

References

1. Jonasch E, Gao J, Rathmell WK. Renal Cell Carcinoma. BMJ (2014) 349:g4797. doi: 10.1136/bmj.g4797

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2021) 71:209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Hsieh JJ, Purdue MP, Signoretti S, Swanton C, Albiges L, Schmidinger M, et al. Renal Cell Carcinoma. Nat Rev Dis Primers (2017) 3:17009. doi: 10.1038/nrdp.2017.9

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Angulo JC, Manini C, Lopez JI, Pueyo A, Colas B, Ropero S. The Role of Epigenetics in the Progression of Clear Cell Renal Cell Carcinoma and the Basis for Future Epigenetic Treatments. Cancers (Basel) (2021) 13:2071. doi: 10.3390/cancers13092071

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kotecha RR, Motzer RJ, Voss MH. Towards Individualized Therapy for Metastatic Renal Cell Carcinoma. Nat Rev Clin Oncol (2019) 16:621–33. doi: 10.1038/s41571-019-0209-1

PubMed Abstract | CrossRef Full Text | Google Scholar

6. The Cancer Genome Atlas Research Network. Comprehensive Molecular Characterization of Clear Cell Renal Cell Carcinoma. Nature (2013) 499:43–9. doi: 10.1038/nature12222

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Latif F, Tory K, Gnarra J, Yao M, Duh FM, Orcutt ML, et al. Identification of the Von Hippel-Lindau Disease Tumor Suppressor Gene. Science (1993) 260:1317–20. doi: 10.1126/science.8493574

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Yao X, Tan J, Lim KJ, Koh J, Ooi WF, Li Z, et al. VHL Deficiency Drives Enhancer Activation of Oncogenes in Clear Cell Renal Cell Carcinoma. Cancer Discov (2017) 7:1284–305. doi: 10.1158/2159-8290.CD-17-0375

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Zhang J, Wu T, Simon J, Takada M, Saito R, Fan C, et al. VHL Substrate Transcription Factor ZHX2 as an Oncogenic Driver in Clear Cell Renal Cell Carcinoma. Science (2018) 361:290–5. doi: 10.1126/science.aap8411

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Choueiri TK, Kaelin WJ. Targeting the HIF2-VEGF Axis in Renal Cell Carcinoma. Nat Med (2020) 26:1519–30. doi: 10.1038/s41591-020-1093-z

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zhang J, Zhang Q. VHL and Hypoxia Signaling: Beyond HIF in Cancer. Biomedicines (2018) 6:35. doi: 10.3390/biomedicines6010035

CrossRef Full Text | Google Scholar

12. Majmundar AJ, Wong WJ, Simon MC. Hypoxia-Inducible Factors and the Response to Hypoxic Stress. Mol Cell (2010) 40:294–309. doi: 10.1016/j.molcel.2010.09.022

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Krieg M, Haas R, Brauch H, Acker T, Flamme I, Plate KH. Up-Regulation of Hypoxia-Inducible Factors HIF-1alpha and HIF-2alpha Under Normoxic Conditions in Renal Carcinoma Cells by Von Hippel-Lindau Tumor Suppressor Gene Loss of Function. Oncogene (2000) 19:5435–43. doi: 10.1038/sj.onc.1203938

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Ferraro E, Germano M, Mollace R, Mollace V, Malara N. HIF-1, the Warburg Effect, and Macrophage/Microglia Polarization Potential Role in COVID-19 Pathogenesis. Oxid Med Cell Longev (2021) 2021:8841911. doi: 10.1155/2021/8841911

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Brugarolas J. PBRM1 and BAP1 as Novel Targets for Renal Cell Carcinoma. Cancer J (2013) 19:324–32. doi: 10.1097/PPO.0b013e3182a102d1

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Hakimi AA, Ostrovnaya I, Reva B, Schultz N, Chen YB, Gonen M, et al. Adverse Outcomes in Clear Cell Renal Cell Carcinoma With Mutations of 3p21 Epigenetic Regulators BAP1 and SETD2: A Report by MSKCC and the KIRC TCGA Research Network. Clin Cancer Res (2013) 19:3259–67. doi: 10.1158/1078-0432.CCR-12-3886

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Krishna C, DiNatale RG, Kuo F, Srivastava RM, Vuong L, Chowell D, et al. Single-Cell Sequencing Links Multiregional Immune Landscapes and Tissue-Resident T Cells in ccRCC to Tumor Topology and Therapy Efficacy. Cancer Cell (2021) 39:662–77. doi: 10.1016/j.ccell.2021.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Du W, Zhang L, Brett-Morris A, Aguila B, Kerner J, Hoppel CL, et al. HIF Drives Lipid Deposition and Cancer in ccRCC via Repression of Fatty Acid Metabolism. Nat Commun (2017) 8:1769. doi: 10.1038/s41467-017-01965-8

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Clark DJ, Dhanasekaran SM, Petralia F, Pan J, Song X, Hu Y, et al. Integrated Proteogenomic Characterization of Clear Cell Renal Cell Carcinoma. Cell (2020) 180:207. doi: 10.1016/j.cell.2019.12.026

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Langfelder P, Horvath S. WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

21. Yu J, Mao W, Sun S, Hu Q, Wang C, Xu Z, et al. Identification of an M6a-Related lncRNA Signature for Predicting the Prognosis in Patients With Kidney Renal Clear Cell Carcinoma. Front Oncol (2021) 11:663263. doi: 10.3389/fonc.2021.663263

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhou Z, Bai J, Zhong S, Zhang R, Kang K, Zhang X, et al. Downregulation of ATP6V1A Involved in Alzheimer's Disease via Synaptic Vesicle Cycle, Phagosome, and Oxidative Phosphorylation. Oxid Med Cell Longev (2021) 2021:5555634. doi: 10.1155/2021/5555634

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zhao Q, Xue J, Hong B, Qian W, Liu T, Fan B, et al. Transcriptomic Characterization and Innovative Molecular Classification of Clear Cell Renal Cell Carcinoma in the Chinese Population. Cancer Cell Int (2020) 20:461. doi: 10.1186/s12935-020-01552-w

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Hao Y, Hao S, Andersen-Nissen E, Mauck WR, Zheng S, Butler A, et al. Integrated Analysis of Multimodal Single-Cell Data. Cell (2021) 184:3573–87. doi: 10.1016/j.cell.2021.04.048

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Liao J, Yu Z, Chen Y, Bao M, Zou C, Zhang H, et al. Single-Cell RNA Sequencing of Human Kidney. Sci Data (2020) 7:4. doi: 10.1038/s41597-019-0351-8

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yu Z, Lu W, Su C, Lv Y, Ye Y, Guo B, et al. Single-Cell RNA-Seq Identification of the Cellular Molecular Characteristics of Sporadic Bilateral Clear Cell Renal Cell Carcinoma. Front Oncol (2021) 11:659251. doi: 10.3389/fonc.2021.659251

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Wang Z, Su G, Dai Z, Meng M, Zhang H, Fan F, et al. Circadian Clock Genes Promote Glioma Progression by Affecting Tumour Immune Infiltration and Tumour Cell Proliferation. Cell Prolif (2021) 54:e12988. doi: 10.1111/cpr.12988

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Kuppe C, Ibrahim MM, Kranz J, Zhang X, Ziegler S, Perales-Paton J, et al. Decoding Myofibroblast Origins in Human Kidney Fibrosis. Nature (2021) 589:281–6. doi: 10.1038/s41586-020-2941-1

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Borcherding N, Vishwakarma A, Voigt AP, Bellizzi A, Kaplan J, Nepple K, et al. Mapping the Immune Environment in Clear Cell Renal Carcinoma by Single-Cell Genomics. Commun Biol (2021) 4:122. doi: 10.1038/s42003-020-01625-6

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: A Web Server for Cancer and Normal Gene Expression Profiling and Interactive Analyses. Nucleic Acids Res (2017) 45:W98–102. doi: 10.1093/nar/gkx247

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yu G, Wang LG, Han Y, He QY. ClusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. Omics (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Otasek D, Morris JH, Boucas J, Pico AR, Demchak B. Cytoscape Automation: Empowering Workflow-Based Network Analysis. Genome Biol (2019) 20:185. doi: 10.1186/s13059-019-1758-4

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Coller HA. Is Cancer a Metabolic Disease? Am J Pathol (2014) 184:4–17. doi: 10.1016/j.ajpath.2013.07.035

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Gebhard RL, Clayman RV, Prigge WF, Figenshau R, Staley NA, Reesey C, et al. Abnormal Cholesterol Metabolism in Renal Clear Cell Carcinoma. J Lipid Res (1987) 28:1177–84. doi: 10.1016/S0022-2275(20)38606-5

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Hakimi AA, Reznik E, Lee CH, Creighton CJ, Brannon AR, Luna A, et al. An Integrated Metabolic Atlas of Clear Cell Renal Cell Carcinoma. Cancer Cell (2016) 29:104–16. doi: 10.1016/j.ccell.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Nickerson ML, Jaeger E, Shi Y, Durocher JA, Mahurkar S, Zaridze D, et al. Improved Identification of Von Hippel-Lindau Gene Alterations in Clear Cell Renal Tumors. Clin Cancer Res (2008) 14:4726–34. doi: 10.1158/1078-0432.CCR-07-4921

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Weidemann A, Johnson RS. Biology of HIF-1alpha. Cell Death Differ (2008) 15:621–7. doi: 10.1038/cdd.2008.12

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Chang YC, Yang YC, Tien CP, Yang CJ, Hsiao M. Roles of Aldolase Family Genes in Human Cancers and Diseases. Trends Endocrinol Metab (2018) 29:549–59. doi: 10.1016/j.tem.2018.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Tao QF, Yuan SX, Yang F, Yang S, Yang Y, Yuan JH, et al. Aldolase B Inhibits Metastasis Through Ten-Eleven Translocation 1 and Serves as a Prognostic Biomarker in Hepatocellular Carcinoma. Mol Cancer (2015) 14:170. doi: 10.1186/s12943-015-0437-7

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Bu P, Chen KY, Xiang K, Johnson C, Crown SB, Rakhilin N, et al. Aldolase B-Mediated Fructose Metabolism Drives Metabolic Reprogramming of Colon Cancer Liver Metastasis. Cell Metab (2018) 27:1249–62. doi: 10.1016/j.cmet.2018.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Wang J, Wu Q, Qiu J. Accumulation of Fructose 1,6-Bisphosphate Protects Clear Cell Renal Cell Carcinoma From Oxidative Stress. Lab Invest (2019) 99:898–908. doi: 10.1038/s41374-019-0203-3

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Hou T, Jian C, Xu J, Huang AY, Xi J, Hu K, et al. Identification of EFHD1 as a Novel Ca(2+) Sensor for Mitoflash Activation. Cell Calcium (2016) 59:262–70. doi: 10.1016/j.ceca.2016.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Stein M, Dutting S, Mougiakakos D, Bosl M, Fritsch K, Reimer D, et al. A Defined Metabolic State in Pre B Cells Governs B-Cell Development and is Counterbalanced by Swiprosin-2/Efhd1. Cell Death Differ (2017) 24:1239–52. doi: 10.1038/cdd.2017.52

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Mun SA, Park J, Park KR, Lee Y, Kang JY, Park T, et al. Structural and Biochemical Characterization of EFhd1/Swiprosin-2, an Actin-Binding Protein in Mitochondria. Front Cell Dev Biol (2020) 8:628222. doi: 10.3389/fcell.2020.628222

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Kang MH, Choi H, Oshima M, Cheong JH, Kim S, Lee JH, et al. Estrogen-Related Receptor Gamma Functions as a Tumor Suppressor in Gastric Cancer. Nat Commun (2018) 9:1920. doi: 10.1038/s41467-018-04244-2

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Tiraby C, Hazen BC, Gantner ML, Kralli A. Estrogen-Related Receptor Gamma Promotes Mesenchymal-to-Epithelial Transition and Suppresses Breast Tumor Growth. Cancer Res (2011) 71:2518–28. doi: 10.1158/0008-5472.CAN-10-1315

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Kang MH, Eyun SI, Park YY. Estrogen-Related Receptor-Gamma Influences Helicobacter Pylori Infection by Regulating TFF1 in Gastric Cancer. Biochem Biophys Res Commun (2021) 563:15–22. doi: 10.1016/j.bbrc.2021.05.076

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Huang W, Wu K, Wu R, Chen Z, Zhai W, Zheng J. Bioinformatic Gene Analysis for Possible Biomarkers and Therapeutic Targets of Hypertension-Related Renal Cell Carcinoma. Transl Androl Urol (2020) 9:2675–87. doi: 10.21037/tau-20-817

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Nam HY, Chandrashekar DS, Kundu A, Shelar S, Kho EY, Sonpavde G, et al. Integrative Epigenetic and Gene Expression Analysis of Renal Tumor Progression to Metastasis. Mol Cancer Res (2019) 17:84–96. doi: 10.1158/1541-7786.MCR-17-0636

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Conway SJ, McConnell R, Simmons O, Snider PL. Armadillo-Like Helical Domain Containing-4 is Dynamically Expressed in Both the First and Second Heart Fields. Gene Expr Patterns (2019) 34:119077. doi: 10.1016/j.gep.2019.119077

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Lee D, Wang YH, Kalaitzidis D, Ramachandran J, Eda H, Sykes DB, et al. Endogenous Transmembrane Protein UT2 Inhibits Pstat3 and Suppresses Hematological Malignancy. J Clin Invest (2016) 126:1300–10. doi: 10.1172/JCI84620

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Wang Y, Chen SX, Rao X, Liu Y. Modulator-Dependent RBPs Changes Alternative Splicing Outcomes in Kidney Cancer. Front Genet (2020) 11:265. doi: 10.3389/fgene.2020.00265

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Bastepe M, Turan S, He Q. Heterotrimeric G Proteins in the Control of Parathyroid Hormone Actions. J Mol Endocrinol (2017) 58:R203–24. doi: 10.1530/JME-16-0221

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Zhao LH, Ma S, Sutkeviciute I, Shen DD, Zhou XE, de Waal PW, et al. Structure and Dynamics of the Active Human Parathyroid Hormone Receptor-1. Science (2019) 364:148–53. doi: 10.1126/science.aav7942

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Wang HJ, Wang L, Song SS, He XL, Pan HY, Hu ZM, et al. Decreased Expression of PTH1R is a Poor Prognosis in Hepatocellular Carcinoma. Cancer Biomark (2018) 21:723–30. doi: 10.3233/CBM-170823

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Li S, Chen P, Zheng K, Wang W, Pei Y, Qiu E, et al. Beta-Alanine Mediated Inhibition of PTHR1suppresses the Proliferation, Invasion and Tumorigenesis in Metastatic Human Osteosarcoma U2OS Cells. Int J Biol Macromol (2018) 111:1255–63. doi: 10.1016/j.ijbiomac.2018.01.106

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Ho PW, Goradia A, Russell MR, Chalk AM, Milley KM, Baker EK, et al. Knockdown of PTHR1 in Osteosarcoma Cells Decreases Invasion and Growth and Increases Tumor Differentiation In Vivo. Oncogene (2015) 34:2922–33. doi: 10.1038/onc.2014.217

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Garcia M, Rodriguez-Hernandez CJ, Mateo-Lozano S, Perez-Jaume S, Goncalves-Alves E, Lavarino C, et al. Parathyroid Hormone-Like Hormone Plays a Dual Role in Neuroblastoma Depending on PTH1R Expression. Mol Oncol (2019) 13:1959–75. doi: 10.1002/1878-0261.12542

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Bush KT, Singh P, Nigam SK. Gut-Derived Uremic Toxin Handling In Vivo Requires OAT-Mediated Tubular Secretion in Chronic Kidney Disease. JCI Insight (2020) 5:e133817. doi: 10.1172/jci.insight.133817

CrossRef Full Text | Google Scholar

61. Janiec A, Halat-Wolska P, Obrycki L, Ciara E, Wojcik M, Pludowski P, et al. Long-Term Outcome of the Survivors of Infantile Hypercalcaemia With CYP24A1 and SLC34A1 Mutations. Nephrol Dial Transplant (2020) 36:1484–92. doi: 10.1093/ndt/gfaa178

CrossRef Full Text | Google Scholar

62. McNeill A, Iovino E, Mansard L, Vache C, Baux D, Bedoukian E, et al. SLC12A2 Variants Cause a Neurodevelopmental Disorder or Cochleovestibular Defect. Brain (2020) 143:2380–7. doi: 10.1093/brain/awaa176

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Preising MN, Gorg B, Friedburg C, Qvartskhava N, Budde BS, Bonus M, et al. Biallelic Mutation of Human SLC6A6 Encoding the Taurine Transporter TAUT is Linked to Early Retinal Degeneration. FASEB J (2019) 33:11507–27. doi: 10.1096/fj.201900914RR

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Hu C, Tao L, Cao X, Chen L. The Solute Carrier Transporters and the Brain: Physiological and Pharmacological Implications. Asian J Pharm Sci (2020) 15:131–44. doi: 10.1016/j.ajps.2019.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Kang W, Zhang M, Wang Q, Gu D, Huang Z, Wang H, et al. The SLC Family Are Candidate Diagnostic and Prognostic Biomarkers in Clear Cell Renal Cell Carcinoma. BioMed Res Int (2020) 2020:1932948. doi: 10.1155/2020/1932948

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: clear cell renal cell carcinoma, weighted gene co-expression network analysis, differentially expressed genes, prognostic genes, single-cell analysis

Citation: Huang H, Zhu L, Huang C, Dong Y, Fan L, Tao L, Peng Z and Xiang R (2021) Identification of Hub Genes Associated With Clear Cell Renal Cell Carcinoma by Integrated Bioinformatics Analysis. Front. Oncol. 11:726655. doi: 10.3389/fonc.2021.726655

Received: 17 June 2021; Accepted: 06 September 2021;
Published: 30 September 2021.

Edited by:

Marco Borghesi, University of Genoa, Italy

Reviewed by:

Shaolong Cao, University of Texas MD Anderson Cancer Center, United States
Pasquale Ditonno, University of Bari, Italy

Copyright © 2021 Huang, Zhu, Huang, Dong, Fan, Tao, Peng and Xiang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhangzhe Peng, pengzhangzhe@csu.edu.cn; Rong Xiang, shirlesmile@csu.edu.cn

Disclaimer: 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.