Functional Assessment of Four Novel Immune-Related Biomarkers in the Pathogenesis of Clear Cell Renal Cell Carcinoma

Background Clear cell renal cell carcinoma (ccRCC) is the most common subtype of renal cell carcinoma whose pathogenesis is not well understood. We aimed at identifying novel immune-related biomarkers that could be valuable in the diagnosis and prognosis of ccRCC. Methods The Robust Rank Aggregation (RRA) method was used to integrate differently expressed genes (DEGs) of 7 Gene Expression Omnibus (GEO) datasets and obtain robust DEGs. Weighted gene co-expression network analyses (WGCNA) were performed to identify hub genes associated with clinical traits in The Cancer Genome Atlas (TCGA) database. Comprehensive bioinformatic analyses were used to explore the role of hub genes in ccRCC. Results Four hub genes IFI16, LMNB1, RHBDF2 and TACC3 were screened by the RRA method and WGCNA. These genes were found to be up-regulated in ccRCC, an upregulation that could be due to their associations with late TNM stages and tumor grades. The Receiver Operating Characteristic (ROC) curve and Kaplan-Meier survival analysis showed that the four hub genes had great diagnostic and prognostic values for ccRCC, while Gene Set Enrichment Analysis (GSEA) showed that they were involved in immune signaling pathways. They were also found to be closely associated with multiple tumor-infiltrating lymphocytes and critical immune checkpoint expressions. The results of Quantitative Real-time PCR (qRT-PCR) and immunohistochemical staining (IHC) analysis were consistent with bioinformatics analysis results. Conclusion The four hub genes were shown to have great diagnostic and prognostic values and played key roles in the tumor microenvironment of ccRCC.


INTRODUCTION
Renal cell carcinoma (RCC) is the third most common malignancy of the urinary system, accounting for approximately 3% of all malignancies (Siegel et al., 2020). According to the global cancer statistics (2018), there were about 403,262 (2.2%) new cases of RCC with a mortality rate of 1.8% (Bray et al., 2018). Clear cell renal cell carcinoma (ccRCC), accounting for over 80% of all RCC cases, is the most common pathological subtype of RCC (Hsieh et al., 2017). An estimated 30% of all ccRCC cases are diagnosed in the metastatic stages with recrudescence occurring in 20 to 30% of patients who have undergone partial or radical nephrectomy (Cairns, 2010;Wolff et al., 2016;Hsieh et al., 2017). Targeted therapies such as sunitinib (Motzer et al., 2006), sorafenib (Hutson et al., 2010) and axitinib (Hutson et al., 2013) are important as first-line ccRCC medications. These therapies, coupled with immunotherapy have a positive prognostic outcome in ccRCC patients (Rini et al., 2019). Despite the advances in the therapeutic management of ccRCC, the recovery rate of these patients is still low (Vera-Badillo et al., 2015;Tsimafeyeu et al., 2017). Therefore, identify more novel immune-related biomarkers that could be vital in the diagnosis, treatment and prognosis of ccRCC is urgent.
The development of high-throughput technologies and bioinformatic advancements have led to the identification of novel ccRCC biomarkers (Mitchell et al., 2018;Linehan and Ricketts, 2019). However, research using small sample sizes and different sequencing platforms has resulted in great variabilities and poor statistical inferences among studies. The Robust Rank Aggregation (RRA) method can integrate differentially expressed gene (DEG) lists of different datasets, thereby overcoming the challenges posed by small sample sizes (Kolde et al., 2012). In addition, the use of DEGs between different samples, while ignoring the internal relationship between genes, can help overcome this problem. Genes with extremely similar expression patterns in different samples can be identified by the weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008). This analysis screens out biomarkers based on internal relationships among genes and correlates gene sets with their clinical traits (Langfelder and Horvath, 2008).
In this study, we aimed at identifying novel immunerelated biomarkers that were significantly associated with the progression and prognosis of ccRCC. Moreover, we investigated the potential molecular mechanisms of these biomarkers, as well as their associations with the tumor microenvironment (TME). Finally, we performed Quantitative Real-Time PCR (qRT-PCR) analysis to detect the expression of novel immunerelated biomarkers in clinical ccRCC samples.

