CBX Family Members in Two Major Subtypes of Renal Cell Carcinoma: A Comparative Bioinformatic Analysis

The biological function and clinical values of Chromobox (CBX) family proteins in renal cell carcinoma (RCC) are still poorly investigated. This study aimed to compare the expression profiles and clinical relevance of CBXs between the two most frequent subtypes of RCC, clear cell renal cell carcinomas (ccRCC) and papillary renal cell carcinomas (pRCC), and to investigate whether CBXs would play a more or less similar role in the pathogenesis and progression of these RCC subtypes. Considering these two RCC populations in the TCGA database, we built a bioinformatics framework by integrating a computational pipeline with several online tools. CBXs showed a similar trend in ccRCC and pRCC tissues but with some features specific for each subtype. Specifically, the relative expressions of CBX3 and CBX2 were, respectively, the highest and lowest among all CBXs in both RCC subtypes. These data also found confirmation in cellular validation. Except for CBX4 and CBX8, all others were deregulated in the ccRCC subtype. CBX1, CBX6, and CBX7 were also significantly associated with the tumor stage. Further, low expression levels of CBX1, CBX5, CBX6, CBX7, and high expression of CBX8 were associated with poor prognosis. Otherwise, in the pRCC subtype, CBX2, CBX3, CBX7, and CBX8 were deregulated, and CBX2, CBX6, and CBX7 were associated with the tumor stage. In addition, in pRCC patients, low expression levels of CBX2, CBX4, and CBX7 were associated with an unfavorable prognosis. Similarly, CBX3, CBX6, and CBX7 presented the highest alteration rate in both subtypes and were found to be functionally related to histone binding, nuclear chromosomes, and heterochromatin. Furthermore, CBX gene expression levels correlated with immune cell infiltration, suggesting that CBXs might reflect the immune status of RCC subtypes. Our results highlight similarities and differences of CBXs within the two major RCC subtypes, providing new insights for future eligible biomarkers or possible molecular therapeutic targets for these diseases.


Introduction
Renal cell carcinoma (RCC) is an insidious malignancy. It is currently the tenth most diagnosed tumor in women and the sixth in men, making up about 4% of all malignant tumors [1]. With a 76% 5-year relative survival rate, RCC is the deadliest urological cancer, although its survival rate has markedly improved since a low of 46.8% in 1977. The survival rate strictly depends on the stage at diagnosis, with a 5-year relative survival of 93% for stage I, localized disease, 72.5% for stage II/III regional disease, and only 12% for stage IV metastatic disease [2]. Although the diagnosis and management of RCC have improved remarkably [3], its incidence is expected to increase worldwide [4]. Therefore, the disease still warrants further investigation for better prevention and management. RCC is a heterogeneous disease comprising a collection of different subtypes, characterized by distinct genetic abnormalities and molecular signatures reflecting the differences in the cell type, biology, and underlying molecular mechanisms [5]. The most frequent renal tumor subtypes are clear cell renal cell carcinomas (ccRCC) and papillary renal cell carcinomas (pRCC), accounting for 70-80 % and 10-15 % of all the RCCs [2], respectively. Although Libraries were prepared by using the Illumina mRNA TruSeq kit following the manufacturer's directions and sequenced on the Illumina HiSeq 2000. mRNA sequence reads were aligned to the reference genome, and gene expression was quantified based on the gene models defined in the TCGA Gene Annotation File (GAF).

Power Analysis
A posteriori power analysis was performed to find out what power would be obtained for the specified effect size and sample size. Considering two groups (normal and tumor), we performed a two-tailed Wilcoxon-Mann-Whitney test with an effect size d = 0.5 and an alpha level of 0.05, obtaining a power of 0.97 for KIRC and 0.74 for KIRP. By considering more than two groups (4 tumor stages), we performed a one-way ANOVA with an effect size d = 0.25, an alpha level of 0.05, and a total sample size of 200, obtaining a power of 0.84.

