- 1Department of Urology, The First Affiliated Hospital of Anhui Medical University, Hefei, China
- 2Institute of Urology, Anhui Medical University, Hefei, China
- 3Anhui Province Key Laboratory of Genitourinary Diseases, Anhui Medical University, Hefei, China
- 4Department of General Surgery, the First Affiliated Hospital of Nanjing Medical University, Nanjing, China
Clear cell renal cell carcinoma (ccRCC) is the most common subtype of renal cancer and is characterized by high rates of metastasis. Cancer stem cell is a vital cause of renal cancer metastasis and recurrence. However, little is known regarding the change and the roles of stem cells during the development of renal cancer. To clarify this problem, we developed a novel stem cell clustering strategy. Based on The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC) genomic datasets, we used 19 stem cell gene sets to classify each dataset. A machine learning method was used to perform the classification. We classified ccRCC into three subtypes—stem cell activated (SC-A), stem cell dormant (SC-D), and stem cell excluded (SC-E)—based on the expressions of stem cell-related genes. Compared with the other subtypes, C2(SC-A) had the highest degree of cancer stem cell concentration, the highest level of immune cell infiltration, a distinct mutation landscape, and the worst prognosis. Moreover, drug sensitivity analysis revealed that subgroup C2(SC-A) had the highest sensitivity to immunotherapy CTLA-4 blockade and the vascular endothelial growth factor receptor (VEGFR) inhibitor sunitinib. The identification of ccRCC subtypes based on cancer stem cell gene sets demonstrated the heterogeneity of ccRCC and provided a new strategy for its treatment.
Introduction
Renal cancer accounts for 3% of all adult malignancies worldwide, and its incidence have been increasing in recent years. Clear cell renal cell carcinoma (ccRCC) is the most common subtype of renal cancer, jeopardizing 70%–80% of renal cell carcinoma (RCC) patients. Over 30% of ccRCC patients have distant metastasis at the time of diagnosis, and one-third of patients with localized ccRCC will develop metastasis after nephrectomy (1–3). The 5-year survival rate of localized ccRCC is about 65%; however, it drops to 10%–20% after cancer metastasis (4). Surgery remains the major approach for the treatment of localized ccRCC; nevertheless, novel therapeutic strategies are urgently needed for metastatic patients.
Significant achievements have been made in treating advanced ccRCC in the last two decades, such as the application of tyrosine kinase and mammalian target of rapamycin (mTOR) inhibitors or monoclonal antibodies against vascular endothelial growth factor (VEGF) (5, 6). Combination therapy with these inhibitors prolonged the life span of patients. However, most tumors will progress within 2 years. Recently, new approaches for boosting the immune response to renal tumors with immune checkpoint inhibitors, which block programmed cell death protein-1 (PD-1) or cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) on T cells, have shown promising effects in a subset of patients (7). Fundamentally, improving the outcome of renal cancer patients will require personalized treatment strategies specific to the biological characteristic of each tumor.
Stem cells are defined as cells with the ability to self-renew and differentiate into mature cells of a particular tissue (8). Cancer stem cells (CSCs) are a subpopulation of cancer cells with a higher self-renewal ability and the ability to reproduce the heterogeneity of tumors (9). CSCs have been characterized in various cancers and have been proven to contribute to drug resistance, tumor recurrence, and distant metastasis; however, the situation in kidney cancer remains obscure (10–12). Some studies have indicated that targeting the Notch, Hedgehog, and Wnt signaling pathways could inhibit the self-renewal and pluripotency ability of ccRCC cancer stem cells (13, 14). DKK3 and Notch3, which are members of the Wnt and Notch pathways, have been proven to be indicators of the prognosis of renal cancer patients (13). IL8/CXCR1 signaling was proven to promote the sphere formation and self-renewal capability of renal tumor cells (15). Both IL8 and CXCR1 are significantly correlated with patient survival. These studies indicated that the genes expressed in renal CSCs (RCSCs) could be effective prognostic factors. However, the functional significance and the prognostic value of stem cell-related genes in ccRCC are still scarcely investigated and need to be further clarified.
In the present study, we aimed to identify the subclasses of ccRCC with different CSC properties based on the expressions of stem cell-related genes. We divided ccRCC into three clusters with different features and prognosis. In addition, we proved the stability and reliability of this clustering with an independent dataset using the unsupervised clustering method. Moreover, we comprehensively analyzed the prognosis of the RCSC subtypes, relationship with immune cells and genes, the sensitivity of the immune checkpoint inhibitor treatment, and potential changes in the biological process. Classification of the stem cell gene-related subtypes may contribute to formulating the optimal treatment for renal cancer patients.
Results
NMF Identifies Three Subclasses in ccRCC
An analysis flowchart was designed to systematically depict our study (Figure 1). Three hundred and ninety-two stem cell-related genes were enrolled for subsequent non-negative matrix factorization (NMF) analysis. The training set comprising 263 ccRCC samples from The Cancer Genome Atlas (TCGA) was clustered based on the expressions of the aforementioned 392 candidate genes using NMF consensus clustering. Cophenetic correlation coefficients, dispersion, and silhouette were calculated to identify the best k value; k = 3 was proven to be the optimal number of clusters (three subclasses were assigned: C1, C2, and C3) (Figure 2A). Based on the present classification, the consensus heatmap showed sharp and crisp boundaries, indicating the applicability and robust clustering of these samples (Figure 2B). To validate the subtype classification, we conducted t-distributed stochastic neighbor embedding (t-SNE) to reduce the dimension of the features and found that samples in the same subclass generally gather in the same region. This indicates that the subclasses were mostly accordant with the t-SNE distribution patterns (Figure 2C). Additionally, to validate this classification, we performed an independent analysis on TCGA testing set and the International Cancer Genome Consortium (ICGC) dataset, the results of which also demonstrated that there were three distinct molecular subclasses. Based on the classification, a significant prognostic difference was observed in both TCGA testing set and the ICGC dataset (Figures 2D–F).
Figure 2 Identification of three stem cell subtypes using non-negative matrix factorization (NMF) consensus clustering in clear cell renal cell carcinoma (ccRCC) The Cancer Genome Atlas (TCGA) training cohort. (A) NMF clustering using 392 stem cell-associated genes. Cophenetic correlation coefficients for k = 2–6 are shown. (B) Consensus heatmap for ccRCC samples when k = 3. (C) t-Distributed stochastic neighbor embedding (t-SNE) analysis supported the stratification into three stem cell subtypes. Dots with different colors represent different samples in the subclasses. (D–F) Overall survival of the three subclasses (C1, C2, C3) in TCGA training and testing sets and the International Cancer Genome Consortium (ICGC) cohort.
Correlation of the ccRCC Subclass With Stem Cell-Related Signatures
Considering that the clustering was based on stem cell-related genes, we further investigated whether the different subclasses had distinct stem cell characteristics. Firstly, 19 stem cell-related biological process scores were calculated using the GSVA R package based on TCGA training cohort. Three subtypes showed significantly different stem cell-related signatures and clinicopathological characteristics (Figure 3A). Similar trends were identified in the ICGC cohort (Supplementary Figure S1). The results showed that C1 was related to the negative regulation of stem cell maintenance (Figure 3B), while C2 and C3 were correlated with positive stem cell maintenance (Figure 3C). Moreover, the scores of the stem cell proliferation signatures were highest in C2 (Figure 3E). Compared with C3, the C2 subtype was characterized by higher stem cell proliferation and differentiation scores (Figures 3D, E). Hence, we defined C1 as the stem cell excluded (SC_E) subclass, C2 as the stem cell activated (SC_A) subclass, and C3 as the stem cell dormant (SC_D) subclass.
Figure 3 Association between the stem cell-associated signatures and the clear cell renal cell carcinoma (ccRCC) subclasses. (A) Heatmap of the specific stem cell-associated signatures. (B–E) Box plot of the signature scores for the stem cell-associated signatures distinguished by different subclasses. Box plot of the stromal (F) and immune (G) scores from estimates of the three subclasses. The p-values are labeled above each box plot with asterisks. ns, no significance. ***p < 0.001, ****p < 0.0001.
To further characterize the subclasses, the expressions of the RCSC surface markers (CD44, CXCRL8, CXCR4, and ENG) and the stem cell-related pathways such as hypoxia-inducible factor (HIF), Notch, Wnt, and Hedgehog were investigated. Subclass C2(SC_A) had the highest expression of RCSC marker genes (Figure 4E). Moreover, subclasses C2(SC_A) and C3(SC_D) had higher normal stem cell and CSC gene markers (Figures 4A, C). Generally, compared to that in subclasses C2(SC_A) and C3(SC_D), HIF, Notch, Wnt, and Hedgehog were less activated in subclass C1(SC_E) (Figures 4B, D, F, G). Besides, the immune and stromal scores of the subclasses were calculated using the ESTIMATE (estimation of stromal and immune cells in malignant tumours using expression data) algorithm. As shown in Figures 3F, G, subclass C2(SC_A) had the highest stromal and immune cell infiltration scores.
Figure 4 Expressions of clear cell renal cell carcinoma (ccRCC) stem cell-related genes and pathways. (A, C, E) Expression levels of cancer stem cell marker genes (A), normal stem cell marker genes (C), and ccRCC cancer stem cell marker genes (E). (B, D, F, G) Box plot of the status of ccRCC cancer stem cell-related pathways. The p-values are labeled above each box plot with asterisks. ns, no significance. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
Correlation of the ccRCC Subclasses With Immune Infiltration
Considering the significant differences in the immune and stromal scores displayed among the subclasses, immune and stromal cell infiltration was explored to depict their microenvironment landscape based on TCGA training cohort. Using the single-sample gene set enrichment analysis (ssGSEA) algorithm, the abundance rates of the 28 immune-related cell types were quantified and presented in a heatmap (Figure 5A). Subclass C2(SC_A) showed a significantly different immune cell infiltration compared with the other two subclasses. Stromal cells, especially fibroblasts, which have an important role in renal cancer progression, showed different infiltration status in the three subtypes (Figures 5B, C). Moreover, we explored the association between subclasses and the expressions of potentially targetable immune checkpoint genes, which have been used for drug inhibitors in clinical trials or approved for specific cancer treatment (Figure 5D). Subclass C2(SC_A) had the highest expressions of most of these genes.
Figure 5 Immune characteristics of the three subclasses in The Cancer Genome Atlas (TCGA) training set. (A) Heatmap describing the abundance of immune cell populations in C1, C2, and C3. Box plot of the abundance of endothelial (B) and fibroblasts (C) distinguished by different subclasses. (D) Expression levels (in fragments per kilobase of transcript per million mapped reads, FPKM) of 17 immune checkpoint genes in the three clear cell renal cell carcinoma (ccRCC) subclasses. ns, no significance. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
Correlation of the ccRCC Subclasses With Mutations and Copy Number Variations
The tumor genomic landscape has been demonstrated to be associated with antitumor immunity. To investigate whether differences exist in the somatic mutation frequencies across the ccRCC subclasses and to observe the different patterns of mutations among the ccRCC clusters, somatic mutation data from TCGA database were analyzed. Figures 6A–C show the landscape of the top 20 mutated genes in the three subtypes. There was no significant difference in the mutation count of the three subclasses (Figure 6E). However, the tumor mutation burden (TMB) (Figure 6D) and the fraction genome altered (Figure 6F) were significantly higher in subclass C1(SC_E) compared to those in subclasses C2(SC_A) and C3(SC_D).
Figure 6 Association between the clear cell renal cell carcinoma (ccRCC) subclasses and mutations. (A–C) Oncoprint of the mutation status of the top 20 genes in subclasses C1(SC_E), C2(SC_A), C3(SC_D). (D–F) The tumor mutation burden (D), mutation count (E), and fraction genome altered (F) in the three subclasses. ns, no significance. *p < 0.05, **p < 0.01, ***p < 0.001.
Transcriptome Features of the ccRCC Subclasses
To better characterize the three ccRCC subclasses, gene differential analysis was conducted. Genes with an adjusted p-value less than 0.01 and an absolute log2 fold change larger than 1 were considered significantly differential. Only those genes which showed significant differences in two possible comparisons were considered subclass-specific genes. Finally, all of the 1,695 subclass-specific genes were identified, including 172 specific genes for C1(SC_E), 1,485 specific genes for C2(SC_A), and 38 specific genes for C3(SC_D). Subsequently, gene ontology (GO) enrichment analysis of the subclass-specific genes was performed using clusterProfiler in the R package. The significantly enriched biological processes are shown in Figure 7. Subclass C2(SC_A), which was enriched in extracellular matrix organization and extracellular structure organization, had significantly different enriched pathways compared with the other subclasses (Figure 7B).
Figure 7 Enrichment analysis of the differentially expressed genes in three different subclasses: C1 (A), C2 (B), and C3 (C).
Prediction of the Therapeutic Response of the ccRCC Stem Cell Subtypes to Immune Checkpoint Inhibitors and Target Therapy
Based on the above results, we further evaluated the response of the three subtypes to immunotherapy. In RCC, the blockade of PD-1 and CTLA-4 has become the new treatment approach in patients with intermediate- and high-risk metastatic tumors, whereas monotherapy with the PD-1 inhibitor nivolumab is the second-line or third-line treatment approach after failure of VEGF tyrosine kinase inhibitors (16). In 2019, the United States Food and Drug Administration (FDA) approved the combination use of PD-1 blockade and anti-angiogenic therapy for the treatment of patients with advanced RCC (17). To predict the sensitivity to immunotherapy of the different clusters, we performed subclass mapping to compare the expression profiles of the three stem cell subtypes with 47 melanoma patients who were treated with immunotherapy (18). The subclass mapping results indicated that the C2(SC_A) subtype might be more sensitive to anti-CTLA-4 treatments (Figure 8A).
Figure 8 Immunotherapy and target therapy response prediction in the clear cell renal cell carcinoma (ccRCC) subclasses. (A) Response of the three ccRCC subclasses to PD-1 and CTLA-4 immunotherapy. (B–F) Sensitivity of the three ccRCC subclasses to sorafenib (B), sunitinib (C), rapamycin (D), pazopanib (E), and axitinib (F). ****p < 0.0001.
Since VEGF receptor (VEGFR) target therapy is a more conventional therapy for patients with advanced ccRCC, we selected five conventional target therapy agents (sunitinib, sorafenib, axitinib, pazopanib, and rapamycin) and evaluated the responses of the three subtypes. We constructed a prediction model on the Genomics of Drug Sensitivity in Cancer (GDSC) cell line dataset using ridge regression and evaluated the prediction accuracy using 10-fold cross-validation. We estimated the half-maximal inhibitory concentration (IC50) of each sample in the training set based on the prediction models for the four target therapy agents. Regarding sunitinib and pazopanib, subclass C2(SC_A) was the most sensitive (Figures 8C, E), while for sorafenib, subclass C3(SC_D) had the worst sensitivity. Subclasses C1(SC_E) and C2(SC_A) showed similar sensitivity values (Figure 8B). Regarding axitinib, subclass C3(SC_D) was the most sensitive (Figure 8F), while for rapamycin, subclass C1(SC_E) was the most sensitive (Figure 8D).
Construction of a Prognostic Model Based on Key Genes
To better characterize the prognosis of each patient, we constructed a risk model based on the key genes. The key genes were extracted from the stem cell-related genes using the NMF package in R. Twenty key genes were identified. Then, we applied least absolute shrinkage and selection operator (LASSO) Cox regression analysis to select the most useful predictive features and identified four genes (A2M, CFL1, FN1, and PSME2) with non-zero regression coefficients (Figures 9A, B). Ultimately, a six-gene risk signature was built, and the risk score of each patient was calculated using the following formula:
Figure 9 Construction of a risk prediction model based on the key genes in the three subclasses. (A) Tuning parameter (λ) screening in the least absolute shrinkage and selection operator (LASSO) regression model. (B) LASSO coefficient profiles of the common genes. (C) From top to bottom are the risk score distribution, survival overview, and heatmap analysis of six genes. (D–F) Kaplan–Meier plots of overall survival (OS) according to the risk scores in The Cancer Genome Atlas (TCGA) training set (D), TCGA testing set (E), and the International Cancer Genome Consortium (ICGC) cohort (F).
Signature risk core = (−0.00205964493004206 × A2M expression) + (0.00157204565454328 × CFL1 expression) + (0.00199772913455084 × FN1 expression) + (0.0181057563160569 × PSME2 expression). Increased expressions of PSME2, FN1, and CFL1 correlated with higher risk scores and worse survival outcomes (Figure 9C). More importantly, the risk score could stratify patients into a high- and a low-risk group with significantly different survival outcomes (Figure 9D). Similar results were found in TCGA testing set and the ICGC dataset (Figures 9E, F).
Discussion
Although a variety of ccRCC classifications based on gene expression have been developed in recent years, a consensus in molecular subtype has not yet been achieved. To identify the ccRCC subgroups associated with CSC and patient prognosis, a ccRCC classification was developed in this study based on 392 genes retrieved from Molecular Signatures Database. Three subclasses of ccRCC with different prognosis were identified. Subsequently, the stem cell signature, immune infiltration, mutation landscape, and clinicopathological characteristics of the subclasses and their sensitivity to immunotherapy were investigated. The results showed that three subclasses were distinct, with significantly different stem cell and immune cell infiltration signatures. Drug sensitivity analysis demonstrated that subclass C2(SC_A) was sensitive toward CTLA-4 inhibitors and sunitinib. In addition, based on the marker genes of each cluster, we constructed a risk model to predict the prognosis of patients. This risk model could stratify patients with different prognosis, and it was validated in an external cohort.
The progression of cancer was accompanied by the gradual loss of differentiation ability and the gain of stem cell-like characteristics (12). Inhibiting the self-renewal capacity and the tumorigenicity of ccRCC significantly suppressed tumor growth and metastasis (19, 20). Given that a relapse in ccRCC has been attributed to the maintenance of ccRCC stemness cells, possessing stem cell properties that lead to therapy resistance (21), there is an urgent need for the development of prognostic biomarkers associated with stem cell properties. Malta et al. developed a novel transcriptome stemness index called mRNAsi (mRNA expression-based stemness index) to evaluate the stemness based on the one-class logistic regression machine learning algorithm (22). However, the mRNAsi was higher in normal renal tissue compared with that in renal tumor and showed no correlation with the survival outcomes of patients. This contradicts usual biological experiment results that demonstrate CSC properties to indicate worse prognosis. A novel stem cell-related prognostic evaluation model is needed.
Based on stem cell-related genes, three subtypes with different survival outcomes were identified. The results showed that subclass C1(SC_E) was distinct with the negative regulation of stem cell maintenance; however, subclasses C2(SC_A) and C3(SC_D) displayed distinct stem cell signatures and were characterized by positive regulation of stem cell maintenance. Moreover, subclass C2(SC_A) was also highly correlated with stem cell division and differentiation. The signature differences between subclasses C2(SC_A) and C3(SC_D) may be caused by the dormancy of CSCs (23, 24). Renal CSC surface markers, such as CD44, CXCR4, and CD105(ENG), were highly expressed in subclasses C2(SC_A) and C3(SC_D) (10, 11). BMP2 and CXCL8, which could promote the self-renewal of RCSCs, also had higher expressions in subclasses C2(SC_A) and C3(SC_D) (15, 25). The HIF, Notch, Wnt, and Hedgehog signaling pathways were reported to promote renal cancer progression by regulating the self-renewal and stemness maintenance of RCSCs (14, 26–28). Hence, we also investigated their status in the different subclasses. Compared to that in subclass C1(SC_E), these pathways were significantly activated in subclasses C2(SC_A) and C3(SC_D). These results demonstrated the validity of this classification. Moreover, these signaling pathways were also differently expressed in subclasses C2(SC_A) and C3(SC_D). These results indicated that these pathways were not only implicated in stemness maintenance but also associated with the proliferation and differentiation of RCSCs.
A recent study has demonstrated that stem cell properties are microenvironment defined during tumor progression (29, 30). Hence, we analyzed the stromal and immune infiltration levels in the three subtypes using R package ESTIMATE. The C2(SC_A) subtype had higher stromal and immune scores. Subsequent analysis further corroborated the C2(SC_A) subtype as possessing distinct stromal and immune features, including high infiltration of fibroblasts, T cells, and macrophages. Heterogeneous stromal cells in the tumor microenvironment can profoundly boost cancer progression (31). Carcinoma-associated fibroblasts (CAFs) form the chief components of the tumor microenvironment in multiple types of malignancies (32). By providing a supporting niche for CSCs, CAFs could facilitate tumor formation and induce chemoresistance (31). This may partly be explained by the concurrence of a high stem cell maintenance score and fibroblast infiltration. Macrophages abundantly exist in the immune milieu, where they share the microenvironment with CSCs. Macrophage-initiated WNT signaling could contribute to the maintenance of stemness, leading to the characteristics of chemoresistance and invasiveness in ovarian cancer (33). A similar phenomenon was found in lung cancer. A positive feedback interaction between macrophages and cancer cells could promote the stemness of cancer cells (34). Tumor-associated macrophages have been demonstrated to contribute to the maintenance of breast CSC populations through triggering the production of the inflammatory cytokines interleukin 1 (IL-1), IL-6, and IL-8, which, in turn, reinforce the CSC states (35). Interestingly, previous studies proved that IL-8 could boost the CSC-like properties of ccRCC (15). Besides, the activation of the Notch signaling pathway could promote the expansion of ccRCC-derived CSCs and induce chemotaxis simultaneously (27). Th17 cells are another type of immune cells that appear to support CSCs. Th17 cell-associated cytokines could transform dormant stem cells into an active state (36). This is consistent with our finding that subclass C2(SC_A) had a higher infiltration of Th17 helper cells. Tumor-specific antigens are usually generated by somatic mutation and can influence the response of patients to immunotherapy (37, 38). Hence, we comprehensively analyzed the mutation status of the three subclasses. Although there was no significant difference in the TMB among the three subtypes, C1(SC_E) had higher mutation counts and fraction genome altered than the other subclasses. These differences might influence their response to immunotherapy. VHL mutation is the most common mutation in ccRCC (39). However, evidence of the relationship between the mutation status of the VHL gene and ccRCC remains few and contradictory. Some studies suggested that VHL mutation could activate effector T cells and promote the secretion of cytokines (40). However, a recent study has proposed that wild-type VHL positively correlated with the expression of programmed death-ligand 1 (PD-L1) (41). In the present study, we found that subclass C2(SC_A) had a relatively lower VHL mutation frequency. PBRM1 is another commonly mutated gene in ccRCC. A previous study indicated that CD8+ T-cell-infiltrated tumors had relatively fewer PBRM1 mutations (42). Similarly, we found that the mutation frequency of PBRM1 was significantly lower in subclass C2(SC_A), which had the highest immune cell infiltration. Moreover, the mutation landscape showed that subclass C1 had the highest mTOR mutation frequency. Correspondingly, patients in subclass C1 were most sensitive to the mTOR inhibitor rapamycin. This demonstrated the reliability of the present study. More studies are needed to investigate the relationship between somatic mutations and immune infiltration.
Considering the close interactions between CSCs and the immune system, developing therapeutic strategies that target immune checkpoints might pave the way to eradicating CSCs. At present, immunotherapy has obtained global attention in cancer management. The efficacy and safety of PD-1 immune checkpoint inhibitors and CTLA-4 inhibitors have been applied clinically and have shown promising outcomes (43, 44). Due to the interaction between PD-1 and PD-L1 in tumor-infiltrating lymphocytes and tumor cells, T-cell exhaustion, tumor-specific T-cell dysfunction, and immune evasion by tumor cells were triggered. Exhausted T cells could produce additional inhibitory molecules to promote the progression of cancer; however, this process could be reversed by a combined PD-1 and CTLA-4 blockade. Treating a mouse tumor model with PD-L1 and CTLA-4 inhibitors could promote the elimination of CSCs (45).
Redirecting immune suppression by targeting checkpoints has brought about clinical response in some RCC patients, and a combination treatment involving checkpoint blockade is now the standard of therapy in advanced RCC patients. However, a substantial subset of patients is not sensitive to checkpoint blockade. The identification of a reliable evaluation system to predict the response to checkpoint blockade is essential to improve the clinical efficacy of these therapies (46). Moreover, the immune checkpoint gene CXCR4, which is also a marker of RCSC, was highly expressed in subclass C2(SC_A). Targeting the tumor microenvironment may provide a promising therapeutic avenue as the eliminated CSCs could be replenished by non-CSCs for the existence of a survival niche (47). The highest expressions of immune checkpoint genes indicated the sensitivity of subclass C2(SC_A) to immunotherapy. The results demonstrated that anti-CTLA-4 therapy and sunitinib are promising for patients in subclass C2(SC_A). The results of our study provide a novel insight into the combination of anti-CTLA-4 therapy and sunitinib, which requires further validation in future research.
Besides, to better characterize the prognosis of each patient, we constructed a risk model based on the key genes (A2M, CFL1, FN1, and PSME2). The risk scores could classify patients into the high- and low-risk groups with significantly different survival outcomes. In addition, the effectiveness of this risk model was validated in TCGA testing set and the ICGC dataset.
In summary, in this study, we explored the stem cell-related process landscape of ccRCC and identified three subclasses with different stem cell activities. We systematically analyzed the differences of these subclasses in the tumor microenvironment, immune cell infiltration, immunotherapy/target therapy response, and the corresponding pathways and constructed a risk model. The results of this study provide a basis and reference for the treatment and prognostic prediction of ccRCC.
With the development of target therapy and immunotherapy, various drugs have been used for the clinical treatment of advanced ccRCC. The VEGFR inhibitor, anti-checkpoint therapy, and the combination of both showed promising efficacy (48, 49). However, few references were established for the selection of a proper treatment plan. The novel developed stem cell subtypes are critical for selecting suitable therapies for ccRCC patients. Patients in the stem cell activated subtype (SC-A) could benefit more from anti-CTLA-4 and sunitinib treatment. For those in the stem cell excluded subtype, an anti-PD-1 therapy might be more suitable. However, there are limitations in the present study. This study was conducted based on an existing public dataset. These findings need to be validated in larger ccRCC patient cohorts with immunotherapy and target therapy experience.
Materials and Methods
Patients and Samples
The RNA sequencing data (raw counts) of 530 and 91 ccRCC patient samples with corresponding clinical information were download from The Cancer Genome Atlas (TCGA; http://cancergenome.nih.gov/) and the International Cancer Genome Consortium (ICGC; www.icgc.org), respectively. Patients with 0 day overall survival (OS) were removed; 526 TCGA samples were retained for further analysis. Subsequently, the dataset from TCGA was randomly divided into a training set and a testing set. The gene expression data of 91 ccRCC samples from the ICGC were used for external validation. In total, 662 ccRCC patients were enrolled in the present study. The gene somatic mutation data (MAF files) of ccRCC were retrieved from TCGA.
Identification of ccRCC Subclasses
Human stem cell-related biological processes were downloaded from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). A total of 392 genes from 19 stem cell-related biological processes were obtained (Supplementary Table S1). These genes were FPKM (fragments per kilobase of transcript per million mapped reads) normalized and used for subsequent non-negative matrix factorization (NMF) analysis using the R NMF package in the training set. NMF is an unsupervised learning technique that has been used to extract meaningful information from high-dimensional data (50). In detail, the best k-value was defined according to the cophenetic correlation coefficients, dispersion, and silhouette. The iteration was set as nrun = 100. This method was also applied to the testing and external validation sets using the same candidate genes. The first value of k for which the cophenetic coefficient starts decreasing was chosen as the optimal number for clustering. A t-distributed stochastic neighbor embedding (t-SNE)-based approach was then used to validate the sample clustering using the mRNA expression data of the above stem cell-related genes.
Gene Set Variation Analysis
Gene set variation analysis (GSVA) is a non-parametric and unsupervised gene set enrichment method that can estimate the scores of certain pathways or signatures based on transcriptomic data. Using the GSVA R package, each sample received 19 scores corresponding to 19 stem cell-related signatures. Subsequently, differences in the signatures of the different clusters were calculated using the t-test in R.
Estimation of Immune Infiltration
The microenvironment cell populations counter (MCP-counter), a method based on transcriptomic data, was used to assess the absolute abundance of two stromal cell populations (endothelial cells and fibroblasts) (51). Furthermore, another approach applied for the qualification of the immune infiltration of 28 immune cells used in this research was single-sample gene set enrichment analysis (ssGSEA), which calculated an enrichment score (ES) representing the variations of the pathway activities within a single sample (51). In addition, immune scores and stromal scores were calculated using the ESTIMATE algorithm, which can reflect the infiltration level of stromal and immune cells.
Enrichment Analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and Gene Ontology (GO) analysis were performed with the R package “clusterProfiler.”
Target Therapy and Immunotherapy Sensitivity Prediction
The sensitivity of patients to target therapy drugs was evaluated based on the GDSC database (https://www.cancerrxgene.org/) (52). IC50 values were estimated using the R package pRRophetic (53). In detail, the IC50 was calculated by ridge regression and the prediction accuracy was evaluated using 10-fold cross-validation based on the GDSC training set. The response to anti-PD-1 and anti-CTLA-4 therapy was predicted by comparing the expression profiles of three subtypes with 47 melanoma patients who respond to the immunotherapy using subclass mapping (https://www.genepattern.org/).
Construction of a Prognostic Model
Univariate Cox regression was used to screen the mRNAs affecting the OS of patients (p < 0.05). Thereafter, survival-related genes were screened with the LASSO multivariate Cox regression algorithm using the R package “glmnet” (version 3.0). Finally, the signature genes and coefficients in the risk score signature were constructed based on the most proper penalty parameter λ. The risk score formula used was:
where Coefi is the coefficient and Expi is the normalized expression of each gene in the signature. The risk score system was constructed using the training set and evaluated in the testing set. Patients were stratified into a high-risk group and a low-risk group based on the median risk score.
Statistical Analysis
All statistical analysis was performed using R programming (https://www.r-project.org/). Unpaired Student’s t-test was used to compare two groups with normally distributed variables. For the comparison of three groups, one-way analysis of variance and the Kruskal–Wallis test were used as parametric and non-parametric methods, respectively. Contingency table variables were analyzed with the chi-square test or Fisher’s exact tests. Survival analysis was carried out using the Kaplan–Meier method and compared using the log-rank test. A p-value less than 0.05 was considered as statistically significant.
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.
Author Contributions
HW designed the study and drafted the manuscript. HW and HX collected, analyzed, and interpreted the data. QC and CL participated in revising the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by Research Fund of Anhui Institute of Translational Medicine (grant number: ZHYX2020A003).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We thank Dr. Jianming Zeng (University of Macau) and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.758989/full#supplementary-material
Abbreviations
ccRCC, clear cell renal cell carcinoma; CSC, cancer stem cell; SC-A, stem cell activated; SC-D, stem cell dormant; SC-E: stem cell excluded; OS: overall survival; TMB: tumor mutation burden.
References
1. Hoffman AM, Cairns P. Epigenetics of Kidney Cancer and Bladder Cancer. Epigenomics (2011) 3:19–34. doi: 10.2217/epi.10.64
2. Linehan WM, Rathmell WK. Kidney Cancer. Urologic Oncol (2012) 30:948–51. doi: 10.1016/j.urolonc.2012.08.021
3. Sudarshan S, Karam JA, Brugarolas J, Thompson RH, Uzzo R, Rini B, et al. Metabolism of Kidney Cancer: From the Lab to Clinical Practice. Eur Urol (2013) 63:244–51. doi: 10.1016/j.eururo.2012.09.054
4. Znaor A, Lortet-Tieulent J, Laversanne M, Jemal A, Bray F. International Variations and Trends in Renal Cell Carcinoma Incidence and Mortality. Eur Urol (2015) 67:519–30. doi: 10.1016/j.eururo.2014.10.002
6. Busch J, Seidel C, Goranova I, Erber B, Peters R, Friedersdorff F, et al. Categories of Response to First Line Vascular Endothelial Growth Factor Receptor Targeted Therapy and Overall Survival in Patients With Metastatic Renal Cell Carcinoma. Eur J Cancer (2014) 50:563–9. doi: 10.1016/j.ejca.2013.10.017
7. Motzer RJ, Escudier B, McDermott DF, George S, Hammers HJ, Srinivas S, et al. Nivolumab Versus Everolimus in Advanced Renal-Cell Carcinoma. N Engl J Med (2015) 373:1803–13. doi: 10.1056/NEJMoa1510665
8. Reya T, Morrison SJ, Clarke MF, Weissman IL. Stem Cells, Cancer, and Cancer Stem Cells. Nature (2001) 414:105–11. doi: 10.1038/35102167
9. Nassar D, Blanpain C. Cancer Stem Cells: Basic Concepts and Therapeutic Implications. Annu Rev Pathol (2016) 11:47–76. doi: 10.1146/annurev-pathol-012615-044438
10. Bussolati B, Bruno S, Grange C, Ferrando U, Camussi G. Identification of a Tumor-Initiating Stem Cell Population in Human Renal Carcinomas. FASEB J (2008) 22:3696–705. doi: 10.1096/fj.08-102590
11. Gassenmaier M, Chen D, Buchner A, Henkel L, Schiemann M, Mack B, et al. CXC Chemokine Receptor 4 Is Essential for Maintenance of Renal Cell Carcinoma-Initiating Cells and Predicts Metastasis. Stem Cells (2013) 31:1467–76. doi: 10.1002/stem.1407
12. Grange C, Tapparo M, Collino F, Vitillo L, Damasco C, Deregibus MC, et al. Microvesicles Released From Human Renal Cancer Stem Cells Stimulate Angiogenesis and Formation of Lung Premetastatic Niche. Cancer Res (2011) 71:5346–56. doi: 10.1158/0008-5472.Can-11-0241
13. Fendler A, Bauer D, Busch J, Jung K, Wulf-Goldenberg A, Kunz S, et al. Inhibiting WNT and NOTCH in Renal Cancer Stem Cells and the Implications for Human Patients. Nat Commun (2020) 11:929. doi: 10.1038/s41467-020-14700-7
14. Qu L, Wu Z, Li Y, Xu Z, Liu B, Liu F, et al. A Feed-Forward Loop Between lncARSR and YAP Activity Promotes Expansion of Renal Tumour-Initiating Cells. Nat Commun (2016) 7:12692. doi: 10.1038/ncomms12692
15. Corrò C, Healy ME, Engler S, Bodenmiller B, Li Z, Schraml P, et al. IL-8 and CXCR1 Expression is Associated With Cancer Stem Cell-Like Properties of Clear Cell Renal Cancer. J Pathol (2019) 248:377–89. doi: 10.1002/path.5267
16. Amin A, Hammers H. The Evolving Landscape of Immunotherapy-Based Combinations for Frontline Treatment of Advanced Renal Cell Carcinoma. Front Immunol (2018) 9:3120. doi: 10.3389/fimmu.2018.03120
17. Escudier B. Combination Therapy as First-Line Treatment in Metastatic Renal-Cell Carcinoma. New Engl J Med (2019) 380:1176–8. doi: 10.1056/NEJMe1900887
18. Roh W, Chen PL, Reuben A, Spencer CN, Prieto PA, Miller JP, et al. Integrated Molecular Analysis of Tumor Biopsies on Sequential CTLA-4 and PD-1 Blockade Reveals Markers of Response and Resistance. Sci Transl Med (2017) 9(379):eaah3560. doi: 10.1126/scitranslmed.aah3560
19. Ge Y, Weygant N, Qu D, May R, Berry WL, Yao J, et al. Alternative Splice Variants of DCLK1 Mark Cancer Stem Cells, Promote Self-Renewal and Drug-Resistance, and can be Targeted to Inhibit Tumorigenesis in Kidney Cancer. Int J Cancer (2018) 143:1162–75. doi: 10.1002/ijc.31400
20. Lawson DA, Bhakta NR, Kessenbrock K, Prummel KD, Yu Y, Takai K, et al. Single-Cell Analysis Reveals a Stem-Cell Program in Human Metastatic Breast Cancer Cells. Nature (2015) 526:131–5. doi: 10.1038/nature15260
21. Hu J, Guan W, Liu P, Dai J, Tang K, Xiao H, et al. Endoglin Is Essential for the Maintenance of Self-Renewal and Chemoresistance in Renal Cancer Stem Cells. Stem Cell Rep (2017) 9:464–77. doi: 10.1016/j.stemcr.2017.07.009
22. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine Learning Identifies Stemness Features Associated With Oncogenic Dedifferentiation. Cell (2018) 173:338–54.e315. doi: 10.1016/j.cell.2018.03.034
23. Talukdar S, Bhoopathi P, Emdad L, Das S, Sarkar D, Fisher PB. Dormancy and Cancer Stem Cells: An Enigma for Cancer Therapeutic Targeting. Adv Cancer Res (2019) 141:43–84. doi: 10.1016/bs.acr.2018.12.002
24. Triana-Martínez F, Loza MI, Domínguez E. Beyond Tumor Suppression: Senescence in Cancer Stemness and Tumor Dormancy. Cells (2020) 9(2):346. doi: 10.3390/cells9020346
25. Wang L, Park P, La Marca F, Than KD, Lin CY. BMP-2 Inhibits Tumor-Initiating Ability in Human Renal Cancer Stem Cells and Induces Bone Formation. J Cancer Res Clin Oncol (2015) 141:1013–24. doi: 10.1007/s00432-014-1883-0
26. Myszczyszyn A, Czarnecka AM, Matak D, Szymanski L, Lian F, Kornakiewicz A, et al. The Role of Hypoxia and Cancer Stem Cells in Renal Cell Carcinoma Pathogenesis. Stem Cell Rev Rep (2015) 11:919–43. doi: 10.1007/s12015-015-9611-y
27. Xiao W, Gao Z, Duan Y, Yuan W, Ke Y. Notch Signaling Plays a Crucial Role in Cancer Stem-Like Cells Maintaining Stemness and Mediating Chemotaxis in Renal Cell Carcinoma. J Exp Clin Cancer Res (2017) 36:41. doi: 10.1186/s13046-017-0507-3
28. Lu J, Wei JH, Feng ZH, Chen ZH, Wang YQ, Huang Y, et al. miR-106b-5p Promotes Renal Cell Carcinoma Aggressiveness and Stem-Cell-Like Phenotype by Activating Wnt/β-Catenin Signalling. Oncotarget (2017) 8:21461–71. doi: 10.18632/oncotarget.15591
29. Lenos KJ, Miedema DM, Lodestijn SC, Nijman LE, van den Bosch T, Romero Ros X, et al. Stem Cell Functionality is Microenvironmentally Defined During Tumour Expansion and Therapy Response in Colon Cancer. Nat Cell Biol (2018) 20:1193–202. doi: 10.1038/s41556-018-0179-z
30. Chang WH, Lai AG. Aberrations in Notch-Hedgehog Signalling Reveal Cancer Stem Cells Harbouring Conserved Oncogenic Properties Associated With Hypoxia and Immunoevasion. Br J Cancer (2019) 121:666–78. doi: 10.1038/s41416-019-0572-9
31. Su S, Chen J, Yao H, Liu J, Yu S, Lao L, et al. CD10(+)GPR77(+) Cancer-Associated Fibroblasts Promote Cancer Formation and Chemoresistance by Sustaining Cancer Stemness. Cell (2018) 172:841–856.e816. doi: 10.1016/j.cell.2018.01.009
32. Kalluri R. The Biology and Function of Fibroblasts in Cancer. Nature Reviews. Cancer (2016) 16:582–98. doi: 10.1038/nrc.2016.73
33. Raghavan S, Mehta P, Xie Y, Lei YL, Mehta G. Ovarian Cancer Stem Cells and Macrophages Reciprocally Interact Through the WNT Pathway to Promote Pro-Tumoral and Malignant Phenotypes in 3D Engineered Microenvironments. J Immunother Cancer (2019) 7:190. doi: 10.1186/s40425-019-0666-1
34. Lu CH, Yeh DW, Lai CY, Liu YL, Huang LR, Lee AY, et al. USP17 Mediates Macrophage-Promoted Inflammation and Stemness in Lung Cancer Cells by Regulating TRAF2/TRAF3 Complex Formation. Oncogene (2018) 37:6327–40. doi: 10.1038/s41388-018-0411-0
35. Lu H, Clauser KR, Tam WL, Fröse J, Ye X, Eaton EN, et al. A Breast Cancer Stem Cell Niche Supported by Juxtacrine Signalling From Monocytes and Macrophages. Nat Cell Biol (2014) 16:1105–17. doi: 10.1038/ncb3041
36. Shahid A, Bharadwaj M. The Connection Between the Th17 Cell Related Cytokines and Cancer Stem Cells in Cancer: Novel Therapeutic Targets. Immunol Lett (2019) 213:9–20. doi: 10.1016/j.imlet.2019.07.001
37. Patel SJ, Sanjana NE, Kishton RJ, Eidizadeh A, Vodnala SK, Cam M, et al. Identification of Essential Genes for Cancer Immunotherapy. Nature (2017) 548:537–42. doi: 10.1038/nature23477
38. Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A, et al. Development of Tumor Mutation Burden as an Immunotherapy Biomarker: Utility for the Oncology Clinic. Ann Oncol (2019) 30:44–56. doi: 10.1093/annonc/mdy495
39. Zhang J, Wu T, Simon J, Takada M, Saito R, Fan C, et al. VHL Substrate Transcription Factor ZHX2 as an Oncogenic Driver in Clear Cell Renal Cell Carcinoma. Science (2018) 361:290–5. doi: 10.1126/science.aap8411
40. Hsieh JJ, Le VH, Oyama T, Ricketts CJ, Ho TH, Cheng EH. Chromosome 3p Loss-Orchestrated VHL, HIF, and Epigenetic Deregulation in Clear Cell Renal Cell Carcinoma. J Clin Oncol (2018) 36:Jco2018792549. doi: 10.1200/jco.2018.79.2549
41. Kammerer-Jacquet SF, Crouzet L, Brunot A, Dagher J, Pladys A, Edeline J, et al. Independent Association of PD-L1 Expression With Noninactivated VHL Clear Cell Renal Cell Carcinoma-A Finding With Therapeutic Potential. Int J Cancer (2017) 140:142–8. doi: 10.1002/ijc.30429
42. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant’ Angelo M, Forman J, et al. Interplay of Somatic Alterations and Immune Infiltration Modulates Response to PD-1 Blockade in Advanced Clear Cell Renal Cell Carcinoma. Nat Med (2020) 26:909–18. doi: 10.1038/s41591-020-0839-y
43. Wei SC, Duffy CR, Allison JP. Fundamental Mechanisms of Immune Checkpoint Blockade Therapy. Cancer Discov (2018) 8:1069–86. doi: 10.1158/2159-8290.Cd-18-0367
44. Killock D. Immunotherapy: Interferon in Anti-CTLA-4 Responses. Nature Reviews. Clin Oncol (2016) 13:652–3. doi: 10.1038/nrclinonc.2016.166
45. Zheng F, Dang J, Zhang H, Xu F, Ba D, Zhang B, et al. Cancer Stem Cell Vaccination With PD-L1 and CTLA-4 Blockades Enhances the Eradication of Melanoma Stem Cells in a Mouse Tumor Model. J Immunother (2018) 41:361–8. doi: 10.1097/cji.0000000000000242
46. Díaz-Montero CM, Rini BI, Finke JH. The Immunology of Renal Cell Carcinoma. Nat Rev Nephrol (2020) 16:721–35. doi: 10.1038/s41581-020-0316-3
47. Plaks V, Kong N, Werb Z. The Cancer Stem Cell Niche: How Essential is the Niche in Regulating Stemness of Tumor Cells? Cell Stem Cell (2015) 16:225–38. doi: 10.1016/j.stem.2015.02.015
48. Motzer RJ, Penkov K, Haanen J, Rini B, Albiges L, Campbell MT, et al. Avelumab Plus Axitinib Versus Sunitinib for Advanced Renal-Cell Carcinoma. N Engl J Med (2019) 380:1103–15. doi: 10.1056/NEJMoa1816047
49. Rini BI, Plimack ER, Stus V, Gafanov R, Hawkins R, Nosov D, et al. Pembrolizumab Plus Axitinib Versus Sunitinib for Advanced Renal-Cell Carcinoma. N Engl J Med (2019) 380:1116–27. doi: 10.1056/NEJMoa1816714
50. Gaujoux R, Seoighe C. A Flexible R Package for Nonnegative Matrix Factorization. BMC Bioinf (2010) 11:367. doi: 10.1186/1471-2105-11-367
51. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression. Genome Biol (2016) 17:218. doi: 10.1186/s13059-016-1070-5
52. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, et al. Genomics of Drug Sensitivity in Cancer (GDSC): A Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res (2013) 41:D955–61. doi: 10.1093/nar/gks1111
Keywords: clear cell renal cell cancer (ccRCC), stem cell, subtype, gene expression, immune microenvironment
Citation: Wang H, Xu H, Cheng Q and Liang C (2021) Identification of a Novel Stem Cell Subtype for Clear Cell Renal Cell Carcinoma Based on Stem Cell Gene Profiling. Front. Oncol. 11:758989. doi: 10.3389/fonc.2021.758989
Received: 15 August 2021; Accepted: 02 November 2021;
Published: 29 November 2021.
Edited by:
Andrea Mari, Careggi University Hospital, ItalyReviewed by:
Gian Maria Busetto, University of Foggia, ItalyMohammed Imran Khan, Western University, Canada
Copyright © 2021 Wang, Xu, Cheng and Liang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chaozhao Liang, bGlhbmdfY2hhb3poYW9AYWhtdS5lZHUuY24=