Data Source and Preprocessing
The workflow of this study is shown in Figure 1. The matrix files of eligible microarray datasets were obtained from the Gene Expression Omnibus (GEO) database 1 . Datasets that had human renal tissue samples and contained at least 10 tumorand non-tumor renal control tissue samples were included in this study. Cell line and xenograft tissues were excluded from this study. Accordingly, 7 GEO datasets were included for DEG analysis, including GSE71963, GSE66270, GSE53757, GSE40435, GSE36895, GSE17895, and GSE16449. The probes were matched to the gene symbols using the annotation files of the respective platforms. Normalized RNA-sequencing data as FPKM (Fragments Per Kilobase Million) and the corresponding pathological information of ccRCC samples were downloaded from The Cancer Genome Atlas (TCGA) database 2 . In total, 517 TCGA-ccRCC and 89 CPTAC-ccRCC samples that had complete clinic-pathological data with follow up time were included into this study.

Identification of Robust DEGs
DEGs between the adjacent normal tissue and ccRCC samples were identified by the 'limma' package (version 3.44.3 3 ), while DEG integration of the 7 microarray datasets to obtain robust DEGs was achieved using the RRA method. In the RRA analysis, | log2-fold change (FC)| and false discovery rate (FDR) < 0.05 were set as the cutoff points for robust DEGs based on the 'RobustRankAggreg' package (version 1.1 4 ). Robust DEGs were further validated between paired ccRCC and adjacent samples in the TCGA-ccRCC database using the 'edgeR' package (version 3.30.0 5 ). Statistical significance was set at | log2FC| > 1 and FDR < 0.05.

Function Enrichment Analysis of Robust DEGs
Gene Ontology (GO) enrichment, including molecular functions (MF), cellular components (CC) and biological processes (BP), and the Kyoto Encyclopedia of Genes and Genome (KEGG) pathway analyses were conducted with the 'clusterprofiler' package (version 3.16.1 6 ). GO terms or KEGG pathways with FDR < 0.05 were visualized by the "GOplot" package.

WGCNA and Identification of the Key Module
The expression data of robust DEGs were retrieved from the TCGA database and used in the WGCNA analysis. The WGCNA method was important in constructing co-expression networks and identifying clinical traits related to DEGs. Pearson's correlations between all pair-wise genes were used to generate the adjacency matrix, while the soft threshold power of β = 5 was used to achieve scale-free topology of the adjacency matrix. The adjacency matrix was then transformed into a topological overlap matrix (TOM). This transformation was done based on the TOM-based dissimilarity measure with a minimum module size of 30 and cut height of 0.25. Robust DEGs with similar expression patterns were classified into the same gene module by average linkage hierarchical clustering. The correlation between module eigengenes (MEs) and clinical traits was calculated to identify

Comprehensive Bioinformatic Analyses of Hub Genes
The TIMER website 7 was used to validate the differences in hub gene expression between pan-cancer and adjacent normal tissues. To identify the diagnostic value of these hub genes in TCGA ccRCC, Receiver Operating Characteristic (ROC) curve analysis was performed and the area under the curve (AUC) was calculated using the "pROC" package (version 1.10.0 8 ).
The 'ggstatsplot' package (version 0.6.1 9 ) was used to assess hub gene expression between different T stages, AJCC stages and tumor grades. Normal tissue protein level and ccRCC data were obtained from the CPTAC database 10 and used to identify the protein level of the hub genes. Moreover, ROC curve analysis was performed to assess the diagnostic value of these proteins.
Exploration of enriched KEGG pathways of the hub genes was achieved using the Gene Set Enrichment Analysis (GSEA) 4.0.1 software. Based on each hub gene's median expression, the 517 ccRCC samples were divided into high-and low-expression groups. "c2.cp.kegg.v7.1.symbols.gmt" as the reference gene set was acquired from the Molecular Signatures Database V7.1 (MSigDB). Statistical significance was set at FDR < 0.05 and | Normalized Enrichment Score (NES)| > 1.

