Nineteen patients with ccRCC, diagnosed at Qilu Hospital of Shandong University, were randomly recruited and their clinical information was listed in Table S1 . Tumors and their matched NT specimens were collected from these patients who underwent nephrectomy. All the samples were stored in nitrogen tanks until use. The study was approved by the Institutional review board of Qilu Hospital of Shandong University (#KYLL-2021(KS)-192) and the signed informed consent was obtained from all patients.

RNA was extracted from primary tissues and cells using a RNAfast2000 kit (Fastagen) and quality control was performed using NanoDrop ND-1000 (Thermo Fisher Scientific). RNA sequencing was performed on 10 paired specimens. Sequencing libraries were generated using NEBNext^R Ultra™ RNA Library Prep Kit (New England Biolabs) according to the manufacturer’s recommendation. RNA sequencing was carried out using Illumina HiSeq 4000 sequencer at Metware Biotechnology (Wuhan, China). Paired-end reads were quality controlled by Q30 and Cutadapt software (v 1.9.3) was used to remove low-quality reads and 3’ adaptor-trimming. Hisat2 (v 2.0.4) was further used to align clean reads from sequencing, and sequencing depth and gene length were adjusted by Fragments Per Kilobase of transcript per Million (TPM) fragments mapped. The sequencing data were deposited in the GEO database (GSE217386).

The qPCR evaluation of CIN25 gene expression was performed on paired specimens from 9 patients with ccRCC. cDNA was synthesized using a PrimeScript™ RT reagent Kit (TAKARA). qPCR was carried out using SYBR Green of RT Master Mix (TAKARA) to assess mRNA levels of target genes based on 2(−ΔΔCT) values. β-actin mRNA levels were used as the internal control for normalization of target gene expression. All the primers were synthesized at Wuhan Genecreate Biotech (Wuhan, China) and primer sequences are listed in Table S2 .

The TCGA cohort of ccRCCs included 530 tumor samples with survival information available and 72 renal NTs. Patient clinical data were summarized in Table S3 (27). Transcriptome, mutation, copy number variations (CNAs) and clinical-pathological data were downloaded from https://gdc.cancer.gov/. One hundred and one patients with ccRCC were in the E-MTAB-1980 cohort (28), and RNA array and clinical information were downloaded from http://www.ebi.ac.uk. Patient clinical characteristics were listed in Table S4 . For RNA sequencing data, mRNA abundances were expressed as TPM. For array results (determined by 4×44K v2 microarray kit) from the E-MTAB-1980 cohort, probe-set values were used to quantify mRNA levels. ccRCC patients receiving Sunitinib treatments were contained in IMmotion150 ( Table S5 ) (29, 30) and IMmotion151 ( Table S6 ) trials (31, 32). Expression differences in CIN25-containing 25 genes were compared between ccRCC tumors and NTs in the TCGA cohort. For RNA expression, log2(TPM+1) based on RNA sequencing data was from https://gdc.cancer.gov/ as stated above. Protein expression data was obtained from Clinical Proteomic Tumor Analysis Consortium (http://ualcan.path.uab.edu/index.html).

The CIN25 gene signature includes the following genes responsible for spindle assembling/checkpoint, DNA damage checkpoint and cell cycle regulation: NCAPD2, ESPL1, CDK1, MELK, PRC1, KIF20A, TOP2A, TTK, TPX2, UBE2C, MCM7, MCM2, RFC4, FEN1, CDC45, FOXM1, RAD51AP1, H2AFZ, MAD2L1, PCNA, RNASEH2A, TGIF2, CCT5, TRIP13 and CCNB2 (22) ( Figure 1A ). The CIN25 score for each sample were expressed as mean Z-scores based on the Z-normalized mRNA level of 25 CIN-related genes above. We also calculated the CIN25 score based on single sample gene set enrichment analysis (ssGSEA) to confirm the accuracy of the Z-score method and other purposes.

Somatic CNAs were downloaded from https://xenabrowser.net/. CNA plots were made using R package ‘oncoPrint’ in ‘ComplexHeatmap’. Aneuploidy scores were the sum total of altered (amplified or deleted) chromosome arms. TMB is defined as the number of non-silent mutations per million bases and the data were downloaded from https://xenabrowser.net/.

Proliferation statuses were estimated using expression levels of Ki-67 mRNA and cell cycle scores, respectively. ccRCC cell cycle, stemness, EMT and angiogenesis signature scores were calculated based on ssGSEA or as the median z-score of genes included in each signature for each sample. These signatures are as follow: Angiogenesis: VEGFA, KDR, ESM1, PECAM1, ANGPTL4 and CD34 (33). Cell Cycle: CDK2, CDK4, CDK6, BUB1B, CCNE1, POLQ, AURKA, KI-67 and CCNB2 (34, 35). EMT: VIM, CDH2, FOXC2, SNAI1, SNAI2, TWIST1, FN1, ITGB6, MMP2, MMP3, MMP9, SOX10, GCS, CDH1, DSP and OCLN (36).

Telomere length data in the TCGA cohort of ccRCCs were from Bartheal et al. (37). Telomerase activity was evaluated using the telomerase score based on expression levels of the following 10 telomerase factors: TERT, TERC, DKC1, NHP2, NOP10, TCAB1, GAR1, NVL, RUVBL1 and RUVBL2 (38).

GSEA (http://www.gsea-msigdb.org/) analyses were performed to enrich KEGG pathways and hallmarks in two CIN25 subtypes of ccRCC tumors. P <0.05 and False discovery rate (FDR) <0.05 was considered statistically significant.

Cox regression analysis was conducted to assess the effect of the CIN25 cluster and clinical variables on survival. Then according to multivariate Cox regression analysis results, we constructed predictive nomograms including CIN25 and stage to predict 1-, 3-, and 5-year OS and/or PFS). Predicted survival of the nomogram against observed ones was plotted using the calibration curve. All nomograms and assessments of their predicative powers were made using R package regplot. Time-dependent Receiver Operator Characteristic (ROC) curves were used to determine sensitivity and specificity of OS and PFS predictions. Time-dependent ROCs and AUCs were made using Rpackage timeROC.

All statistical analyses were carried out using R package version 4.0.5. Wilcox and K-W sum tests were used for analysis of differences between two groups and among multi groups, respectively. Spearman’s Rank-Order Correlation coefficient was applied to determine correlation coefficients r between two variables. Survival analyses were made using log-rank test. The Survival and Survminer packages were employed to draw Kaplan–Meier survival curves for visualization of OS and PFS. Univariate and multivariate Cox regression analyses were used to determine the effect (HR and 95% CI) of various quantitative predictor variables on OS and PFS. P < 0.05 were considered as statistically significant.

Although aneuploidies and SCNAs have been well documented in ccRCCs, it remains unclear whether there exists the CIN25 signature as identified in other tumor types. We thus probed this issue first. RNA sequencing was performed on primary ccRCC tumors and their matched NTs from 10 patients who underwent nephrectomy. Expression levels of 25 genes in the CIN25 signature were evaluated in both tumors and NTs. As shown in Figure 1B , tumors exhibited significantly upregulated expression of 21/25 genes. The analysis of CIN25 ssGSEA in these samples further unraveled enhanced CIN25 levels in tumors ( Figure 1C ). For validation, we did qPCR-based expression analyses of these 25 genes in paired tumors and NTs from 9 patients, and largely similar results were obtained ( Figure 1D ). To confirm this finding obtained from our small patient cohorts, we analyzed the TCGA ccRCC sequencing data for their CIN25 signature expression. The comparison between 530 tumors and 72 NTs revealed significantly higher mRNA levels of 22/25 genes in tumors than in NTs ( Figure 1E ). Protein information was available in 20 of 25 genes, and protein levels were similarly higher in tumors, which is consistent with the transcriptomic analysis data ( Figure S1 ).

The results above demonstrate highly upregulated expression of almost all CIN25 genes in ccRCC tumors, however, a significant heterogeneity was observed among them. To determine whether ccRCC tumors could be classified based on the CIN25 expression score, we performed consensus cluster analyses of the TCGA cohort. Nonnegative matrix factorization clustering of CIN25 mRNA data showed consistency K = 2, indicating that a two-cluster classification was optimal ( Figure 2A ). In a total of 530 tumors, CIN25-cluster 1 (CIN25-C1, low CIN level) and cluster 2 (CIN25-C2, high CIN level) were 350 (66%) and 180 (34%), respectively ( Figure 2B ). Because the CIN phenotype is characterized by the presence of aneuploidy, we further compared global CNAs, and calculated aneuploidy, amplified and deleted scores between two CIN25 clusters ( Figures 2C, D ). Indeed, the aneuploidy score was significantly higher in CIN25-C2 tumors (CIN25-C1 vs CIN25-C2, P = 1.78E-04) ( Figures 2C, D ). Interestingly, the amplified score was more robustly higher in the CIN25-C2 tumors than in CIN25-C1 ones (CIN25-C1 vs CIN25-C2, the amplified and deleted scores, P = 2.86E-18 and 4.95E-02, respectively) ( Figure 2D ). Moreover, we also calculated CIN25 ssGSEA score of each tumor based on the expression of 25 genes and observed a drastically higher CIN25 ssGSEA score in CIN25-C2 tumors ( Figure 2E ). To validate the CIN25 clustering classification of ccRCC tumors, we carried out the same analysis of the E-MTAB1980 ccRCC cohort, and tumors were readily categorized into two distinct CIN25 clusters, with higher CIN25 ssGSEA scores in CIN25-C2 tumors ( Figures 2F, G ).

We next determined the potential association between CIN25 subtypes and clinic-pathological variables in ccRCC tumors. We first examined the distribution of two clusters between two genders and different age groups (≥60 and <60 years) in the TCGA cohort and did not observe significant differences, although male patients had a slightly higher frequency of CIN25-C2 than did females (38.6% vs 29.4%, P = 0.055) ( Figure 3A ). CIN25-C2 was more frequently observed in higher-stage (P = 5E-06) and higher-grade tumors (P = 0.007) ( Figure 3B ). Very similar results were obtained from the analysis of the E-MTAB1980 cohort ( Figures 3C, D ).

We further performed the same analysis of 10 ccRCC patients whose tumors were with transcriptomic profiling. Because 10 tumors were too few to make a CIN classification, we calculated ssGSEA score to express CIN25 levels in each tumor. The CIN25 ssGSEA score was significantly increased in higher-stage (III/IV vs I/II, P = 0.019) and grade (III/IV vs I/II, P = 0.032) tumors ( Figure 3E ), which was consistent with the result obtained from the TCGA patient analysis.

It is well established that telomere dysfunction drives CIN in oncogenesis (39). We thus sought to determine whether telomere length was altered in the TCGA ccRCC cohort. Telomeres were significantly shorter in tumors than in matched NTs ( Figure 4A ), but there was no statistically significant difference in telomere length between CIN25-C1 and C2 subtypes ( Figure 4A ). Because telomeric DNA is synthesized by telomerase, while telomerase activity is primarily governed by its catalytic component telomerase reverse transcriptase (TERT) (40), we further compared TERT expression and telomerase activity between CIN-C1 and C2 tumors. As shown in Figure 4B , TERT mRNA levels were significantly higher coupled with the increased frequency of TERT copy number gain in CIN25-C2 tumors (C2 vs C1 for TERT mRNA and copies: P = 1.84E-08, and 0.018, respectively). Telomerase activity, as determined using telomerase score (38), increased markedly in the CIN25-C2 tumors compared with that in CIN25-C1 tumors (P = 2.15E-05) ( Figure 4C ). Moreover, there was a significantly positive correlation between telomerase and CIN25 ssGSEA scores (R = 0.43, P <2.22E-16) ( Figure 4C ).

CIN is one subtype of genomic instability, whereas the later also includes several other forms of genomic aberrations such as nucleotide instability (NIN), microsatellite instability (MSI), homologous recombination deficiency (HRD), etc. (41). Thus, we further addressed the association of CIN25 clusters with the following important alterations: (i) Tumor mutation burden (TMB) (P = 0.034) ( Figure 4D ). Moreover, we compared the top 10 mutated genes between two subtypes. As expected, VHL, PBRM1, BAP1, MTOR and SETD2 are among the top mutated genes in both subtypes, however, significantly higher mutated frequencies of BAP1 and SETD2 were observed in the CIN25-C2 tumors (BAP1 and SETD2: P = 0.0003 and 0.018, respectively ( Figure S3 ). In addition, KDM5C mutation was more frequent in the CIN25-C1 tumors (P = 0.04). (ii) Intratumor genetic heterogeneity (P = 0.01) ( Figure 4E ). (iii) HRD (P = 0.0002) ( Figure 4F ). (iv) MSI (P = 0.00004) ( Figure 4G ). CIN25-C2 tumors exhibited significantly higher levels or frequencies of all the aberrations analyzed above.

We then wanted to assess whether this CIN25 classification system could predict patient survival. The TCGA cohort of 530 ccRCC patients was first evaluated as the discovery set. Log-rank test analysis unravelled that those patients in the CIN25-C2 group had significantly shorter OS and PFS, as shown by Kaplan–Meier survival curves (P = 7.57E-06 and 4.83-07 for OS and PS, respectively) ( Figure 5A ). We further performed univariate COX regression analyses by including patient age, gender, stage, and grade together with the CIN25 clustering system. Advanced Stages, higher grades and CIN25-C2 were all associated with shorter OS and PFS ( Figures 5B, C ). Multivariate COX regression analyses showed that all three of them were independent prognostic factors for shorter OS and PFS ( Figures 5B, C ).

The E-MTAB-1980 ccRCC cohort as the validation set were further analyzed in the same manner. There was no PFS information available, and we only evaluated OS. Kaplan–Meier survival analysis showed that CIN25-C2 was associated with significantly shorter OS (P = 0.0003) ( Figure 5D ), and the CIN25 subtype and stages were independent OS predictors, as assessed using univariate ( Figure 5E ) and multivariate COX regression analyses ( Figure 5F ).

The data above consistently show that CIN25-C2 subtype and advanced stages are independent prognostic variables for OS and/or PFS in both TCGA and E-MTAB-1980 cohorts. We thus established a prognostic nomogram composed of CIN25 subtypes and stages. For the TCGA cohort, the nomograms exhibited a highly accurate estimation of OS and PFS possibilities at 1-, 3- and 5-years ( Figures S2A, B ). Similar results were obtained for OS prediction in the E-MTAB-1980 cohort ( Figure S2C ). To further evaluate the sensitivity and specificity of their prediction, we did time-dependent ROC analyses. In the TCGA cohort, Area under ROC curves (AUCs) for 1-, 3- and 5-year OS were 0.799, 0.767 and 0.740, respectively, while the AUCs for 1-, 3- and 5-year PFS were 0.825, 0.797 and 0.798, respectively ( Figures S2D, E ). For 1-, 3- and 5-year OS in the E-MTAB-1980 cohort, AUCs were 0.886, 0.871 and 0.838, respectively ( Figure S2F ). Separate analyses of these two variables showed that AUCs were largely between 0.6 and 0.7, between 0.7 and 0.8 for all CIN25- and stage-based predictions of 1-, 3- and 5-year survival ( Figures S2G, H ). AUCs obtained from stage-prediction were bigger in all the estimations.

Sunitinib has long been applied for advanced ccRCC treatment as the first line drug (14), however, reliable biomarkers to predict its efficacy or patient response are few (6, 8). We thus determined whether the CIN25 cluster classification could help distinguish responders from non-responders in patients treated with Sunitinib. Toward this end, the IMmotion151 cohort of 416 ccRCC patients treated with Sunitinib was first analyzed as the discovery set (31, 32). Patient response to Sunitinib was categorized into complete remission (CR), partial remission (PR), stable disease (SD) and progressive disease (PD), respectively. A total of 416 patients were classified into CIN25-C1 (273 patients) and C2 (143 patients) groups. The CR and PR (CRPR) rate was 42% and 26% in CIN25-C1 and C2 groups, respectively (P = 0.0004) ( Figure 6A ). The disease progression during the Sunitinib treatment occurred in 14.9% and 30.9% for CIN25-C1 and C2 patients, respectively. The median PFS for C1 and C2 patients was 5.6 and 11.2 months, respectively (P = 7.78E-08; HR, 1.90 (95% CI: 1.45 – 2.47) ( Figure 6B ). We then analyzed the IMmotion150 cohort (29, 30) to validate the findings obtained from IMmotion151 cohort. In a total of 85 available patients, CIN25-C1 and C2 were 58 and 27, respectively. The total CRPR rate was 41.4% and 14.8% in CIN25-C1 and C2 groups, respectively (P = 0.002) ( Figure 6C ). Almost the half of CIN25-C2 patients (48.1%) underwent progression during the treatment, while only 13.8% of CIN25-C1 patients did so ( Figure 6C ). Higher CRPR rates in CIN25-C1 group led to longer PFS, and the median PFS for C1 and C2 patients was 4.4 and 9.8 months, respectively (P = 0.002; HR, 2.13 (95% CI: 1.18 – 3.84) ( Figure 6D ).

We next performed the GSEA analysis to probe differences in signaling pathways between two tumor groups. Figures 7A, B showed significantly enriched KEGG and hallmark pathways in CIN25-C2 tumors, and almost all of them are oncogenic and play key parts in ccRCC development and progression, such as G2/M checkpoint, E2F and MYC targets, IL6-JAK-STAT3, glycolysis, EMT and others ( Figure 7C ). Consistent with these enriched pathways, CIN25-C2 tumors had robustly strong proliferation activity compared to CIN25-C1 tumors, as assessed using proliferation marker Ki-67 and cell cycle score, and stemness score ( Figure 7D ); furthermore, an established EMT 16 gene signature (36) was further used to evaluate EMT between CIN25-C1 and C2 tumors and significantly increased EMT scores were observed in the CIN25-C2 group ( Figure 7D ) (P = 0.035).

EZH2, a histone methyltransferase catalyzing H3K27 trimethylation (H3K27me3), has been shown to promote stemness, EMT and Sunitinib resistance in ccRCC and other tumors (42–44). Given the results above, we set out to determine whether EZH2 expression differed between CIN25-C1 and C2 tumors. The analysis of both TCGA and E-MTAB1980 ccRCC cohorts showed robustly higher EZH mRNA levels in CIN25-C2 than in C1 tumors (CIN25-C1 vs C2: P = 2.21E-38 and 3.12E-06, respectively) ( Figure 8A ). In the Sunitinib-treated IMmotion150 and IMmotion151 cohorts, similar results were obtained (CIN25-C1 vs C2: P = 9.40E-08 and 1.71E-27 for IMmotion150 and 151, respectively) ( Figure 8B ). We further compared differences in EZH2 expression between responders and non-responders to Sunitinib. As expected, tumors from resistant patients expressed significantly higher levels of EZH2 than did those from responders (P = 0.021 and 0.004, respectively) ( Figure 8C ).

In addition to higher EZH2 expression, poor angiogenesis is also a well characterized predictor for Sunitinib resistance (34, 35), and we thus analyzed the angiogenesis score in CIN25 subtypes of ccRCC tumors. As shown in Figure 8D , a significantly lower angiogenesis score was observed in the CIN-C2 tumors from the IMmotion151 (CIN25-C1 vs C2: P = 9.78E-17) and IMmotion150 (P = 0.0096) cohorts of ccRCC patients treated with Sunitinib. The TCGA and E-MTAB1980 ccRCC analyses showed similar results, which validated the observations above ( Figure 8E ).