Stemness-associated senescence genes as potential novel risk factors for papillary renal cell carcinoma.

BACKGROUND: Papillary renal cell carcinoma (PRCC) is the 2nd most common type of renal carcinoma; however, there is limited data about PRCC, and strategies for the diagnosis and treatment of PRCC need to be identified.
METHODS: In this study, the stemness-associated senescence (SAS) phenotype of PRCC was obtained by a bioinformatics analysis. We acquired the gene expression profiles of patients with PRCC and calculated the PRCC messenger ribonucleic acid stemness index (mRNAsi). We then screened the SAS genes from the GenAge database. A least absolute shrinkage and selection operator-Cox regression was conducted to examine correlations between risk signatures and the abundance of the SAS genes in the PRCC samples. Functional enrichment analyses were then performed via molecular co-expression studies of mRNAsi, and the risk scores of PRCC patients were calculated.
RESULTS: We identified the following 8 SAS signatures that were strongly associated with prognosis in PRCC patients: cyclin-dependent kinase 1, heat shock protein family D member 1, platelet-derived growth factor receptor A, cyclin-dependent kinase inhibitor 2B, pyrroline-5-carboxylate reductase 1, sequestosome-1, sirtuin-3, and cyclin-dependent kinase inhibitor 1A. The SAS signatures were significantly associated with the stage and type of PRCC. The calculated risk scores can be used to divide PRCC patients into low- and high-risk groups, and provide guidance in determining treatment plans.
CONCLUSIONS: We have developed a reliable prognostic tool to predict the clinical outcomes of PRCC patients. This tool could improve treatment decisions regarding drug therapy, surgery, and conservative options. 2021 Translational Andrology and Urology. All rights reserved.

Introduction

Renal cell carcinoma (RCC) is the most common cancer affecting the kidney, and 2–3% of adults have this malignancy (1). RCC types can be classified based on their genetic and histological features. Papillary RCC (PRCC) accounts for 10–20% of all RCCs (2); however, data on and evidence about the diagnosis and treatment of PRCC are lacking. The early diagnosis of PRCC is uncommon. Indeed, the majority of patients who are diagnosed with PRCC are already in an advanced stage (3). Surgery is still the primary approach for the treatment of localized diseases; however, up to 40% of surgically resected patients are prone to recurrence or metastasis with a poor prognosis (4).

PRCC is usually classified according to its pathohistological features: a single-layered small cell and scanty cytoplasm are regarded as the signature of type 1 PRCC, while large pseudostratified cells and eosinophilic cytoplasm are regarded as the signature of type 2 PRCC (5). Subtyping is a significant clinical parameter, but the prognosis evaluation value of PRCC morphological subtypes remains controversial. The results of a previous multivariate analysis of 43 PRCC cases revealed that morphological type was associated with PRCC survival (6). However, other studies have reported that the histological subtyping of PRCC has no effect on oncologic outcomes after nephron sparing surgery (7). Thus, prognostic molecular models need to be developed that predict clinical outcomes and contribute to decisions about pharmacotherapy, surgery, and conservative options.

Cancer progression involves transcriptional and epigenetic dysregulation, which frequently leads to oncogenic de-differentiation and the acquisition of stemness features; thus, stemness features indicate greater aggression and resistance to therapy (8). About 6% of RCC cells are regarded as a so-called “side population” that exhibits sphere-forming capacity in vitro and also expresses a set of cancer stem cell markers, including prominin-1 (CD133), paired box gene 2, and Notch Receptor 1. PRCC shares many similarities to proximal tubular progenitor cells, and this type of tumor can develop in a multi-step fashion from benign renal adenomas (9). There were also some researches about the potential novel risk factors for clear cell renal cell carcinoma (ccRCC) (10,11). However, currently, renal carcinoma stem cells are not fully understood, and the putative links between carcinomas and renal stem cells have yet to be elucidated. Focusing on the stemness of PRCC will help to reveal the molecular mechanisms of PRCC initiation and development. Therefore, we used an innovative one-class logistic regression machine learning algorithm (OCLR) to extract transcriptomic and epigenetic feature sets and then obtained the messenger ribonucleic acid stemness index (mRNAsi) for assessing the degree of PRCC stemness based on previous studies (12,13).