Tumor Immunity Analysis of Hub Genes
The estimate, stromal and immune scores of each TCGA ccRCC sample were downloaded from the ESTIMATE bioinformatics website 11 and used to determine the association between hub genes and tumor purity and the association between the infiltration level of immune cells and the level of stromal cells in ccRCC tissues. The ESTIMATE algorithm is based on a single sample GSEA to evaluate tumor purity.
The TIMER website 7 was used to explore the relationships between hub gene expression and abundance of tumorinfiltrating lymphocytes (TILs) such as CD8 + T cells, CD4 + T cells, B cells, dendritic cells, macrophages and neutrophils (Li et al., 2017). Estimation of 22 TIL compositions from bulk tissues based on their gene expression profiles was performed using the CIBERSORT method (Chen et al., 2018). The LM22 signature matrix was used to identify the 22 TILs containing seven T cell types, natural killer (NK) cells, myeloid subsets, plasma cells and naive and memory B cells. To further identify the relationships between hub genes and TILs, the CIBERSORT website 12 in combination with the FPKM data of TCGA ccRCC and the LM22 signature matrix was used to estimate the TIL fractions. The sum of the 22 estimated TIL fractions in each sample is equal to 1. Spearman rank correlation analysis was used to evaluate the relationships between TILs. Moreover, the TISIDB website 13 was used to explore the associations between hub genes and critical immune checkpoint inhibitors (ICIs: CTLA4, HAVCR2, LAG3, PDCD1, and TIGIT).

RNA Extraction and qRT-PCR
Total RNA was isolated with TRIzol reagent (TaKaRa Biomedical Technology, Dalian, China) according to the manufacturer's instructions. Complementary DNA was reverse-transcribed using the Prime Script RT reagent Kit (TaKaRa). the qRT-PCR analysis was conducted using TB R Green PCR Master Mix (TaKaRa). The specific primers set for IFI16, LMNB1, RHBDF2, TACC3, CD4, CD8 and GAPDH are presented in Supplementary Table 1. All data analyses were managed using RocheLightCycler480. Data were calculated from three biological and technical replicates then normalized to GAPDH expression levels using the 2 − Ct method.

Patients and Clinical Samples
Primary ccRCC patients who underwent radical surgery without any preoperative chemotherapy or radiotherapy at the First Affiliated Hospital of Guangzhou Medical University between 2016 and 2019 were enrolled in the study. As for Formalin fixed paraffin-embedded ccRCC specimens, a total of 150 patients diagnosed with primary ccRCC who underwent operation at the Department of Urology of the First Affiliated Hospital of Guangzhou Medical University (Guangzhou, China) and Nanfang hospital (Guangzhou, China) between February 2008 to August 2015 were enrolled in our study. The follow up of participants (n = 150) were gotten through phone calls until death or the cut-off date of August 1, 2015. The mean follow-up time was 68 months (from 4.0 to 90.0 months). All the deaths were ascribed to ccRCC. Pathological TNM staging was reassessed in accordance with the American Joint Committee on Cancer (AJCC). Histological and pathological diagnoses of the specimens was assigned basing on the 2007 World Health Organization (WHO) Consensus Classification and Staging System of Renal Tumor and Fuhrman grade by two experienced pathologists. Written informed consent was obtained for all patients before specimen collection, following the ethical protocols of the Ethics Committee of the First Affiliated Hospital of Guangzhou Medical University. All the study protocol was approved by the Ethics Committee of the First Affiliated Hospital of Guangzhou Medical University and Southern Medical University Institutional Board.