GEPIA2
Gene Expression Profiling Interactive Analysis (http://gepia2.cancer-pku.cn/ (accessed on 20 December 2021)) is a user-friendly and interactive web resource for analyzing the RNA sequencing expression data of 9736 tumors and 8587 normal samples from the TCGA and the GTEx projects, using a standard processing pipeline [23]. We used GEPIA2 to analyze the tissue-wise expression of eight CBX family members in ccRCC and pRCC between tumor/normal tissues. The TCGA dataset was available within GEPIA2 (Kidney renal clear cell carcinoma, KIRC: 523 and Kidney renal papillary cell carcinoma, KIRP: 286). We selected only TCGA normal for differential analysis and plotting. We also profiled the expression of CBX genes based on the patient pathological stage. Statistical significance for differential gene analysis was estimated by one-way ANOVA and setting p value < 0.05. We also performed a multiple-gene comparison and obtained an expression matrix plot of CBX genes in two RCC subtypes.
We also used GEPIA2 to perform survival analysis by Mantel-Cox test to compare the survival contribution, in terms of overall survival (OS) and disease-free survival (DFS), of CBX genes in RCC subtypes. The median expression threshold was considered for splitting the high-expression and low-expression cohorts and a p value cut-off of 0.05 to determine significance. For CBXs, which gave significant results, we also performed survival analysis based on the expression status of CBX genes and plotted a Kaplan-Meier curve.

TIMER2.0
Correlations between tumor-infiltrating immune cells (TIICs) and CBXs were assessed by the Tumor Immune Estimation Resource platform (TIMER2.0, http://timer.cistrome.org/ (accessed on 20 January 2022)), which provides dedicated functions for analysis and visualization of TIICs across multiple cancer types associated with their clinical impact [24,25]. We selected CBXs as genes of interest via the "Gene module" to generate a functional heatmap that reported the correlation between CBXs expression values, from KIRC and KIRP TCGA datasets, with the immune infiltration cells. TIICs included B cells, CD4 + T cells, CD8 + T cells, dendritic cells, macrophages, and neutrophils. "Purity Adjustment" based on the partial Spearman's correlation was selected to perform this association analysis.

cBioPortal
The cBioPortal for Cancer Genomics Portal (https://www.cbioportal.org/ (accessed on 21 January 2022)) is a free tool for interactive exploration of multidimensional cancer genomics datasets. We used it to analyze the frequency of CBX gene alteration and coexpression of eight CBX family members. KIRC-TCGA PanCancer (in total 512 samples) and KIRP-TCGA PanCancer (in total 283 samples) data were selected for analysis. A z-score with the range ±2 between tumor and normal conditions was chosen for mRNA expression values. For the correlation analysis of all possible combinations for the CBXs family, Spearman's correlation coefficient >0.3 and the p value < 0.05 were considered for significance.

Protein-Protein Interaction (PPI) Network and Functional Enrichment Analysis
The protein-protein interaction network was generated with the GeneMania Cytoscape plugin [26]. As CBX3/6/7 were found to be the genes with the highest alteration rates in both RCC subtypes, we decided to focus network analysis only on these genes. Physical and co-expression interactions were selected, and the enrichment analysis (FDR ≤ 0.05) based on Gene Ontology (GO) Biological Process was performed. Pathway analysis was performed by using Reactome (REAC) database implemented in g:Profiler [27]. The threshold of 0.05 was considered for the statistical significance after the Bonferroni correction.

Cell Lines Validation
To validate our results in cell lines, we downloaded expression data from the Cancer Cell Line Encyclopedia (CCLE) (https://depmap.org/portal/ (accessed on 22 August 2022)), an extensive collection of expression and genetic data for human cancer cell models [28]. The portal contains expression data from non-cancerous cell lines, clear cell renal carcinoma, and renal cell carcinoma. Each cell line is identified inside this data base with DepMap ID (for details see Table S1). Expression data are reported as mean ± standard deviation.

CBX Family mRNA Transcriptional Levels and Clinicopathological Parameters
GEPIA2 was used to analyze the transcriptional levels of eight CBX genes in ccRCC and pRCC tissues compared with the normal kidney tissues. Clinical data of the selected populations, based on the TCGA database (KIRC and KIRP TCGA datasets), are summarized in Table 1. In ccRCC tissues compared with normal samples, the CBX1/2/5/6/7 were significantly downregulated in expression level, the mRNA expression of CBX3 was upregulated considerably, and no significant expression difference was found for CBX4/8 ( Figure 1). In pRCC tissues, the CBX2/7 were significantly downregulated in expression level, the mRNA expression of CBX3/8 was significantly upregulated, and no significant expression difference was found for other CBX members ( Figure 1). We also evaluated the relative expression levels of CBXs in comparing the two pathological conditions. We found that the CBXs relative expression followed a similar trend for both tumor subtypes: CBX3 was the highest among all CBXs, while CBX2 was the lowest (Figure 2). Additionally, the association of the expression of CBX genes with the tumor stage was analyzed. Specifically, the median expression of CBX2 showed an evident increase with more advanced tumor stages in pRCC, while the median expression of CBX7 showed decreasing trends in both RCC subtypes ( Figure 3). Although without a net trend, the mRNA expression distribution of CBX1 in clear cell subtypes and CBX6 in both subtypes were significantly different among tumor stages ( Figure 3).

Prognostic Value of CBXs Family mRNA Expression Levels in RCC Patients
To screen for the prognostic impact of CBXs between RCC subtypes, the "Survival Map" module of GEPIA was used. In this way, we compared the survival contribution of CBX family members in ccRCC and pRCC in terms of OS and DFS. As shown in Figure 4, we denoted higher risk associated with CBX8 and lower risk with CBX1/5/6/7 in ccRCC. Similarly, CBX7 was associated with lower risk also in pRCC, while CBX2 and CBX4 were associated with a higher risk in the papillary subtype. Specifically, the increased CBX8 and the decreased CBX1/5/6/7 expression values were significantly associated with shorter OS. Moreover, the low expression of CBX1 and CBX7 also had significantly unfavorable effects on DFS. Otherwise, in the pRCC group, the increased CBX4 and the decreased CBX7 expression levels were significantly associated with shorter OS. Furthermore, the high expression of CBX2 and CBX4 had significantly unfavorable effects on DFS.

Genetic Alteration and Co-Expression of CBX Family in RCC
Using cBioPortal, we analyzed genetic alterations of the eight CBX proteins and their correlations with each other in ccRCC and pRCC populations. Results showed that the CBXs were altered in 76% of patients affected by ccRCC and in 99% of patients affected by pRCC ( Figure 5A,B, respectively). Among eight CBXs, CBX3/6/7 were most altered in both RCC subtypes. Moreover, in both subtypes, the primary alteration type for CBX3 regarded high mRNA expression, while for CBX6/7 low mRNA expression. Diagnostics 2022, 12, x FOR PEER REVIEW 6 of 21        Using cBioPortal, we analyzed genetic alterations of the eight CBX proteins and their correlations with each other in ccRCC and pRCC populations. Results showed that the CBXs were altered in 76% of patients affected by ccRCC and in 99% of patients affected by pRCC ( Figure 5A,B, respectively). Among eight CBXs, CBX3/6/7 were most altered in both RCC subtypes. Moreover, in both subtypes, the primary alteration type for CBX3 regarded high mRNA expression, while for CBX6/7 low mRNA expression.

Immune Cell Infiltration of CBXs in RCC Patients
To explore whether the expression of CBX genes was associated with the immune infiltration of RCC, we used the TIMER2.0 database ( Table 2

Immune Cell Infiltration of CBXs in RCC Patients
To explore whether the expression of CBX genes was associated with the immune infiltration of RCC, we used the TIMER2.0 database ( Table 2)   Legend: The color pink indicates a significant positive correlation (p < 0.05, ρ > 0), the color blue indicates a significant negative correlation (p < 0.05, ρ < 0), and the color grey indicates no significant correlation (p > 0.05).

Protein-Protein Interaction (PPI) Network and Functional Enrichment Analysis
We decided to construct a functional protein-association network focusing on a subset of three CBX family members, specifically CBX3/6/7. For both assessed RCC subtypes, these three CBXs presented a similar profile: CBX3/6/7 showed the highest rates of genetic alteration, CBX3 mRNA level was upregulated, CBX6/7 associated with tumor stage, and CBX7 also associated with better prognosis. We first explored the top 20 genes linked to CBX members and then submitted the network for enrichment analysis. As shown in Figure 7, the PPI network was composed of 26 nodes and 138 edges. The top 20 functional partners were as follows: KMT5C, BMI1, RNF2, PHC2, PHC3, RING1, PCGF6, PCGF3, SYAP1, H3C13, PCGF5, LBR, H3-3A, FAF1, GPR173, PHC1, CBX8, H3C1, MACROH2A1, SP100. The GeneMANIA results also revealed that this network was most enriched for biological processes in protein ubiquitination, protein sumoylation, negative regulation of transcription, chromatin and nucleosome assembly, and G0 to G1 transition, while for cellular components, it was significantly enriched in PcG protein complex, chromatin, and nucleosome organization. From the Reactome analysis, we found that six pathways, including HSA-8939243 (RUNX1 interacts with co-factors whose precise effect on RUNX1 targets is not known),

Cell Lines Validation
To validate our results, we used public kidney cell line data in CCLE. We evaluated the expression data from two non-cancerous cell lines, ten clear cell renal carcinoma and twenty-two renal cell carcinoma (Table S1). We observed that the expression trend of the tissues was also confirmed in cell lines (Figure 8). Specifically, we found that CBXs expression followed a similar trend for both tumor cell types (clear cell renal carcinoma and renal cell carcinoma). We also confirmed that CBX3 showed the highest expression values among all CBXs. In contrast, CBX2 showed the lowest ones ( Figure 8).

Cell Lines Validation
To validate our results, we used public kidney cell line data in CCLE. We evaluated the expression data from two non-cancerous cell lines, ten clear cell renal carcinoma and twenty-two renal cell carcinoma (Table S1). We observed that the expression trend of the tissues was also confirmed in cell lines (Figure 8). Specifically, we found that CBXs expression followed a similar trend for both tumor cell types (clear cell renal carcinoma and renal cell carcinoma). We also confirmed that CBX3 showed the highest expression values among all CBXs. In contrast, CBX2 showed the lowest ones ( Figure 8).  Table S1). For each cell type, mean and standard deviation are also reported.

Discussion
The critical role of epigenetic alterations in the development and progression of RCCs has been depicted [22,29,30]. However, epigenetic biomarkers for RCC are still not used in clinical practice, and targeted epigenetic therapies are under investigation [29]. As the critical components of the epigenetic regulatory complexes, CBX family proteins are involved in the development and progression of multiple cancers [18,31,32]. However, the distinct roles of the CBXs family in RCC have yet to be clarified. Recently, two papers have explored the clinical significance and biological function of CBXs in ccRCC. A first paper assessed the molecular mechanism of CBX4 on the development and progression of clear cell subtype in vitro and in vivo models [33], and a second used a bioinformatic approach to investigate the expression and prognostic value of CBXs [34]. To the best of our knowledge, no study has shown whether such similarities and differences exist in RCC subtypes. Based on the rapid development of second-generation sequencing technology and the establishment of a large number of databases, a comprehensive study of the CBXs family in RCC subtypes could help discover new biomarkers with diagnostic and prognostic potential or new therapeutic targets for this deadly disease. Therefore, in this work, we performed an integrated in silico approach using several bioinformatics tools to provide a comprehensive and comparative analysis of CBXs in the two most frequent RCC subtypes.
Overall, the relative expression levels of CBXs in ccRCC and pRCC subtypes showed a similar trend, with CBX3 and CBX2 having, respectively, the highest and lowest among all CBXs. However, the relevance and involvement of CBX members in the clinical evolution of the disease did not always show as much concordance in the two investigated tumor subtypes.
CBX1 has been reported to be upregulated in several cancers, such as breast [35], colorectal [36], and gastric cancer [37]. It has been detected to have higher expression in hepatocellular carcinoma (HCC), associated with worse clinical outcomes, and functions as an oncogene by interacting with the transcription factor HMGA2 to activate the Wnt/β-catenin signaling pathway in hepatocellular carcinoma [38]. Contrary to the protumor effect assessed in HCC, here, we disclosed a potential antitumor effect of CBX1 in ccRCC, consistent with recent results reported by Zhu and colleagues [34]. Its significant deregulation was related to tumor stage and worse prognosis in ccRCC patients. In the papillary subtype, CBX1 did not seem to play a specific role in the disease because its expression profile appeared similar in tumoral and normal tissues and was not associated with progression or the prognosis of patients.
The overexpression of CBX2 has been reported in breast cancer, where it promoted the progression of the disease by the PI3K/AKT signal pathway [39]. Again, CBX2 has been found to regulate proliferation and apoptosis by phosphorylating YAS in hepatocellular cells [40]; it promoted the proliferation, invasion, and migration of gastric cancer by activating the YAP/beta-catenin pathway [41]. Our results highlighted CBX2 as the most down-expressed gene among all CBX family members in both RCC subtypes. Additionally, CBX2 provided prognostic information in the pRCC population, as it was associated with tumor stage and worse clinical outcomes.
Many studies have reported elevated expression of CBX3 in different human cancer tissues and proposed it as a predictor of unfavorable prognosis for malignancies [42][43][44]. High CBX3 expression promoted glioma proliferation by targeting CDKN1A [44]. Moreover, CBX3 binding activity to methylated histone H3K9 was required to promote the proliferation of colorectal cancer [45] and lung cancer [46]. Recently, some researchers found that high expression of CBX3 promoted ovarian cancer cell proliferation and was significantly related to unfavorable progression and chemoresistance, impacting the treatment outcomes of OV patients [47]. Our results showed similar biological behavior of CBX3 in both RCC subtypes. CBX3 emerged as the most overexpressed gene among all CBX family members in both RCC subtypes in ccRCC and pRCC. However, its overexpression did not correlate with cancer stages or patient prognosis for both subtypes.
CBX4 was overexpressed in breast cancer [48] and osteosarcoma [49]. Functionally, it promoted hepatocellular carcinoma via SUMOylating HIF-1α to upregulate vascular endothelial growth factor [50]. Silencing of CBX4 inhibited cell growth and metastasis in lung cancer by regulating the BMI-1 pathway [51]. More recently, the oncogenic activity of CBX4 has been demonstrated in ccRCC by inhibiting tumor suppressor KLF6 via interaction with HDAC1 [33]. CBX4 overexpression promoted tumor growth and metastasis, whereas CBX4 silencing resulted in the opposite phenotypes. In our study, CBX4 transcriptional level did not show significant changes between tumor and normal conditions in RCC patients. However, in the pRCC population, decreased CBX4 correlated with poor prognosis (shorter OS and DFS).
Regarding CBX5, its upregulation has been associated with increased cell proliferation and poor clinical prognosis [52,53]; furthermore, its downregulation was linked to a higher invasive potential of cancer cells in metastatic cells of colon cancer and thyroid carcinoma compared with non-metastatic cells [53]. Therefore, CBX5 deregulation can play dual mechanistic functions for cancer cell proliferation and metastasis suppression, and the underlying cellular mechanisms are not yet comprehensively known. Here, CBX5 mRNA was downregulated in ccRCC and was associated with a shorter survival time. In the papillary subtype, CBX5 transcriptional level was not significantly deregulated in tumor tissues and did not associate with other clinical parameters.
Our findings supported a different potential involvement of CBX6 in ccRCC and pRCC, and a similar tumor suppressive role for CBX7 in both assessed RCC subtypes. However, to date, the roles of these two CBX members appear unclear: depending on the tumor type or cellular context, they can act as oncogenes or tumor suppressors. For instance, significant downregulation of CBX6 was reported in breast cancer, where ectopic overexpression inhibited tumor cancer progression [54]. Otherwise, in liver cancer, high expression of CBX6 was associated with a worse prognosis [55]. In addition, regarding CBX7, experimental evidence provides contradictory results. CBX7 decreased expression has been reported for many carcinomas, such as breast [56], lung [57], colorectal [58], and thyroid [59]. In these malignancies, downregulation of CBX7 correlated with lower cancer aggressiveness and favorable prognosis, suggesting an oncosuppressor role of CBX7. Conversely, overexpression of CBX7 and its oncogenic effects have been reported in lymphoma [60], prostate [61], and ovarian cancers [62,63].
In addition, CBX8 has been hypothesized to have a contradictory role in various tumors. It was upregulated [64,65], associated with advanced tumor stages, and promoted cancer cell proliferation by inhibiting the p53 pathway in bladder cancer tissues [65]. Similarly, it promoted proliferation in colon cancer, but low expression was associated with poor prognosis of patients [66]. CBX8 has been reported to promote tumorigenesis in esophageal squamous cell carcinoma and to suppress metastasis by repressing Snail [67]. Here, CBX8 mRNA showed a significant increase between tumor and normal conditions in pRCC but not in ccRCC. However, only in this last subtype did CBX8 high expression have a worse prognosis.
Genetic alteration is a common phenomenon in various tumors, including RCC. All eight CBX family members were altered in RCC patients, and the total genetic alteration rate was 76% and 99% for ccRCC and pRCC, respectively. Several clues reveal the involvement of CBX methylation in tumors. For instance, CBX1/3/5 were found to act as methyl readers, which are involved in interpreting H3K9me3 signatures induced by H3K9 methyltransferases [17]. Furthermore, differentially expressed CBXs were found to have a strong association with SUMOylation of DNA methylation proteins in colorectal cancer [36]. In our analysis, CBX3/6/7 resulted in being the genes with the highest alteration rates in both RCC subtypes and their functions primarily related to the SUMOylation of DNA methylation proteins, chromatin organization proteins, RNA binding proteins, and signaling pathways that regulate pluripotency of stem cells, and transcriptional regulation by AKT, PTEN, and RUNX1 pathways. Of note, in agreement with these results, the involvement of the PTEN/Akt pathway in the tumor suppressive activity of CBX7 has been demonstrated in pancreatic cancer cells [68]. Furthermore, it has been reported that this protein regulates stem cell-like properties of gastric cancer cells via the activation of the AKT pathway [69].
RCC is considered an immunogenic tumor, and increasing findings support that immune cell infiltration is closely associated with clinical outcomes in RCC [70][71][72]. In this context, epigenetic regulation of innate immune response is an emerging field. CBX7 knockdown led to apoptosis of CD4 + T cells by hyperdemethylation of the FasL gene promoter and increased expression of FasL [73]. Moreover, experimental evidence has identified Cbx2 as a positive mediator in antiviral immunity, demonstrating its mechanistic involvement as an epigenetic modifier in the innate immune response. Precisely, it has been elucidated that CBX2 can bind to and recruit Jmjd3 to the Ifnb promoter in macrophages, leading to demethylation of H3K27me3 and increasing transcription of IFN-β, an essential factor involved in the innate antiviral immunity [74].
In our study, the expression of CBXs for both RCC subtypes correlated with the infiltration of different immune cells. Particularly, CBXs were associated mainly with neutrophils, CD4 + T cells, and B cells in ccRCC and with dendritic cells, neutrophils, and CD8 + T cells in the papillary subtype. Hence, we hypothesized that CBXs might regulate RCC tumor immunity through multiple immune cell populations. Therefore, in the future, it could be interesting to underpin the role of misregulated CBX1/4/5/7 in ccRCC and CBX4/7 in pRCC in affecting the immune status of RCC. We did not investigate the relationship between CBX members and mast cells. Mast cells are a kind of immune cell that can secrete diverse, active compounds to participate in the immune response. However, recent findings suggested that biological processes such as mast cell activation and mast cell-mediated immunity were regulated by the CBX family members in ovarian cancer [47].
We are aware that our study has several shortcomings. Firstly, the clinicopathological information of the RCC patients was obtained from public databases; secondarily, the study is fully bioinformatics-based and lacks further experimental validation in vitro and in vivo. Furthermore, we lacked research on the detailed mechanisms of individual CBX members in RCC. However, the concordance of expression trends between tissues and cell lines that emerged from in silico validation encourages future functional studies aimed at a deeper understanding of the involvement of this protein family in kidney cancer. Further studies are required to explore the specific mechanism between each CBXs member and kidney cancer.

Conclusions
In conclusion, we analyzed the expression, prognostic values, gene alteration, immune infiltration, and functional protein-association networks of different CBX family members in the two most frequent subtypes of RCC. In ccRCC, as in pRCC, increased expression of CBX3 and lower expression of CBX2 were detected. In addition, decreased expression levels of CBX6/7 were associated with tumor stage for both RCC subtypes. Increased expression values of CBX8 and decreased CBX1/5/6/7 were significantly associated with worse OS in ccRCC patients. Otherwise, in pRCC patients, increased expression levels of CBX4 and decreased CBX7 were significantly associated with shorter OS. Furthermore, a high genetic alteration rate of the CBXs family was observed in RCC. Moreover, individual CBX members were associated with different degrees of immune cell infiltration. Overall, these findings highlighted the interesting value of the CBXs family in the initiation and progression of renal cancer. Our results may provide new cues to select prognostic biomarkers or molecular targets among CBX family members in two major subtypes of RCC.