RNA sequencing data of 72 normal renal samples and 539 ccRCC samples with corresponding somatic mutation data and clinical data were obtained from The Cancer Genome Atlas (TCGA) database (https://cancerg enome.nih.gov/). Microarray data sets (GSE46699 (63 normal renal and 67 ccRCC samples), GSE36895 (23 normal renal and 29 ccRCC samples), and GSE16449 (18 normal renal and 52 ccRCC samples)) and clinical data (GSE3, GSE22541 and GSE29609) were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/).

Clinical and pathological data of 59 patients undergoing renal cancer resection in Beijing Hospital were collected retrospectively and further screening and analysis were performed. The inclusion criteria were as follows: (Siegel et al., 2020): patients who underwent nephrectomy and pathological diagnosis of ccRCC in Beijing Hospital; (Bray et al., 2018); surgery occurred between October 2019 and June 2021; and (Bhatt and Finelli, 2014) patients were ≥18 years old. The exclusion criterion was: (Siegel et al., 2020) no fresh tissue samples available. All fresh tissue samples were collected immediately after surgical resection, quickly frozen in liquid nitrogen and stored at −80 °C until RNA was extracted from the samples. Reverse transcription polymerase chain reaction (RT-PCR) and western blotting (WB) were performed on the fresh ccRCC samples to assess CGN expression. All protocols were approved by the Research Ethics Committee of Beijing Hospital.

To clarify the association between genes, a weighted gene co-expression network was constructed for the hierarchical clustering of all genes subjected to co-expression analysis. The protocol used was as follows: (Siegel et al., 2020): The Pearson correlation coefficient between genes was determined. (Bray et al., 2018). A weighted adjacency matrix α _mn = |c_mn|^β was constructed, where α _mn represents the adjacency matrix between gene m and gene n; c_mn represents the Pearson correlation coefficient between gene m and gene n; and β is the soft power value. The pickSoftThreshold function of the WGCNA software package was used to select the appropriate soft threshold power β. An appropriate soft power value β can ensure that the network is in accordance with a standard non-scale network to achieve a scale-free topology. (Bhatt and Finelli, 2014). The adjacency matrix was converted into a topological overlap matrix (TOM), and the dissimilarity matrix between genes dissTOM = 1-TOM was calculated. (Linehan, 2012). Hierarchical clustering was performed using dissTOM such that genes with similar expression patterns are placed in the same gene module. (Lee et al., 2002). The minimum number of module genes was set at 50, and the dynamic hybrid cutting algorithm was used to obtain the gene modules and merge modules that were highly similar. (Patard et al., 2003). The Pearson correlation coefficient between each module and the disease traits was determined, and the p value was used to determine the hub module. The genes in the gene module with the highest association coefficient were used as candidate prognostic molecular markers and were included in subsequent analyses.

The DEGs were screened from the normal and tumor tissue groups. The screening criteria were false discovery rate (FDR) < 0.05 and │log₂ FC│≥1.0, where FC is the fold change that is the multiple of the differential expression levels between the two groups. Subsequently, the genes of the hub module and the DEGs in the weighted co-expression network were intersected to identify the differentially co-expressed genes, which were visualized using the R package VennDiagram. A hypergeometric test was used to test the statistical significance of Venn diagram, and p < 0.05 was considered statistically significant.

The clusterProfiler package in the R software was used to perform a GO enrichment analysis on the intersecting genes. A recently developed algorithm based on the pathway topology to identify the functional roles of pathway components was used for pathway annotations (Sorokin et al., 2021).

R language survival package was used to perform a survival analysis on the intersecting genes. From the TCGA database, univariate Cox survival analysis-screened genes with a FDR<0.05 were included in the multivariate Cox survival analysis to obtain independent prognostic genes. A Kaplan–Meier (KM) survival analysis was then performed to screen for the most relevant genes for prognosis. Finally, TCGA, GSE3, GSE22541 and GSE29609 datasets were integrated with the application of the “ComBat” algorithm to eliminate batch effect, and the prognostic effect of CGN was verified.

To explore the relationship between CGN and clinical traits, CGN expression level was divided into two groups according to the median value. The Wilcox test was used to explore the relationship between CGN and age, sex, T stage, N stage, M stage, tumor grade, and tumor stage. Statistical significance was set at p < 0.05.

The “CIBERSORT” R package was used to analyze the infiltrating immune components of each sample, and the Pearson correlation analysis was used to determine the linear relationship between the 22 immune cell types (naïve B cells, memory B cells, plasma cells, CD8 T cells, CD4 naïve T cells, CD4 memory resting T cells, CD4 memory activated T cells, follicular helper T cells, regulatory T cells, gamma delta T cells, resting NK cells, activated NK cells, Monocytes, M0 macrophages, M1 macrophages, M2 macrophages, resting dendritic cells, activated dendritic cells, resting mast cells, activated mast cells, eosinophils, and neutrophils) and CGN. We further divided the expression of CGN into high and low expression groups according to the median value, and the Wilcoxon test was used to detect the immune cells related to CGN. Finally, the intersection of the two results identified the immune cells that are most closely related to CGN.

The linear relationship between CGN and TMB in the TCGA dataset was obtained by the Pearson correlation analysis. Then, CGN was divided into high and low expression groups based on the median value. The Wilcoxon test was used to compare the relationship between these two groups and TMB. TMB was divided into two groups based on the median value and a KM survival analysis was performed. The KM method and log-rank test shows the survival curves of CGN stratification and TMB stratification. The Wilcoxon test was used to compare the relationship between high and low expressions CGN and PD-L1 of TCGA and GSE16449.

The linear relationship between CGN and existing ccRCC biomarkers (BAP1, BIRC5, CXCR4, and SETD2) in the TCGA dataset was obtained by the Pearson correlation analysis. Statistical significance was set at p < 0.05.

Protein expression of the CGN gene in the tumor and normal tissues was obtained using the HPA database (https://www.proteinatlas.org/).

The mRNA transcription level of CGN in the tumor and normal tissues of the ccRCC patients that were preserved in our hospital was analyzed. TRI Reagent (Sigma) was used to extract the total RNA from the specimens. A NanoDrop^® ND-1000 spectrophotometer was used to determine the RNA concentration and purity. Total RNA was reverse-transcribed into cDNA using SuperScriptTM III Reverse transcriptase (Invitrogen), and the cDNA was used as a template to detect the expression of each gene by RT-PCR. The primer sequences used were as follows:

CGN forward: 5′-CAG​GGC​ATT​GGC​AGA​GTA​TGT-3′;

CGN reverse: 5′-CCT​CAA​CCT​GGC​GAG​TAT​CT-3′; β-actin forward: 5′-GTG​GCC​GAG​GAC​TTT​GAT​TG-3′;

β-actin reverse: 5′-CCT​GTA​ACA​ACG​CAT​CTC​ATA​TT-3′.

Total protein was extracted with pre-chilled RIPA lysis buffer and was quantified using a BCA protein assay kit (Cwbiotech), according to the manufacturer’s guidelines. SDS-PAGE was performed on the protein samples and blocked at room temperature after transfer to PVDF (Millipore) membranes with added anti-CGN (Sigma) primary antibody diluted 1:1000 and incubated overnight. A secondary antibody (goat anti-rabbit IgG (H + L), HRP 1:10,000, Jackson) was added, and the samples were incubated for 1 h, washed with TBST washing solution, reacted with ECL reagent (Millipore), and visualized. The software used for grayscale analysis of the image was Gel Image System ver.4.00 (Tanon, China).

Wilcox test is a non-parametric statistical hypothesis test mainly used for comparison between two groups. Pearson correlation analysis is used to compare the association between two parameters. Univariate and multivariate Cox regression analysis were used to study the relationship between gene expression and the overall survival (OS) rate of patients. The expression level of CGN was divided into a high-risk group and a low-risk group according to the median, and the OS of the patients was analyzed by the Kaplan-Meier method. p < 0.05 was considered statistically significant, and all statistical analyses were performed using the R software (version 4.1.0).

The TCGA-KIRC, GSE46699, GSE36895, and GSE16449 datasets were used to construct a gene co-expression network, and 12, 11, 10, and 11 modules, respectively, were identified in these datasets (Supplementary Figure S1A–D). The module-trait relationship heat map (Figures 1A–D) shows the relationship between each module and the clinical characteristics. The results showed that the salmon module in TCGA-KIRC (R = 0.85, P = 7e-175), the blue module in GSE46699 (R = 0.86, P = 1e-39), the turquoise module in GSE36895 (R = 0.96, P = 2e-28), and the turquoise module in GSE16449 (R = 0.9, P = 3e-26) had the strongest association with normal tissues.

Considering |log₂FC| ≥1.0 and FDR <0.05 as the threshold values, in the TCGA-KIRC dataset, there were 1902 down-regulated genes and 5467 up-regulated genes; in the GSE46699 dataset, there were 522 down-regulated genes and 484 up-regulated genes; in the GSE36895 dataset, there were 868 down-regulated genes and 728 up-regulated genes; in the GSE16449 dataset, there were 1777 down-regulated genes and 1646 up-regulated genes (Figures 2A–D). Subsequently, the intersection of genes was taken from the modules with the highest association obtained from the WGCNA and DEGs analysis. As shown in Figures 3A–C, there were 124 overlapping genes between the differential genes and the co-expression modules.

To determine the biological functions of the intersecting genes, we performed a GO functional enrichment analysis (Supplementary Figure S2) and a pathway analysis (Supplementary Datasheet S1) on the intersecting genes. In the univariate Cox survival analysis, four overlapping genes were observed to be significantly correlated with the OS (FDR <0.05) (Supplementary Table S1). Multivariate Cox survival analysis also confirmed that the expression of the four genes CGN, FECH, UCHL1, and WT1 were independently correlated with the OS of ccRCC patients.

In the KM survival analysis, gene expression was divided into two groups, high and low, according to the median value. The results showed that low expression of CGN had the strongest correlation with a poor prognosis and low OS in ccRCC patients. (Supplementary Table S2 and Figure 4). The integrated datasets of TCGA, GSE3, GSE22541 and GSE29609 also showed a significant correlation between CGN and prognosis (Supplementary Figure S3). Therefore, the CGN gene was selected for further analysis and verification.

We analyzed the relationship between CGN expression and age, gender, tumor grade, tumor stage, T stage, N stage, and M stage. The results showed that CGN expression decreased significantly with age, advanced tumor stage, high grade, and advanced T, N, and M stages (Figure 5).

The CIBERSORT algorithm was used to analyze the proportion of tumor-infiltrating immune subgroups, and 22 immune cell maps of the ccRCC samples and the association matrix of immune cells were constructed (Supplementary Figure S4). To confirm the association between the expression of CGN and immune cells, the Pearson correlation analysis (Figures 6A–G) and the Wilcoxon test (Figure 6H) were performed and the results corresponding to the intersection of these two analyses were used to identify the immune cell types that are most closely related to CGN. (Figure 6I). The results showed that six types of TICs are related to the expression of CGN. Among them, activated NK cells, resting dendritic cells, resting monocytes, and resting mast cells were positively correlated with CGN expression, while CD4⁺ memory activated T cells, regulatory T cells, and M0 macrophages were negatively correlated with CGN expression. These results further support the idea that changes in the expression levels of CGN can affect the immune activity of the tumor microenvironment (TME).

The results showed that the group with a high CGN expression had lower TMB (p < 0.001; Supplementary Figure S5A). Furthermore, a negative association between CGN expression and TMB (R = −0.21, p < 0.001) was also identified (Supplementary Figure S5B). The KM survival curve results for patients with low TMB showed significant clinical benefits and significantly prolonged rates of survival (p < 0.001; Supplementary Figure S5C). After CGN and TMB were grouped together, the ccRCC patients with high TMB and low CGN had the worst prognosis (p < 0.001; Figure 7A). In addition, we divided CGN expression into high and low groups according to the median value, and determined the relationship between these two groups and PD-L1 expression. Figures 7B–D shows that the expression of PD-L1 in patients with high CGN values is significantly increased (p < 0.05).

As Supplementary Figures S6A,B showed, CGN expression was positively association with BAP1 and SETD2 (p < 0.05). In addition, CGN expression was negatively correlated with BIRC5 and CXCR4 (Supplementary Figures S6C,D).

According to the HPA database, the CGN protein level in the tumor tissues was significantly lower than that in the normal tissues (Figures 8A,B). To better study the expression levels of CGN in the normal tissues and tumor tissues of ccRCC patients, 59 samples (Supplementary Tables S3, S4) of each of these tissue types were collected. When compared with the normal tissues, the expression levels of CGN in the tumor tissues were significantly reduced (p < 0.001; Figure 8C). WB results also showed that protein level of CGN in the tumor tissues was significantly lower than that in the normal tissues (p < 0.001; Figure 8D).