Cellular senescence can be described as the stable arrest of the cell cycle coupled with stereotypical phenotypic changes that can contribute to stem cell dysfunction, senescence-related aging phenotypes, chronic diseases, a loss of resilience, and geriatric syndromes (14). Stemness appears to be regulated by senescence via overlapping signaling networks. A previous study showed that the key signaling regulators of the senescence machinery, including p16INK4a, protein 53, p21CIP1 and H3K9me3, operate as critical regulators of stem cell functions (15). Milanovic et al. revealed that senescence evoked cancer cell-intrinsic reprograming into a stem-like state and the acquisition of tumor-initiating potential after their spontaneous escape or forced release from a chemotherapy induced cell-cycle arrest (16). Previous studies have revealed that several gene loci, including tumor suppressor genes, become hypermethylated with age in human and mouse models. The epigenetic defects that accumulate during senescence may specifically affect the function and behavior of stem cells (17). Thus, senescence is not only a fundamental process in PRCC progression, but might provide a mechanistic link between stemness and the malignant progression of PRCC. Based on this, we got aging-related genes from GenAge Database and screened the stemness-associated senescence (SAS) genes by correlation analysis between the expression of aging-related genes in the PRCC samples and their mRNAsi.

We hypothesized that the SAS phenotype of PRCC may confer malignant properties to tumors, resulting in a poor prognosis. Based on this, we sought to clarify the correlation between risk signatures and the abundance of SASs in PRCC.

