Integrated the analysis of single-cell and bulk RNA-sequencing to predict renal cell carcinoma by identifying signature based on Mitophagy-related Genes

doi:10.21203/rs.3.rs-4930192/v1

Download PDF

Research Article

Integrated the analysis of single-cell and bulk RNA-sequencing to predict renal cell carcinoma by identifying signature based on Mitophagy-related Genes

https://doi.org/10.21203/rs.3.rs-4930192/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

The features of resistance to mitophagy contribute significantly to invasion, malignancy and cell survival. But the mechanism of mitophagy in clear cell renal cell carcinoma (ccRCC) remains unclear It is valuable to estimate mitophagy molecular characters as a clinical factor for prognosis and immune phenotypes in ccRCC.

Methods

Clinical data of ccRCC patients, including genome and transcriptome data, were downloaded from The Cancer Gene Atlas (TCGA) and International Cancer Genome Consortium (ICGC) database. The differentially expressed genes (DEGs) of patient clusters determined by mitophagy gene expression and univariate Cox regression analysis were identified and used to classify patient clusters for constructing mitophagy scores via PCA analysis. Immune cell infiltration and immune cell function were analyzed by ssGSEA algorithm and TIDE algorithm.

Results

Based on the expression of mitophagy marker genes, ccRCC patients were divided into three mitophagy clusters with different gene expression patterns, prognosis and immune niches. 1,356 DEGs of mitophagy clusters related with prognosis were screened out for building mitophagy score. ccRCC patients with high mitophagy scores have better prognosis. Meanwhile, lower mitophagy patients with high expressed several immune-checkpoint proteins and had high immunophenoscore after immune-checkpoint blockers treatment, indicating better responsiveness to immune therapy.

Conclusions

mitophagy features are tightly correlated with ccRCC prognosis and immune responsiveness. mitophagy score built here is able to predict the prognosis and immune features of ccRCC patients and be indicative for immunotherapy.

renal cell carcinoma

single cell

mitophagy

TCGA

ICGC

Renal carcinoma incidence is on the rise by year, especially tend to be younger population[1, 2]. Recent years, there were more than 431,000 new cases of renal carcinoma and 179,000 deaths in globally[3]. ccRCC is the most common histologic subtype of kidney cancer, with extensive tumor heterogeneity[4]. ccRCC is frequently genetically varied, and the somatic mutations such as VHL, PBRM1, SETD2, BAP1, KDM5C and so on[5]. Although nephrectomy has shown well efficacy in treating localized ccRCC, more than 30% of patients go through advanced disease progression and 25% cancer recurrence[6].

Mitochondrial autophagy is a subtype of autophagy that occurs in mitochondria, which specifically removes damaged and aging mitochondria, thus guaranteeing the quantity and quality of mitochondrial populations in the cells[7, 8]. When mitochondria are defective, some special receptors recognize them, subsequently creating a separator/phagosome around the defective mitochondria for selective degradation as autophagosomes[9]. Many studies suggested that it is commonly associated with many diseases due to the protective mechanism of mitophagy[10–12]. However, mitochondrial autophagy exists "double-edged sword" effect in many cancers, which can both promote and inhibit tumor progression. For one thing, mitophagy clears dysfunctional mitochondria and protects cells from intracellular oxidation and genotoxic stress, by preventing tumorigenesis in early-stage tumors. For another, when tumor formed, tumor cells undergo mitosis to reduce oxidative stress and provide a matrix for circulating cell proliferation and survival, which is more common in higher grade tumors or malignancies. Mitochondrial autophagy has been shown to be associated with metabolic reprogramming, anti-cancer therapy resistance, and stem cell sexuality of cancer cells[13, 14]. Therefore, based on the role of mitophagy in cancer, screening mitophagic biomarkers will help target tumors, which may be an alternative and precise method with higher anti-cancer efficiency and fewer side effects.

In this study, we performed comprehensive bioinformatic analysis of molecular characters mitophagy in ccRCC patients and revealed the relationship between mitophagy features and tumor prognosis as well as immune phenotypes. By constructing a scoring system based on the mitophagy related gene expression pattern, we suggested molecular characters of mitophagy are reliable prognostic signatures and provided a novel tool to evaluate ccRCC and immune properties.

2.1 Data collection and processing