Tissue Microarray Construction and Immunohistochemistry (IHC)
Tissue microarray (TMA) was established from 150 formalinfixed paraffin-embedded human ccRCC tissues block according to the standard methods. IFI16, RHBDF2, TACC3 protein expression was confirmed using an immunoperoxidase method. The tissue array was deparaffinized, rehydrated, and inhibited Frontiers in Cell and Developmental Biology | www.frontiersin.org endogenous peroxidase activities by 0.3% H2O2 for 30 min. For antigen retrieval, the slides were boiled in sodium citrate buffer (10 mM, pH 6.0) in a pressure cooker for 7 min.
Afterward, non-specific binding was blocked with 5% normal goat serum, and then incubated with primary antibody (IFI16, 1:200, ab13454, Abcam Inc., Cambridge, MA, United States). (RHBDF2, 1:500, abs10947, Absin, Shanghai, China), (TACC3, 1:100, 14970s, Cell Signaling Technology, Inc., United States). Sequentially tissue array was incubated with polyperoxidaseanti-mouse IgG (Zhongshan Biotech. China). Diaminobenzidine (DAB) was visualized as a chromogen substrate. The nucleus was counterstained with hematoxylin. IFI16, RHBDF2, TACC3 staining in nuclear was reckoned as detectable immunoreactions. To evaluate the consequences of IFI16, RHBDF2, TACC3 staining, the intensity and percentage of cells in cancerous and non-cancerous tissues were scored by two pathologists independently. The intensity of staining was determined in accordance with the following scale: 0 (no staining); 1 (weak staining, light yellow); 2 (moderate staining, yellowish brown) and 3 (strong staining, brown). Based on the percentage of positively stained tumor cells, the score of staining extent was denoted on 4point scale as follows: 0 (less than 5%); 1 (5 to 25%); 2 (25 to 50%); 3 (more than 75%). The final scores were then calculated according to score × proportion of positive tumor cells for IFI16, RHBDF2, TACC3 expression (range from 0 to 9). Tumor tissues with scores of 0-1 was recognized as low expression because approximately 90% of normal kidney tissues expressed a low level of IFI16, RHBDF2, TACC3 with an IHC score of ≤ 1 in our preliminary study. Then we defined 0-1 as low expression and 2-9 as high expression.

Statistical Analysis
The hazard ratio (HR) and 95% confidence interval (CI) were calculated by Cox regression analysis and Kaplan-Meier survival analysis. Spearman correlation analysis was used to evaluate the correlation between two continuous variables. The Kruskal-Wallis or student's t-test was used to compare between groups for continuous variables. When multiple comparisons were performed, p-values were corrected according to the FDR method. The FDR method was used to control for multiple testing that could lead to a false positive. All experiments were repeated thrice and the data were presented as means and standard deviation (SD) in all plots shown in this study unless differently stated in the legend. All statistical analyses were conducted by R software (version 3.6.2 16 ) and all P < 0.05 (2-sided) were considered statistically significant.

Functional Enrichment Analysis of Robust DEGs
Significantly enriched BP of robust DEGs was identified, including monovalent inorganic cation homeostasis, small molecule catabolic process, carboxylic acid biosynthetic process and organic acid biosynthetic process (Figure 2A). Several CC GO terms were detected, including apical part of the cell, apical plasma membrane, basolateral plasma membrane and blood microparticle ( Figure 2B). In GO terms of MF, secondary active transmembrane transporter activity, active transmembrane transporter activity and cofactor binding were significantly enriched terms ( Figure 2C). Based on KEGG pathway analysis, glycolysis/gluconeogenesis, PPAR signaling pathway and collecting duct acid secretion were mostly associated with the robust DEGs ( Figure 2D).

WGCNA and Identification of Key Modules and Hub Genes
To identify key modules significantly related to ccRCC clinical traits, WGCNA was performed on the TCGA-ccRCC dataset incorporating 841 robust DEGs derived from the previous analysis (Figure 3). Clinical information of ccRCC patients such as age, overall survival status (OSS), overall survival time (OST), disease-free status (DFS), disease-free time (DFT), T stage, N stage, M stage, AJCC stage and tumor grade were retrieved from TCGA ( Figure 3A). By setting the cut height at 0.25 and β at 5 (scale-free R2 = 0.85), 841 robust DEGs were divided into six modules (Figures 3B-E). As shown in the heatmap of the module-trait relationship, the brown module was significantly related to clinical traits ( Figure 3F). Gene significance > 0.3 of AJCC stage and tumor grade and MM > 0.7 were selected as cutoff points (Figures 3G,H). We identified 4 hub genes from the brown module: IFI16, LMNB1, RHBDF2 and TACC3.