We present the following article in accordance with the TRIPOD reporting checklist (available at https://dx.doi.org/10.21037/tau-21-913).

Methods

RNA sequence data processing

We obtained the ribonucleic sequencing (RNAseq)-HTSeq-FPKM data of 321 PRCC samples and clinical characteristics of patients with PRCC from the UCSC Xena Database (https://xenabrowser.net/datapages/), which makes comprehensive information and retrieval services accessible from The Cancer Genome Atlas (TCGA) data (18). The clinical features of the PRCC patients are presented in . Two hundred seventy-one PRCC patients’ clinical data and transcriptome profiles were ultimately obtained. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Table 1

Clinical characteristics of PRCC patients

Characteristics	Patients (n=271)
Age
<60	108 (39.85%)
≥60	161 (59.41%)
NA	2 (0.74%)
Gender
Male	199 (73.43%)
Female	72 (26.57%)
BMI
<18.5	5 (1.85%)
18.5–24.9	43 (15.87%)
>24.9	154 (56.83%)
NA	69 (25.46%)
Type
I	72 (26.57%)
II	83 (30.63%)
NA	116 (42.80%)
Disease stage
I	170 (62.73%)
II	19 (7.01%)
III	50 (18.45%)
IV	14 (5.17%)
NA	20 (7.38%)

PRCC, papillary renal cell carcinoma; BMI, body mass index; NA, not available.

Calculation of the stemness index

The messenger ribonucleic acid stemness index (mRNAsi) was obtained from previous studies (12,13). Pluripotent stem cell data were obtained from the Progenitor Cell Biology Consortium (PCBC) data set (https://www.synapse.org/#!Synapse:syn1773109/wiki/54962). We then built a predictive model using a one-class logistic regression (OCLR) algorithm based on the data set. The workflow can be found at https://bioinformaticsfmrp.github.io/PanCanStem_Web/. We used the stem index model to score PRCC samples from TCGA data sets.

Screening SAS genes

Aging-related genes were accessed via the GenAge Database (available at https://genomics.senescence.info/genes/) (19). A correlation analysis between the expression of aging-related genes in the PRCC samples and their mRNAsi was calculated using Pearson’s correlation coefficients; SAS genes were screened using a P value <0.05.

Identification of the prognostic SAS signature

We conducted a univariate Cox regression model to detect prognostic SASs. SASs with P values <0.05 were identified. These SASs were then analyzed in R statistical software (version 3.5.1) using the “glmnet” package (R Foundation for Statistical Computing, Vienna, Austria) and a least absolute shrinkage and selection operator (LASSO)–Cox regression model analysis was conducted (20). A list of prognostic SAS genes with related coefficients was accessed from the expression profile, and PRCC patients’ overall survival (OS) was determined according to the best lambda value.

Risk score of every PRCC patient based on prognostic SAS signatures

The risk score of each PRCC patient was calculated from the data of the expression level of each prognostic SAS and its corresponding coefficient. For prognostic comparisons, OS was analyzed, and we included 321 PRCC samples with risk scores. We used the “surv_cutpoint” function of the R package “survminer” to obtain the best cut-off points for the continuous variables. The PRCC patients were allocated to the high- and low-risk groups using the best cut-off value of the risk score.

Weighted co-expression network construction and functional enrichment analysis associated with PRCC

A co-expression molecular analysis of the mRNAsi and the risk scores of PRCC patients was conducted using the “WGCNA” package in R (21). To investigate the biological functions of the gene modules and key genes, the Cytoscape plug-in ClueGO (v2.5.2) was used to perform the enrichment analysis based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) and WikiPathways databases. The threshold values were defined as follows: FDR <0.05 and P<0.01. In our study, Cytoscape software was used to visualize the significant functional enrichment results.

Statistical analysis

All statistical analyses were conducted using R (version 3.5.1). The predictive accuracy of the SAS risk score were determined by Kaplan-Meier curves. For Kaplan-Meier curves, P values and hazard ratio (HR) with 95% confidence interval (CI) were generated by log-rank tests and univariate Cox proportional hazards regression. When P value was <0.05, the statistical analyses were defined as statistically significant. Pearson’s chi-square tests were executed for the correlation analysis.

Results

Clinical characteristics of PRCC patients

A flowchart detailing the study’s data preparation, processing, and analysis approach is shown in . A total of 271 PRCC patients were enrolled in this study. The detailed clinical characteristics of these patients are summarized in .

Figure 1

The flowchart of our study.

Stemness varies across PRCC and predicts patient survival

We used the stem index model to score the mRNAsi of the PRCC samples (see https://cdn.amegroups.cn/static/public/tau-21-913-1.xlsx). The “Survminer” R package was used to determine the median cut-off point of 0.63 (see ). The PRCC patients were then allocated to the mRNAsi high-risk group (n=236) or the mRNAsi low-risk group (n=35). The results of the Cox regression analysis suggested that patients in the high-risk group had lower survival rates than those in the low-risk group (P=0.002; see ).

Figure 2

Prognostic value of the mRNAsi in PRCC. (A) PRCC patients were allocated to the low-risk mRNAsi group (n=236) or the high-risk mRNAsi group (n=35) based on the median cut-off point of 0.63. (B) Kaplan-Meier OS curves for PRCC patients. PRCC, papillary renal cell carcinoma; OS, overall survival.

Generation of prognostic SAS genes

We accessed aging-related genes from the GenAge Database. Given the relationship between stemness and aging, we then analyzed aging-related gene expression in PRCC samples. The SAS genes are listed in https://cdn.amegroups.cn/static/public/tau-21-913-2.xlsx (P<0.05, n=295). Using a univariate Cox regression model, we identified a set of 30 prognostic SAS genes (see ; P<0.05). A LASSO-Cox regression model was further used to obtain the optimal lambda value from the minimum partial likelihood deviance, which was related to the 8 SAS signatures and was strongly correlated with the prognosis of PRCC patients (see ). The expression levels of the 2 SAS signatures of sirtuin-3 (SIRT3) and cyclin-dependent kinase inhibitor 1A (CDKN1A) were revealed by negative coefficients; these 2 SAS signatures were associated with better survival. Conversely, cyclin-dependent kinase 1 (CDK1), heat shock protein family D member 1 (HSPD1), platelet-derived growth factor receptor A (PDGFRA), cyclin-dependent kinase inhibitor 2B (CDKN2B), pyrroline-5-carboxylate reductase 1 (PYCR1), and Sequestosome 1 (SQSTM1) were negatively correlated with better OS. shows the survival status of all the PRCC patients and the heatmap of the 8 SAS signature prognostic genes.

Figure 3

Development of the prognostic signatures based on the 8 SAS signatures. (A) Univariate Cox regression analyses of SAS signatures associated with OS. (B,C) The LASSO regression identified 8 prognostic SAS signatures. (D) Survival status of patients in different groups (up), distribution of groups based on the signature (middle), heatmap of expression profiles of included SAS (down). LASSO, least absolute shrinkage and selection operator; SAS, stemness-associated senescence; OS, overall survival.

SAS signatures and the patient survival

Based on the expression of these 8 SAS signatures for OS prediction, we established the following formula to calculate the risk score:

Risk score = (0.771443736798139 × CDK1 expression) + (0.296896595577103 × CDKN2B expression) + (0.274905845547291 × PDGFRA expression) + (0.211375135068766 × PYCR1 expression) + (0.0892428222185716 × SQSTM1 expression) + (0.0470091220909639 × HSPD1 expression) + (–0.13579723424102 × SIRT3 expression) + –0.00446166774920266 × CDKN1A expression).

Further, we calculated the SAS signature risk score for every PRCC patient. The patients were partitioned into low- (n=236) and high-risk (n=35) groups based on the median cut-off point of 3.08 (see ). The Cox survival analysis showed that being a high-risk group patient was correlated with a lower survival rate [hazards ratio (HR) =0.0502, 95% confidence interval (CI), 0.0253–0.0997, P=1.21e–17; see ).

Figure 4

Prognostic value of the risk score in PRCC. (A) PRCC patients were allocated to a low-risk score group (n=236) or a high-risk score group (n=35) based on the median cut-off point of 3.08. (B) Kaplan-Meier OS curves for patients allocated to the high- and low-risk score groups. PRCC, papillary renal cell carcinoma.

The “Survival ROC” R package was used to generate time-dependent receiver operating characteristic (ROC) curves. shows that the areas under the curve for 1-, 3-, and 5-year OS were 0.949, 0.839, and 0.802, respectively. Given the above, the SAS prognostic signatures performed well in predicting the OS of PRCC patients.

Figure 5

Survival analysis of the signatures to define the risk groups. Time-dependent ROC analysis of the 8 SAS signatures. SAS, stemness-associated senescence.

SAS signature association with PRCC type and stage

Next, we examined whether the mRNAsi and risk score were affected by PRCC type or stage. In PRCC, stemness varied markedly across type, such that type 1 had a lower risk score than expected, while type 2 had a higher risk score than expected (see ). As for stage, the highest risk score was observed in stage IV and the lowest risk score was observed in stage I (see ). Similar to the risk score, the mRNAsi was positively correlated to PRCC type and tumor stage; thus, similar results were found for the mRNAsi (see ). In summary, we observed positive correlations between stemness or risk score and PRCC type and tumor stage.

Figure 6

SAS signature association with PRCC type and stage. Differences in tumor-type (A,C) or stage (B,D) among the mRNAsi (A,B), or risk score (C,D), and high- or low-risk groups as defined by the SAS prognostic signatures. SAS, stemness-associated senescence; PRCC, papillary renal cell carcinoma.

Weighted gene co-expression network analysis and biological pathway identification

The PRCC samples were clustered using Pearson correlation coefficients and an average linkage (default) algorithm (see ). We conducted a weighted correlation network analysis (WGCNA) to define molecular co-expression across the mRNAsi and risk scores of PRCC patients and construct the co-expression network. The soft-thresholding parameter (β) was 12, and it was selected that resulted in a scale-free R2 fit of 0.9 (see ). Next, 17 modules were identified via average linkage hierarchical clustering—the yellow and green-yellow modules had the highest association with the mRNAsi, while the ink and magenta modules were found to be highly related to risk score (see ). These modules were selected for further analysis.

Figure 7

WGCNA analysis. (A) Analysis of the scale-free fit index for various soft-thresholding parameters (β). (B) Analysis of the mean connectivity for various soft-thresholding powers. (C) Assessment of the scale-free topology when β=12. (D) Heatmap of correlations between clinical traits and module eigengenes. WGCNA, weighted correlation network analysis.

Further, we performed a KEGG enrichment analysis to identify the biological function of the SAS prognostic genes, and found that mRNAsi-related genes were mainly involved in spliceosome, mRNA processing, and the mRNA surveillance pathways (see ). Risk score-related genes were mainly involved in the Ataxia Telangiectasia Mutated Protein (ATM) signaling network in the development of disease, base excision repair, and the cell cycle (see ).

Figure 8

Biological pathway identification. Enriched pathways of the mRNAsi (A,B) or risk score (C,D) associated module genes. The network was generated via the plug-in ClueGO and CluePedia of Cytoscape. The size of the nodes represent the term significance after Bonferroni correction.

Discussion

PRCC has 2 different subtypes: type 1 and type 2. Type 1 is often related to a higher grade, a later stage, poorer outcomes and a more aggressive disease; however, sensitive and specific indicators have yet to be identified (22,23). Some studies have shown that PRCC histological type is not associated with multifocality, and it seems not to be an argument against conservative surgery (24). Histological type was previously identified as a prognostic factor in a univariate analysis but was not retained in the multivariate analysis (24). Thus, the exploration and establishment of sensitive and reliable diagnostic or typing methods for PRCC is of great importance.

Recent advancements in bioinformatics analyses have raised hopes for the investigation of molecular mechanisms and a wide spectrum of disease predictions. In this study, the risk scores were calculated based on the Cox coefficients of the SAS genes. Kaplan-Meier survival analyses and time-dependent ROCs were used to evaluate SAS genes’ prognostic power. Finally, functional enrichment analyses were performed via molecular co-expression of the mRNAsi and risk scores of the PRCC patients.

The following 8 SAS signatures were found to be highly associated with prognosis in PRCC patients: CDK1, HSPD1, PDGFRA, CDKN2B, PYCR1, SQSTM1, SIRT3, and CDKN1A. Among these SAS signatures, CDK1 is upregulated in PRCC patients and represented the largest proportion of the risk score. Previous study had shown CDK1-related genes could effectively predict the overall survival in a TCGA PRCC cohort (n=287, 741.6 versus 1,623.7 days, P<0.0001) (25). CDK1 drives cell division and there is extensive evidence that CDK1 up-regulation promotes cancer proliferation, migration, and invasion. Our results also showed that the risk scores of related genes were involved in the cell cycle. Additionally, CDK1 is also crucial in controlling stemness. Embryonic stem cells (ESCs) display a particularly high level of CDK1 activity that maintains the epigenetic identity of ESCs. CDK1 phosphorylation modifies many epigenetic regulation proteins and the inhibition of CDK1 phosphorylation could change the histone-modification status of ESCs (26). CDK1 maintains the stemness features of bladder cancers, and potentiates the tumorigenesis and aggressive behavior of bladder cancers via a CDK1-dependent mechanism. Moreover, CDKs activity promotes timely onset of senescence (27). CDK1 stimulates pluripotency and regulates the cell-cycle progression of ESCs by the Thr177‐phosphorylation of TFCP2L1 (28). Consistent with the results of our study, Liu et al. showed that CDK1 was significantly associated with a poor prognosis in lung adenocarcinoma (29).

PYCR1 was also identified as a SAS signature that is upregulated in PRCC patients. Feng et al. confirmed PYCR1 expression had a close relationship with the cancer malignant progression in non-small-cell lung carcinoma, lung cancer, and prostate cancer (30,31). A previous study showed that PYCR1 was over-expressed in PRCC tissues and highly regulated PYCR1 expression led to poor prognosis in PRCC patients (32). PYCR1 has been shown to be related to the aging process. Mutations in the PYCR1 gene cause cutis laxa with progeroid features and mental retardation (33). SIRT3 is localized in mitochondria and plays an important role in mitochondrial biology. The knockout of SIRT3 in nude mice led to the development of several diseases related to aging at an accelerated rate, such as cancer, cardiovascular disease, neurodegenerative disease, and metabolic syndrome (34). In our study, we found that high SIRT3 expression levels in PRCC patients were correlated with better survival. Senescence-associated proliferation arrest is characterized by the activity of the TP53-CDKN1A/p21WAF1 and CDKN2A/p16INK4a-RB pathways (35). In relation to SAS signatures, the results of previous studies and the present study were similar in relation to their basic biological functions, which indicates that we calculated and analyzed an accurate prognostic model for PRCC.

However, our study had several limitations. First, our study only extracted mRNA data from TCGA. The 8 SAS prognostic signatures should be further analyzed by experiments to examine the relationship between their protein expression levels and PRCC stage and type. Second, our study provides evidence that the 8 SAS prognostic signatures were significantly correlated to the survival of PRCC patients; however, the analysis was only conducted via a database; thus, further functional experiments of PRCC need to be conducted. Because of the lack of in vitro and in vivo model in PRCC, the experimental data couldn’t be included in the present study.

Conclusions

In conclusion, we identified 8 SAS prognostic signatures and investigated their prognostic value in PRCC patients. We found that SAS signatures are significantly correlated with stage and type in PRCC. These 8 SAS signatures could be expected to be applied for diagnostic biomarkers for PRCC. Furthermore, the risk scores can be calculated to classify PRCC patients into high- and low-risk groups, and provide better reference values to guide therapy than any single conventional clinical parameter.

The article’s supplementary files as