In this study, single-cell RNA sequencing data of two ccRCC samples were collected from GSE152938 by Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. Bulk transcriptome data from a total of 621 ccRCC samples was collected from The Cancer Gene Atlas (TCGA, 530 samples) from University of California Santa Cruz (UCSC, https://xenabrowser.net/datapages/) Xena browser and International Cancer Genome Consortium (ICGC, 91 samples, https://dcc.icgc.mitophagy/). The TPM values of TCGA and ICGC were transformed from the FPKM values. The datasets of TCGA and ICGC were combined as a merge matrix and “ComBat” algorithm of R package was used to eliminate the batch effect from non-biological technical biases of each dataset[15]. Clinical information of samples including age, gender, survival status, tumor type, tumor grade, which obtained from TCGA and ICGC clinical information 2.2 Single cell RNA sequence data processing.

scRNA-seq data is filtered and calibrated using the R 4.4.0 "Seurat" and "SingleR" packages. We screened cells with unique feature counts > 4500 or < 200 and mitochondria cell counts > 10%. Total expression normalization of the feature expression measurements for each cell was performed using the default parameters of Seurat's NormalizeData function. Using “Harmony” package to transfer all cell data to the composite seurat object. Use the runtsne function to scale and calculate the main component of the gene (min.dist = 0.2, neighbors = 20). Selecting the "FindClusters" function (resolution = 0.2) and significant principal components for TSNE analysis and cluster analysis.

2.3 Clusters of cell annotation

To identify the cell type, we performed a total of two annotation modes. Automatic annotation: scMaymoMap is a way to automatically annotate scRNAseq data [16]. Given a sample reference dataset with a known label (single cell or batch size), it labels new cells in the dataset based on their similarity to the reference. Thus, for reference datasets, the burden of manually interpreting clustering and defining tagged genes only needs to be done once, and this biological knowledge can be propagated to new datasets in an automated manner. According to “limma” package, DEGs of single cell were screened and criteria is |logFC|>0.25, P < 0.05.

2.4 Cell interaction analysis

CellChat is a database containing information on ligands, receptors and their interactions. The database can be used for comparative reasoning analysis and quantitative description between communication network units [17]. Cell-cell communication analysis using the R "Cellchat" package. The reference human ligand receptor database is CellChatDB. Human cell- cell communication is determined by assessing ligand and receptor expression in CellChatDB (R-pack CellChat 0.0.2). We studied the interactions between different cell types, screening for pathways where the number of cells is less than 10. We next used CellphoneDB package[18] to further analysis the interaction between cancer cell and T cell.

2.5 Mitophagy pathway related gene collected

22 mitophagy pathway related genes were downloaded from Gene Set Enrichment Analysis (GSEA, https://www.gsea-msigdb.org/). According to Venn diagram, the overlapped mitophagy genes were flittered from DEGs of single cell and mitophagy pathway.

2.6 Unsupervised clustering analysis of ccRCC based on overlapped genes expression

The expression of mitophagy marker genes was extract from the merge matrix. According to the expression, we executed hierarchical agglomerative clustering of each sample. R package “ConsensusClusterPlus” was used to determine the number of clusters and stability, which is based on unsupervised clustering “Pam” method according to Euclidean and Ward’s linkage, and 500 times repeats were conducted to make sure the clustering stability[19].

2.7 Identification of DEGs associated with the mitophagy-related Phenotype

All patients were classified into different mitophagy-clusters based on mitophagy marker genes expression to identify DEGs among distinct mitophagy clusters. R package “limma” was used to identify the DEGs among mitophagy clusters and the cutoff criteria was determined as adjusted p < 0.05 and |Log₂ foldchange| > 1[20].

2.8 Mitophagy score construction

Univariate Cox regression model was performed for the prognostic analysis for each DEGs. The genes with significant prognosis relevance were extracted as selected DEGs for further analysis. We then conducted principal component analysis (PCA) based on selected DEGs for ccRCC patients. Both principal component 1 and 2 were selected to depict the mitophagy signature for each patient and both principal component 1 and principal component 2 of each patient from PCA were extracted for the mitophagy score construction. Finally, we defined the mitophagy score for each patient as follow[21]:

$$\:Mitophagy\:score=\:\sum\:PCA1+\:\sum\:PCA2$$

2.9 Somatic gene mutation analysis

The corresponding somatic mutation data of TCGA was obtained from UCSC Xena browser (https://xenabrowser.net/datapages/). Tumor mutation information of each patient were collected. To determine the relationship between somatic gene mutation and mitophagy score, we grouped patients into low and high mitophagy score subgroups by R package “Maftool”[22]. The top 20 driver mutation genes were further analyzed.

2.10 Immune Cell Infiltration Analysis and relation with mitophagy score

R algorithm “ssGSEA” was employed via R packages (“GSVA”, “GSEABase”, and “limma”) to comprehensively assess the immunologic characteristics of each patient[23]. Pearson correlation analysis was further exploited to calculate the correlation between immune cells and mitophagy score. The expression of immune checkpoint genes was extracted from ccRCC gene expression matrix and R package “limma” was used to analyze correlation of expression of immune checkpoint genes with mitophagy scores.

2.11 Genomic data and clinical information with Immunotherapy

To further validate the prognostic significance of the mitophagy score, we applied the Immunophenotype scores (IPS) to predict response to immunotherapy by immune checkpoint inhibitors. IPS of patients were collected from the Cancer Immunome Atlas (TCIA, https://tcia.at/home)[24].

2.12 Cell culture

KIRC cell line SW839 and HEK-293T were all obtained from American Type Culture Collection (ATCC, Rockville, MD, USA). Cells were cultured in RPMI-1640 containing 10% fetal bovine serum (FBS, Gibco, Carlsbad, CA, USA), 1% penicillin-streptomycin (HyClone, Logan, UT, USA), and then preserved in the atmosphere with 5% CO₂ at 37°C

2.13 Real-time quantitative PCR

RNA was extracted from 293T and SW893 using TRIzol (Thermo Fisher Scientific). The purity and concentration of the extracted RNA were determined. In addition, 500 mg RNA was reversetranscribed into complementary DNA (cDNA) using PrimeScript TR premix reagent (RR036A, TaKaRa). The reaction mixture was then prepared for quantitative PCR (qPCR) using cDNA, primer (Sangon Biotech) and PowerUp SYBR Green (A25742, Thermo).

2.14 Cell proliferation assay

In total, 2,000 cells/well were seeded in 96-well plates. After culturing the cells for a period of time, the cells were incubated solution for 1 h. The growth of the cells was detected by measuring the OD value at 480-nm wavelength.

2.15 5-Ethynyl-20-deoxyuridine (EdU) assay

EdU incorporation assay was performed using Clict-EdU according to the manufacturer’s instructions. The cells were fixed in 4% paraformaldehyde, incubated with 20 uM EdU for 2 h, and again incubated in PBS + 3% BSA. Then, the click-iT reaction reagents were added to the cells followed by 0.5-h incubation. The samples were blocked using PBST blocking buffer for 1 h at room temperature. The cells were incubated overnight with an anti-biotin primary antibody diluted in a blocking solution.

2.16 Survival Analysis

The survival analysis in patients was further assessed by the Kaplan-Meier plotter (KM plotter, http://www.kmplot.com/analysis/). About all samples were collect, and their KM survival curves were conducted.

2.17 Statistical analysis

R version 4.2.1 is used for all statistical analyses. The two groups were compared using the Wilson test. The Kruskal-Wallis test was used for more than 3 groups. Kaplan-Meier (K-M) plots were used to show survival differences between patients and log-rank tests. The correlation coefficients were analyzed by Pearson correlation analysis. A P-value of less than 0.05 indicates statistical significance. χ2 test was used to analyze the cluster of high and low mitochondrial autophagy score and somatic mutation frequency, and the correlation was analyzed. The correlation coefficients were calculated by Pearson correlation analysis. A two-tailed P < 0.05 was considered statistically significant.

3.1 Clustering of ccRCC single cells

After data processing and filtrating, ccRCC single-cell data were grouped into 12 clusters, and the cell types were annotated using classical signature genes as markers (Fig. 1A), which were annotated as a total of 8 cell types, including dendritic cell, endothelial cell, cancer cell, monocyte, natural killer cell, pericyte Proximal tubule brush border cell (PTBBC) and T cells (Fig. 1B-1C). The heatmap of annotated cell performed in Fig. 1D and the bubble plots was visualized the expression of classical marker genes which shown in Fig. 1E. The expression of ccRCC specific markers (HILPDA, ANGPTL4, NNMT, SLC17A3, NDUFA4L2 and CA9) were high expression in cancer cell that further verified the annotation of cancer cell (Fig. 1F-1K). The correlation among all cell types were show in Fig. 2A, and we found cancer cell and T cell were significantly associated, and the interaction may be closely related to the immune microenvironment of the tumor. After that, DEGs among all cell types performed in Fig. 2B. Total of 2,271 DEGs were screened in cancer cell. The most upregulated gene including NNMT and MT1E and the most downregulated genes including HLA-DRA, S100A9 and APOE (Fig. 2B).

3.2 Communication analysis of single cell types

We performed cell communication analysis to further investigate the interactions between cells. We found that there is a high intensity interaction between each cell. We found strong communication between cancer cell and T cells, and receptor-ligand interactions analysis shown TNFSF9-TNFrSF9, NECTIN2-TIGIT, NECTIN2-NECTIN3, SIRPA-CD47.1, SIRPA-CD47, CD70-CD27, CD58-CD2, CXCL14-CXCR4 (Fig. 2C). These data suggest that there are many immunobiological processes in cancer cells and T cells. Through signal pathway analysis, it was found that PD-L1 signaling pathway mainly originates from endothelial cells and T cells, and has significant correlation with T cells, dendritic cells and cancer cells (Fig. 2D). The PD-L2 signaling pathway is mainly derived from dendritic cells and monocytes, and is significantly correlated with PTBBC and NK cells (Fig. 2E). The CEACAM signaling pathway showed that cancer cells were positively associated with other cells, especially dendritic cells and monocytes (Fig. 2F). The CD39 pathway was mainly derived from cancer cells, NK, PTBBC and T cell and was significantly correlated with each cell types (Fig. 2G). The CD48 pathway was mainly derived from cancer cells, epithelial cell and PTBBC and was significantly correlated with each cell types (Fig. 2H). The CD96 signaling pathway was mainly derived from dendritic cell, endothelial cell, cancer cells, monocyte, pericyte and PTBBC and was significantly correlated with each cell (Fig. 2I). The CD45 signaling pathway is mainly derived from endothelial cell, cancer cells and PTBBC and has a significant correlation with each cell types expect dendritic cell (Fig. 2J). Based on further study of cell-cell interactions. We found that there is a high intensity interaction between each cell. Interestingly, we found strong communication between cancer cells, endothelial cells, dendritic cells, T cells, PTBBC, NK and monocytes (Fig. 3).

3.3 Expression and landscape of mitophagy marker genes in ccRCC

According to GSEA database, we collected 39 related mitophagy gene and screened 17 overlapped gene from DGEs of scRNA (Fig. 4A). In order to unveil the mitophagy characters in ccRCC patients, the 17 overlapped genes were used to represent mitophagy feature and detect their expression level to reflect mitophagy properties in ccRCC. These marker genes comprised IRGM, VDAC1, OPTN, SLC25A5, LRBA, RETREG1, PHB2, HDAC6, HTRA2, SLC25A4, VPS13C, CERS1, HUWE1, PINK1 MAP1LC3B, BECN1 and CDC37. We firstly checked bulk expression of 17 overlapped genes between ccRCC and normal tissue. The results shown PINK1, PHB2, LRBA, SLC25A4, SLC25A5, HUWE1, BECN1, VPS13C and HDAC6 were high expression in normal tissues than ccRCC. CDC37, CERS1, VDAC1, OPTN, IRGM and RETREG1 upregulated in ccRCC than normal tissue (Fig. 4B). Based on scRNA, PINK1, PHB2, MAP1LC3B, BECN1, SLC25A4, CDC37, SLC25A5, HUWE1, RETREG1, VDAC1, OPTN, VPS13C exist high ontology expression and LRBA, HDAC6, CERS1, HTRA2 have low ontology expression in RCC (Figure S1). All of the marker genes have copy number variation (CNV) in ccRCC. For IRGM and VDAC1, the frequency of copy amplification is obviously higher than copy deletion, which is in contrast to PINK1 and SLC25A4 (Fig. 4C-4D). The genomic alterations of mitophagy markers in ccRCC were performed in Fig. 4E. Somatic mutation analysis revealed mutation frequency and mutation type of each mitophagy marker genes. HUWE1 and LRBA have the highest mutation frequency. Then the overall survival of mitophagy marker genes which exist prognostic to ccRCC were shown in Figure S2. All of these data suggested mitophagy characters are highly possible to corelate with ccRCC progress and be used to depict ccRCC malignancy.

3.4 Classification of ccRCC patients based on mitophagy characters

The expression and mutation pattern of mitophagy marker genes varied among ccRCC patients and were related to tumor prognosis. Therefore, we tried to make subtype classification for ccRCC according expression of 17 mitophagy marker genes to further study impacts of mitophagy features on tumor pathology. The ccRCC patients were divided into 3 clusters (we name them mitophagy cluster A, B, C respectively, Fig. 4F-4H). Each mitophagy cluster has a specific expression pattern of mitophagy markers (Fig. 4I). Principle component analysis (PCA) showed that patients within the same mitophagy cluster share similar transcriptional profiles, indicating mitophagy characters are compatible to be applied to distinguish ccRCC patients with different pathological features (Fig. 4J). And We further verified the overall survival of 3 cluster, and find patients in cluster exist worse survival rate (Fig. 4K). Next, we made pairwise comparisons of gene expression profiles among mitophagy clusters for identification of differential expression genes (DEGs). After that Gene Set Variation Analysis (GSVA) analysis pathways was conducted on each set of mitophagy cluster to pinpoint out the most discrepant molecular networks (Fig. 5A-5C). We also analyzed the immune cell infiltration (ICI) in ccRCC samples of three mitophagy clusters and the correlation of each ICI type were performed in Fig. 5D. We Next found that the infiltration of all the ICI types varies significantly among mitophagy clusters, indicating each mitophagy cluster has a specific immune microenvironment (Fig. 5E).

3.5 Reconstruction of tumor clusters in ccRCC

To further explore mitophagy characters in ccRCC progression and construction new model for ccRCC prognosis prediction, we focused on the intersection of DEGs between mitophagy cluster A and B, mitophagy cluster B and C, mitophagy cluster A and C. This gene set contains 2,413 genes which are enriched in diverse gene ontological terms (Fig. 6A). After correlation analysis with ccRCC prognosis, we screened out 1949 genes, and we named this gene set “mitophagy related gene set”. The GO analysis of 1,949 gene shown in Fig. 6B. Subsequently, we reclassified ccRCC patients according to expression of mitophagy related gene set into three gene clusters, called gene cluster A, B and C respectively (Fig. 6B-6D). Like classification by mitophagy marker genes, patients from the same gene cluster have a similar gene expression pattern of mitophagy related gene set in PCA analysis. What’s more, Kaplan-Meier analysis revealed considerable difference of survival probabilities in patients of three gene clusters, represented as patients in gene cluster A with best survival probabilities and patients in gene cluster C with worst survival probabilities (Fig. 6E). The ICI analysis performed in Fig. 6F.

3.6 Construction of mitophagy scoring model for quantification of mitophagy characters and prognosis evaluation in ccRCC

We calculated mitophagy scores for individual ccRCC patient and classify patients into high mitophagy score group and low mitophagy score group referring to the median of all mitophagy scores and Sankey plot elucidated the attribution change of patients across mitophagy clusters, gene clusters, mitophagy scores groups and survival states (Fig. 7A). To quantify mitophagy characters in each ccRCC patient, we set up mitophagy scoring system based on PCA with mitophagy related gene set. Then we tested the mitophagy scoring in mitophagy clusters and gene clusters and found there are significant differences in mitophagy scores of clusters divided by both methods, indicating the level of mitophagy scores can indicate mitophagy characters in ccRCC clearly and comprehensively (Fig. 7B-7C). The prognostic of mitophagy show patient with high score exist poor survival rate (Fig. 7D). Tumor mutation burden (TMB) is a prognostic factor that tumor with high TMB is associated with poor survival[25]. We analyzed TMB overall survival and find high TMB value exist poor prognosis (Fig. 7E). Furthermore, the we found high score with high TMB value tend to worse survival (Fig. 7F), which coincide formal survival analysis. Gene Set Enrichment Analysis (GSEA) revealed significant enrichment of signal pathways related to hallmarks of cancer in patients with low mitophagy scores over patients with high mitophagy scores (Fig. 7G). Wilcoxon test was performed, and the results revealed a significant increase in the expression of CTLA4, IFNG, GZMB, PDCD1, CD8A, TBX2, LAG3, GZMA and PRF1 in the low-ICI-score group compared with patients in the high-ICI-score group, and HAVCR2, and CD274 were high expressed in patients with high mitophagy score (Fig. 7H).

As to investigate the genomic alteration profiles corresponding to mitophagy scores, we examine the mutation status of the top 20 most frequently mutated genes in the high mitophagy score group and the low mitophagy score group (Fig. 7I-7J). Despite the overall higher mutation rate in low mitophagy score group, mutation rate of some genes which are identified as most frequently mutated genes in cancer mutation analysis in the high mitophagy score group is different with the low mitophagy score groups[26], such as PBRM1, TTN and SETD2.

3.7 The relationship of clinicopathologic features and mitophagy score

As above data, we confirmed that mitophagy score act as independent factor to predict ccRCC prognosis. For furthering explore the guide of clinical significance, we analyzed the clinicopathologic features of mitophagy score. TNM staging is the most significant guide to judgment the grade of renal carcinoma, the higher level of TNM staging predict higher malignancy of tumor[27]. According to this analysis, we find patients with elevated mitophagy tend to link with higher TNM stage (Fig. 8A-8C). Meanwhile, we also discover that high mitophagy score correlated with high stage and grade of ccRCC (Fig. 8D-8E). Furthermore, high mitophagy patients tend to high death accidents (Fig. 8F).

Immune checkpoint inhibition has been at the center of cancer immune therapy and is involved extensively in clinical trials for ccRCC treatment[28]. Thereafter, we explored the potential of mitophagy characters as a predictive factor for therapies incorporating immune checkpoint inhibitors. We inspected the expression of typical immune checkpoint related genes and immune-activation genes in the high mitophagy score group and the low mitophagy score group and all of these genes have higher expression level in the high mitophagy score group, implying patients in high mitophagy score group could respond to immune checkpoint blocking with better performance (Fig. 8G-8I). Actually, predictions by immunophenoscore (IPS) confirmed that patients in high mitophagy score groups have better response to anti-CTLA4 antibody monotherapy than combinative treatment embracing both anti-PD-1 antibody and anti-CTLA-4 antibody. Thus, mitophagy characters are magnificently informative for immune phenotypes of ccRCC, which is suggestive for management of glioma immune therapies.

3.8 Significance of mitophagy score for patient prognosis

In view of the significant significance of mitochondrial autophagy score in predicting the prognosis of ccRCC patients, this paper aims to establish a clinical prediction model to predict the survival probability of ccRCC patients based on clinical characteristics. Firstly, univariate and multivariate COX regression models were used to screen out 6 clinical features including mitochondrial autophagy score, age, grade, class T, class M and stage classification (Fig. 9A-9B). Then we constructed a nomograph based on four clinical features which were easily accessible (Fig. 9C). The Cindexes indicated that the nomogram had a good predictive value (Figs. 9D-9F).

3.9 HSPA8 act as antioncogenic effect in RCC

HSPA8 a hubgene of PPI network with significant differences between high mitophagy score group and low mitophagy score group (Fig. 10A), to validate the rationality of the scoring system. The expression level of HSPA8 in tumor samples was significantly lower than that in normal by TCGA (Fig. 10B) and patients with high level of HSPA8 tend to have better survival (Fig. 10C). According to scRNA of ccRCC, we found HSAP8 high expression in ccRCC cell (Fig. 10D-10E). And according to IHC and qPCR assay, HSPA8 were low expression in tumor (Fig. 10F-10G). Furthermore, knockdown HSPA8 decreased the proliferation of RCC (Fig. 10H) and reduce RCC cell cycle (Fig. 10I). Final, after knockdowning HSPA8, the RCC cell growth reduced (Fig. 10J)

Mitochondrial autophagy is a key mechanism of metabolic reprogramming and aerobic glycolysis regulation in tumor cells, and it has broad prospects as a new method for the treatment of various cancers [29]. Mitochondrial autophagy plays a dual role in PC progression, depending on its clinical features, tumor microenvironment, and mutation status [30]. Even though the role of mitochondrial autophagy in ccRCC tumorigenesis has been well documented, the extent to which mitochondrial autophagy influences ccRCC prognosis and treatment strategies is unclear. In this study, we comprehensively analyzed the value of mitochondrial autophagy associated genes in risk stratification, immune activity, gene mutation, and chemotherapy response in ccRCC patients. Firstly, based on scRNA, we screened DGEs of RCC cancer cell type and overlapped with the genes of mitophagy pathway. We found that overlapped genes could serve as a ccRCC classifier. Our category-based mitochondrial autophagy gene profile was developed to stratify patients with ccRCC. The high-risk group had low survival rate, high mutation load and high sensitivity to treatment. Our prognostic model consists of genes associated with mitochondrial autophagy (IRGM, VDAC1, OPTN, SLC25A5, LRBA, RETREG1, PHB2, HDAC6, HTRA2, SLC25A4, VPS13C, CERS1, HUWE1, PINK1 MAP1LC3B, BECN1 and CDC37). We also find overlapped mitophagy gene could act as intendent factors to predict ccRCC prognosis. Accordingly, these mitophagy-related signature genes play a vital role during ccRCC tumorigenesis and might have clinical value as potential biomarkers.

Development of mitophagy scoring model and its application for ccRCC patients’ classification in our work could be a helpful complement to traditional classification for ccRCC diagnosis and treatment, providing a new paradigm beyond genetic characterization of ccRCC, which could lead to benefits to clinical practice. Our work uncovered the relations between ccRCC pathology and mitophagy characters by identifying genomic mutation features corresponding to the high mitophagy score group or the low mitophagy score group. Such effort deepened the understanding of the molecular pathology of ccRCC and shed light on the identification of new therapeutic targets with specificity against ccRCC patients with different cellular hierarchy and genetic background. Our analysis delineated the clear correlation of mitophagy characters with ccRCC immune phenotypes, highlighting the importance of mitophagy model as a subtype in ccRCC immunology and bringing insights on explorations of contributions of mitophagy heterogeneity to immunosuppression of ccRCC. mitophagy related pathways have been reported to participate in ccRCC cell survival and suggested to be drug targets for ccRCC therapy[31]. The knowledge about molecular and cellular mechanisms of mitophagy pathways in ccRCC pathology is limited. Our work for the first the time couple’s endosome pathways with a certain subtype of ccRCC mitophagy subtype and highlights its functions in ccRCC progression, which is enlightening for researches about subtype specific functions of signal pathways in tumor progenitor cells and cancer biology.

Mutation identification is also a key step in cancer risk assessment, which may help develop prevention or treatment strategies for ccRCC. In our study, we further explored the mutational landscape of mitochondrial autophagy associated genetic signatures. Somatic mutation count exhibited different between high mitophagy group and low mitophagy. Among them, the incidence of BPRM1 mutation in the low mitophagy group was higher than the high mitophagy group among them (29% vs 38%). Research report that BPRM1 mutations may associated with a more favorable course[32], which may contribute poor survival rate of patients with high mitophagy score. In addition, mutations in SETD2 are associated with increased DNA methylation loss in non-promoter regions. This may be due to reduced trimethylation (H3K36me3) of histone H3 (H3K36), which is involved in maintaining the heterochromatic state[33]. Former research suggested mutations to SETD2 occur in 10% of ccRCC. SETD2-mutated ccRCCs have been related to worse overall survival[34]. According to our result we explore that SETD2 high mutated in samples with high mitophagy score, which may result in poor prognosis of this group.

HSPA8 (member of the Heat Shock Protein Family A (Hsp70) 8), known as a homologous protein of the HSPA/HSP70 family, is a conserved and constitutively expressed chaperone in the cell membrane, cytosol, and nucleus. It regulates cellular protein homeostasis by mediating protein assembly, renaturation, transport and degradation[35, 36]. HSPA8 play a vital function of regulating autophagy in cancers[37, 38]. Whereas, the regulation role of HSPA8 to mitophagy in unclear. In this research, we find HSPA8 is the hub key gene within high and low mito-score group. And HSPA8 contribute RCC cell proliferation and cell cycle.

In conclusion, we developed and validated a mitochondrial autophagy associated gene model for prognostic stratification in patients with ccRCC. In addition, this new profile can be used to predict survival and treatment response during treatment in patients with ccRCC.

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.

Data availability statement

In this study, data were generated at the TCGA, GEO and ICGC database with the cancer type of prostate cancer. Derived data supporting the findings are available from the corresponding author on reasonable request.

Acknowledgments

In this study, results shown here are based on data publicly available at the TCGA database. We thank all patients and all investigators who participated in TCGA dataset for supplying data.

Authors’ Contributions

All authors contributed to the conception and design of the study. Guanbao Tang contributed to the acquisition of samples and data analysis. Jiwei Yang and Dongdong Chen analyzed data and wrote the manuscript. Xuewen Guo, Hao Chen, lizhun Wang and Tongyi Men critically revised the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

In this study, all human data was acquired from public database of TCGA, GEO, ICGC and HPA. This study does not involve human experiments, nor does it involve clinical studies. Human Ethics and Consent to Participate declarations: not applicable.

Funding

No Funding in this study

Siegel, R.L., et al., Cancer statistics, 2022. CA Cancer J Clin, 2022. 72(1): p. 7-33.
Miller, K.D., et al., Cancer statistics for adolescents and young adults, 2020. CA Cancer J Clin, 2020. 70(6): p. 443-459.
Sung, H., et al., Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin, 2021. 71(3): p. 209-249.
Argentiero, A., et al., Anti-angiogenesis and Immunotherapy: Novel Paradigms to Envision Tailored Approaches in Renal Cell-Carcinoma. J Clin Med, 2020. 9(5).
Lin, P.H., et al., Genomic characterization of clear cell renal cell carcinoma using targeted gene sequencing. Oncol Lett, 2021. 21(2): p. 169.
Kotecha, R.R., R.J. Motzer and M.H. Voss, Towards individualized therapy for metastatic renal cell carcinoma. Nat Rev Clin Oncol, 2019. 16(10): p. 621-633.
Onishi, M., et al., Molecular mechanisms and physiological functions of mitophagy. EMBO J, 2021. 40(3): p. e104705.
Poole, L.P. and K.F. Macleod, Mitophagy in tumorigenesis and metastasis. Cell Mol Life Sci, 2021. 78(8): p. 3817-3851.
Panigrahi, D.P., et al., The emerging, multifaceted role of mitophagy in cancer and cancer therapeutics. Semin Cancer Biol, 2020. 66: p. 45-58.
Ma, X., et al., Role and Mechanisms of Mitophagy in Liver Diseases. Cells, 2020. 9(4).
Thangaraj, A., et al., Mitigation of cocaine-mediated mitochondrial damage, defective mitophagy and microglial activation by superoxide dismutase mimetics. Autophagy, 2020. 16(2): p. 289-312.
Han, S., et al., Mitophagy coordination with retrograde transport ensures the integrity of synaptic mitochondria. Autophagy, 2020. 16(10): p. 1925-1927.
Vara-Perez, M., B. Felipe-Abrio and P. Agostinis, Mitophagy in Cancer: A Tale of Adaptation. Cells, 2019. 8(5).
Kulikov, A.V., et al., Mitophagy: Link to cancer development and therapy. Biochem Biophys Res Commun, 2017. 482(3): p. 432-439.
Johnson, W.E., C. Li and A. Rabinovic, Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics, 2007. 8(1): p. 118-27.
Yang, L., et al., Single-cell Mayo Map (scMayoMap): an easy-to-use tool for cell type annotation in single-cell RNA-sequencing data analysis. BMC Biol, 2023. 21(1): p. 223.
Jin, S., et al., Inference and analysis of cell-cell communication using CellChat. Nat Commun, 2021. 12(1): p. 1088.
Efremova, M., et al., CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc, 2020. 15(4): p. 1484-1506.
Wilkerson, M.D. and D.N. Hayes, ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics, 2010. 26(12): p. 1572-3.
Ritchie, M.E., et al., limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res, 2015. 43(7): p. e47.
Sotiriou, C., et al., Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst, 2006. 98(4): p. 262-72.
Mayakonda, A., et al., Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res, 2018. 28(11): p. 1747-1756.
Hänzelmann, S., R. Castelo and J. Guinney, GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics, 2013. 14: p. 7.
Jiang, P., et al., Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med, 2018. 24(10): p. 1550-1558.
Valero, C., et al., The association between tumor mutational burden and prognosis is dependent on treatment context. Nat Genet, 2021. 53(1): p. 11-15.
Mendiratta, G., et al., Cancer gene mutation frequencies for the U.S. population. Nat Commun, 2021. 12(1): p. 5961.
Rini, B.I., S.C. Campbell and B. Escudier, Renal cell carcinoma. Lancet, 2009. 373(9669): p. 1119-32.
Miao, D., et al., Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma. Science, 2018. 359(6377): p. 801-806.
Liu, T., et al., GPD1L inhibits renal cell carcinoma progression by regulating PINK1/Parkin-mediated mitophagy. J Cell Mol Med, 2023. 27(16): p. 2328-2339.
Xiang, J., et al., Deciphering the implications of mitophagy-related signatures in clinical outcomes and microenvironment heterogeneity of clear cell renal cell carcinoma. J Cancer Res Clin Oncol, 2023.
Yan, M., et al., Knockdown of NR3C1 inhibits the proliferation and migration of clear cell renal cell carcinoma through activating endoplasmic reticulum stress-mitophagy. J Transl Med, 2023. 21(1): p. 701.
Carlo, M.I., et al., Genomic Alterations and Outcomes with VEGF-Targeted Therapy in Patients with Clear Cell Renal Cell Carcinoma. Kidney Cancer, 2017. 1(1): p. 49-56.
Comprehensive molecular characterization of clear cell renal cell carcinoma. Nature, 2013. 499(7456): p. 43-9.
Piva, F., et al., BAP1, PBRM1 and SETD2 in clear-cell renal cell carcinoma: molecular diagnostics and possible targets for personalized therapies. Expert Rev Mol Diagn, 2015. 15(9): p. 1201-10.
Radons, J., The human HSP70 family of chaperones: where do we stand? Cell Stress Chaperones, 2016. 21(3): p. 379-404.
Stricher, F., et al., HSPA8/HSC70 chaperone protein: structure, function, and chemical targeting. Autophagy, 2013. 9(12): p. 1937-54.
Losmanova, T., et al., Chaperone-Mediated Autophagy Markers LAMP2A and HSPA8 in Advanced Non-Small Cell Lung Cancer after Neoadjuvant Therapy. Cells, 2021. 10(10).
Tian, Y., et al., Maslinic acid induces autophagy by down-regulating HSPA8 in pancreatic cancer cells. Phytother Res, 2018. 32(7): p. 1320-1331.

No competing interests reported.

FigureS1.tif
Figure S1 (A-P) All overlapped mitophagy-related genes expression in RCC scRNA
FigureS2.tif
Figure S2 (A) Kaplan-Meier curve of mitophagy signature gene with prognostic significant.

Download PDF

Version 1

posted

You are reading this latest preprint version

Integrated the analysis of single-cell and bulk RNA-sequencing to predict renal cell carcinoma by identifying signature based on Mitophagy-related Genes

Status:

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

1 Introduction

2 Material and methods

2.1 Data collection and processing

2.3 Clusters of cell annotation

2.4 Cell interaction analysis

2.5 Mitophagy pathway related gene collected

2.6 Unsupervised clustering analysis of ccRCC based on overlapped genes expression

2.7 Identification of DEGs associated with the mitophagy-related Phenotype

2.8 Mitophagy score construction

2.9 Somatic gene mutation analysis

2.10 Immune Cell Infiltration Analysis and relation with mitophagy score

2.11 Genomic data and clinical information with Immunotherapy

2.12 Cell culture

2.13 Real-time quantitative PCR

2.14 Cell proliferation assay

2.15 5-Ethynyl-20-deoxyuridine (EdU) assay

2.16 Survival Analysis

2.17 Statistical analysis

3 Results

3.1 Clustering of ccRCC single cells

3.2 Communication analysis of single cell types

3.3 Expression and landscape of mitophagy marker genes in ccRCC

3.4 Classification of ccRCC patients based on mitophagy characters

3.5 Reconstruction of tumor clusters in ccRCC

3.6 Construction of mitophagy scoring model for quantification of mitophagy characters and prognosis evaluation in ccRCC

3.7 The relationship of clinicopathologic features and mitophagy score

3.8 Significance of mitophagy score for patient prognosis

3.9 HSPA8 act as antioncogenic effect in RCC

4 Discussion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1