Comprehensive Bioinformatic Analyses of Hub Genes
As shown in the pan-cancer view (Supplementary Figure 2), the 4 hub genes (IFI16, LMNB1, RHBDF2, and TACC3) were significantly up-regulated in ccRCC samples and other cancer types when compared to adjacent normal tissues (p < 0.001). Furthermore, the ROC curve analysis showed that these hub genes had a high diagnostic value as biomarkers for TCGA ccRCC (IFI16 AUC: 0.921, LMNB1 AUC: 0.87, RHBDF2 AUC: 0.957, TACC3 AUC: 0.896; Supplementary Figure 3A). These genes were significantly differentially expressed in ccRCC samples with different T stages, AJCC stages and tumor grades (p < 0.001). Higher expression levels were an indication of advanced T stages, AJCC stages and tumor grades (Figures 4A-C).
The protein levels of the 4 hub genes were significantly up-regulated in ccRCC samples compared to normal tissues (p < 0.001) (Figure 4D). The ROC curve analysis revealed that these genes had good efficacies in the diagnosis of ccRCC (IFI16 protein AUC: 0.955, LMNB1 protein AUC: 0.959, RHBDF2 protein AUC: 0.94, TACC3 protein AUC: 0.915; Supplementary Figure 3B).
Further investigations of enriched KEGG pathways of IFI16, LMNB1, RHBDF2 and TACC3 in ccRCC showed that highly expressed genes [IFI16 (A), LMNB1 (B), RHBDF2 (C), and TACC3 (D)] were all enriched in T cell receptor signaling pathway, NK cell-mediated cytotoxicity, antigen processing and presentation and primary immunodeficiency. Meanwhile, the NOD-like receptor signaling, cytosolic DNA sensing and Toll-like receptor signaling pathways were enriched in the highexpression groups of IFI16, LMNB1 and RHBDF2. The cell cycle and B cell receptor signaling pathway were enriched in LMNB1 and IFI16 high-expression groups, respectively. The GSEA showed that LMNB1, IFI16, RHBDF2 and TACC3 were closely associated with immune signaling pathways (Figures 5A-D).

Tumor Immunity Analysis of Hub Genes
IFI16, LMNB1, RHBDF2 and TACC3 were positively correlated with the estimate, stromal, and immune scores (Supplementary  Figure 4). These results suggest that the 4 hub genes were negatively correlated with the tumor purity of ccRCC and were up-regulated in the TME. The expression levels of IFI16, LMNB1, RHBDF2 and TACC3 positively correlated with infiltration levels of the six TILs, including CD8 + T cells, CD4 + T cells, B cells, dendritic cells, macrophages and neutrophils (Supplementary Figure 5). This shows that these genes play a key role in the immune infiltration of ccRCC. As shown in Figure 6A, abundant fractions of 22 TILs were different in each ccRCC sample. This explains the tumor heterogeneity among different individuals. In addition, different TIL subpopulation ratios were weakly to moderately correlated ( Figure 6B). The analysis showed that these genes were positively correlated with multiple TILs, especially activated CD4 + memory T cells, CD8 + T cells, regulatory T cells (Treg) and follicular helper T (Tfh) cells, but were negatively correlated with resting mast cells, resting NK cells and activated NK cells (Figure 6C). There was also a positive correlation between the 4 hub genes and the expression levels of TIGIT, HAVCR2, CTLA4, PDCD1 and LAG3 in ccRCC, revealing that these genes might be associated with the immunosuppressive microenvironment (Figure 7).  (Figures 8A-D), which were consistent with the Kaplan-Meier survival analysis results of CPTAC-ccRCC cohort (Supplementary Figure S6). In the univariate Cox proportional hazards regression analysis, 11 variables, included T stage, M stage, N stage, age, tumor grade and AJCC stage, showed statistical significance with hub gene expression. Based on the multivariate Cox proportional hazards regression analysis, the 4 hub genes were regarded as independent prognostic factors for ccRCC ( Table 2).

Validation of the Expression of Hub Genes in Clinical ccRCC Specimens
To detect the expression of the 4 hub genes (IFI16, LMNB1, RHBDF2, and TACC3) in ccRCC, we performed the qRT-PCR analysis in clinical specimens. The clinicopathological information of 15 ccRCC patients is shown in Supplementary  Table 2. The mRNA expression of all the hub genes was significantly higher in ccRCC tissues when compared with adjacent normal tissues (Figure 9). This is consistent with the results of our bioinformatics analysis. These findings suggested that the expression of the hub genes may act as a promising biomarker for ccRCC. However, we don't found closely correlation between the 4 signatures and the immune cells biomarkers (CD4 and CD8) in ccRCC sample (Supplementary Figure 7). To further confirm the correlation of hub genes expression with tumor progression in ccRCC. We performed immunohistochemical staining (IHC) in human tissue samples to detected the protein expression of IFI16, RHBDF2, TACC3. Patients' characteristics were retrospectively collected from the review of medical records and detailed in Supplementary  Table 3. Results demonstrated that IFI16 was increased significantly in most of the paired ccRCC tissues compared with adjacent normal kidney tissues ( Figure 10A and Table 3, p < 0.0001). However, RHBDF2 was decreased in the paired ccRCC tissues than adjacent normal kidney tissues (Supplementary Figure 8 and Supplementary Table 4, p < 0.0001) and TACC3 was negative in most tissues (Supplementary Figure 8 and Supplementary Table 5). Subsequently, we analyzed the clinical correlation of IFI16 and clinicopathological characteristics. As shown in Table 3, upregulation of IFI16 was significantly associated with pathology grade (p < 0.0001). However, no significant correlation was found between IFI16 protein expression with other clinical features. Kaplan-Meier analysis for 150 patients with follow-up data suggested that patients with higher levels of IFI16 presented significantly lower overall survival rates than those with low levels of IFI16 expression ( Figure 10B, Log rank, p = 0.046). Furthermore, the multivariate analysis of the Cox regression FIGURE 9 | The expression of these Hub genes in human ccRCC specimens and adjacent normal tissues (ANT). (A-D) q-RTPCR analysis of IFI16 (A), TACC3 (B), LMNB1 (C), and RHBDF2 (D) in paired ccRCC tissues (n = 15). GAPDH was used as a loading control. *p < 0.05, **p < 0.01, ***p < 0.001. model, IFI16 expression (p = 0.0226, HR = 5.474) was confirmed to be independent prognosis factors for ccRCC patients (Table 4). To further investigate the expression profile of IFI16 in human ccRCC, we detected IFI16 expression in four human ccRCC cell lines (i.e., SW839, OSRC-2, SM12-PN6, Caki-1, and 786-O) and found that the mRNA levels of IFI16 was relatively higher in most of the ccRCC cells ( Figure 10C). These results suggested that IFI16 may be a probable independent predictor in patients with ccRCC.

DISCUSSION
Clear cell renal cell carcinoma is a complex and highly heterogeneous disease whose pathogenesis remains unclear (Hsieh et al., 2017). Therefore, understanding the potential molecular mechanisms of ccRCC is crucial for better diagnosis, treatment and prognostic predictions (Linehan and Ricketts, 2019). Although previous studies have used high-throughput technologies and advanced bioinformatics to find novel biomarkers and therapeutic targets for ccRCC, there are inconsistencies among the DEGs analyzed in different studies (Barbieri et al., 2017).
In this study, 7 GEO datasets were integrated using RRA to minimize inconsistencies and identify robust DEGs. Enrichment analyses to explore potential biological functions of robust DEGs in ccRCC were also performed, after which hub genes associated with the pathogenesis of ccRCC were screened. Four hub genes were finally screened from the two phenotypes. Comprehensive bioinformatic analyses of the four hub genes showed that they were closely associated with the pathogenesis of ccRCC. Their protein levels were significantly up-regulated in ccRCC samples when compared to normal tissues. This high expression coincided with a poor OS, hence, these genes could be regarded as independent ccRCC prognostic factors. The GSEA analysis showed that these genes were closely related to immune signaling pathways. To the best of our knowledge, this is the first study combining RRA with WGCNA to explore hub genes involved in ccRCC pathogenesis.
Robust DEGs, such as AQP9 (Xu et al., 2019) and SULT (Li et al., 2019), are biomarkers of ccRCC and play a key role in its pathogenesis. Based on GO enrichment results of robust DEGs, studies have proved that small molecule catabolic processes and carboxylic acid biosynthetic processes are significantly associated with the initiation and progression of cancer (Sciacovelli and Frezza, 2016;DiNatale et al., 2020). Selvakumar et al. (2014) reported that active transmembrane transporter activity plays a key role in ccRCC, which was an enriched GO term of robust DEGs. Enrichment of robust DEGs in some KEGG pathways such as glycolysis/gluconeogenesis (Massari et al., 2015;Ciccarese et al., 2016) and PPAR signaling (Chang and Lai, 2019) also prove their relevance in ccRCC pathogenesis. Documented evidence suggests that ccRCC cells use the glycolytic pathway for energy production in the presence of oxygen, a phenomenon known as the Warburg Effect (Massari et al., 2015;Ciccarese et al., 2016). Dysregulated PPAR signaling pathway in pancancer results in dysregulated metabolism and is associated with immunosuppression (Chang and Lai, 2019). We noted multiple molecular mechanisms of robust DEGs that were closely associated with ccRCC pathogenesis.
Four hub genes (IFI16, LMNB1, RHBDF2, and TACC3) were identified in this study. Interferon-inducible 16 (IFI16), an innate immune sensor for DNA in cells, can recruit STING after DNA stimulation (Unterholzner et al., 2010). This interferon can also activate the STING signaling pathway that plays a key role in the immune escape, thereby promoting tumor progression (Lemos et al., 2016;He et al., 2017). Cai et al. reported that IFI16 promotes the progression of cervical cancer by up-regulating PD-L1 in TME through the STING-TBK1-NF-kB pathway (Cai et al., 2020). Lamin B1 (LMNB1) is associated with various cellular physiological activities, including nuclear autophagy, DNA replication and transcription, nuclear migration, DNA repair pathways, etc. (Barascu et al., 2012;Yang et al., 2019). Previous studies showed that LMNB1 was overexpressed in pancreatic cancer, liver cancer and prostate cancer (Butin-Israeli et al., 2012;Irianto et al., 2016) and its overexpression was associated with poor clinical outcomes in the cervical (Yang et al., 2019) and colon (Izdebska et al., 2018) cancers. These findings were consistent with our results. Rhomboid 5 homolog 2 (RHBDF2) induces gastric cancer cell invasiveness by regulating Transforming Growth Factor Beta 1 (TGFB1) signaling (Ishimoto et al., 2017), a finding that corresponds with our results. Mutations in RHBDF2 accelerate tumorigenesis by activating epidermal growth factor receptor (EGFR) signaling (Hosur et al., 2014) and are associated with tylosis esophageal cancer (Blaydon et al., 2012). Transforming acidic coiled-coil protein 3 (TACC3) is overexpressed in RCC cells and can promote proliferation, invasion and migration of RCC cells (Guo and Liu, 2018). Overexpression of TACC3 is associated with poor prognosis in the breast (Song et al., 2018), prostate (Qie et al., 2020) and colorectal (Du et al., 2016) cancers. Commonly occurring gene fusions such as FGFR3-TACC3 fusions are potent oncogenes that rely on mitochondrial respiration (Frattini et al., 2018). This finding is similar to the finding in our study. However, there is no documented evidence on the roles of the four hub genes (IFI16, LMNB1, RHBDF2, and TACC3) in ccRCC pathogenesis. Furthermore, we performed the qRT-PCR analysis in clinical samples and found that the mRNA expression of the four hub genes was significantly higher in ccRCC tissues when compared with adjacent normal tissues. This study highlights the roles of the four hub genes in ccRCC; however, more in vivo and in vitro experiments are needed to authenticate these findings.
Clear cell renal cell carcinoma is a highly immune-infiltrated tumor and its pathogenesis was closely associated with TME (Galon and Bruni, 2019). Immunotherapy plus targeted therapy are the new conventional approaches for systemic treatment of metastatic ccRCC (Rini et al., 2019). Our study shows that the four hub genes are involved in multiple immune-related signaling pathways and they positively correlate with Estimate scores. This shows that they were overexpressed in the TME. We hypothesized that the effect of these genes on tumor prognosis is related to tumor immunity. Based on this aspect, we found that these genes were positively correlated with multiple TILs, especially Treg cells, activated CD4 + memory T cells and CD8 + T cells. Multiple pieces of research documented that the high abundance of CD4 + T cells and CD8 + T cells in ccRCC were related to its pathogenesis and poor prognosis (Nakano et al., 2001;Remark et al., 2013). Increased Treg cell numbers can suppress anti-tumor immune responses and are correlated with poor ccRCC prognosis (Liotta et al., 2011;Kang et al., 2013;Polimeno et al., 2013). There was a negative association between resting NK and activated NK cells and the four hub genes. Low NK-cell densities were related to a worse prognosis in ccRCC (Remark et al., 2013). These results reveal that TACC3, RHBDF2, LMNB1 and IFI16 may promote tumor progression by regulating TILs in ccRCC. However, exploring the molecular mechanisms by which these genes regulate TILs in ccRCC will be significant.
Among the ICIs, cytotoxic T-lymphocyte-associated antigen 4 (CTLA4), the godfather of checkpoint inhibitors, can dampen early activation and differentiation of T cells and actively send inhibitory signals to T cells (Pardoll, 2012;Buchbinder and Desai, 2016). Programmed cell death protein 1 (PD1; also known as PDCD1) is highly expressed on Treg cells of various cancers, thereby suppressing T cell effector functions (Pardoll, 2012;Hellmann et al., 2016). Lymphocyte activation gene 3 (LAG3) and T cell membrane protein 3 (TIM3; also known as HAVCR2) have been linked to the inhibition of lymphocyte activity (Hellmann et al., 2016;Du et al., 2017). Blocking these receptors has been shown to strengthen anti-tumor immunity in tumor animal models (Pardoll, 2012). TIGIT, T cell immunoglobulin and ITIM domain, is an inhibitory immunoreceptor and an interesting cancer immunotherapy target (Manieri et al., 2017). We report that IFI16, LMNB1, RHBDF2 and TACC3 expression was positively associated with ICIs. A high expression of these genes coincided with poor ccRCC prognosis, as they likely promote ccRCC progression by tumor immune escape.
Our study provides new insights into immunotherapy and TME in ccRCC. However, there were some limitations in this study. First, retrospective study designs induce heterogeneity in results, thus, more in vivo and in vitro experiments should be performed to validate our findings. Second, the biological mechanisms of TACC3, LMNB1, RHBDF2, and IFI16 identified in this study warrant further investigation.
In conclusion, there are robust DEGs and several gene modules that are associated with the clinical pathogenesis of ccRCC. Four hub genes (IFI16, LMNB1, RHBDF2 and TACC3) were up-regulated in ccRCC tissues and correlated with ccRCC progression. These genes were associated with poor prognosis of ccRCC and may play key roles in the TME of ccRCC by regulating TILs or ICIs. However, in vivo and in vitro experiments are needed to validate the contribution of these genes in the pathogenesis of ccRCC. Additionally, more studies should be conducted to understand the molecular mechanisms of IFI16, LMNB1, RHBDF2 and TACC3 in the pathogenesis of ccRCC.

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 author/s.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by First Affiliated Hospital of Guangzhou Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
DL, XW, MW, WC, and SY performed data analyses. DL performed the experiments. MW, WC, and SY assisted in data collection. XW and DL wrote the manuscript. Funding was obtained by GZ, YL, and DL. GZ, YL, and DG supervised the study and revised the manuscript. All authors participated in preparing the manuscript and approved the final version of the manuscript.