Quantitative Proteomics Identifies Secreted Diagnostic Biomarkers as well as Tumor-Dependent Prognostic Targets for Clear Cell Renal Cell Carcinoma

Renal cell carcinoma (RCC) is among the 10 most common cancers in both men and women and causes more than 140,000 deaths worldwide every year (1). One third of the patients already display metastasis at the time of diagnosis and are expected to have a median survival of around 26 months (2, 3). Clear cell RCC (ccRCC) is the most common histologic subtype accounting for more than 70% of the cases (3). Traditional treatment with chemotherapy or radiotherapy fail to treat this malignancy, limiting the treatment options primarily to surgical excision, together with targeted therapy and/or immunotherapy in suitable cases (1). Biomarkers can be used as indicators of tumor progression and patient outcome. In recent years, their discovery has largely benefitted from technological advances in high-throughput mass spectrometry (MS)-based proteomics. However, despite their urgent need, no universal biomarkers of clinical utility are currently known for ccRCC (4). This is due to previous studies lacking in coverage depth with generally less than 2,000 protein identifications and presenting insufficient validation of the candidates by multiple independent methods or in independent cohorts (5–11). Also, most candidates are restricted to the tissue context, however, biomarkers for easily available and noninvasive biospecimens such as body fluids are of more practical utility in the clinical setting (12, 13). Urine has the advantage over blood to be collectable in large amounts, to be less prone to proteolytic degradation and to be less complex, however, reliable urine biomarkers are still scarce (12, 14–16).

Over the past years, it also became clear that a “one-size-fits-all” treatment strategy for patients with RCC is insufficient for the effective management of the disease and the prevention of recurrence. This is particularly due to genetic and epigenetic intrapatient tumor heterogeneity caused by branched clonal evolution of the cancer cells leading to poor treatment response or to treatment resistance (4, 17). Another crucial factor causing tumor diversity is the composition of the tumor microenvironment, which consists of the extracellular matrix and different types of stromal and immune cells (18). In recent years, patient classification based on in-depth molecular expression profiles rather than histologic parameters have grown in importance. Drug treatment according to individual molecular aberrations opens new avenues for a tumor-oriented intervention. Especially in the case of grade 2 and 3 tumors, the disease progress and outcome can vary dramatically (19), and therefore be difficult to predict, which further necessitates appropriate tumor stratification. In the past years, several genomics- and transcriptomics-based classification systems have been proposed for ccRCC (20), such as the widely investigated microarray-derived ccA/ccB clustering signature, consisting of 34 marker genes (ClearCode34; refs. 19, 21). This clustering signature revealed noticeable prognostic difference between the low-risk ccA tumors, characterized particularly by fatty acid metabolism and angiogenesis, and the high-risk ccB tumors, associated with epithelial-to-mesenchymal transition (EMT) and Wnt signaling. In contrast, classification systems of RCC based on protein expression data have rarely been described previously (20, 22). Only one proteomics study describes a clustering signature which distinguishes ccRCC tissues into three clusters (ccRCC1–3) with diverse malignancy and immune cell infiltration (22).

Here, we performed global quantitative shotgun proteomics with ccRCC tissues and validated potential biomarker candidates by taking orthogonal approaches. We also discovered a novel secreted biomarker in urine samples. We further developed a two-cluster based classification system for the risk stratification of patients with ccRCC using the generated proteome profiles.

The study proposal was approved by the Ethics Committee of Koc University (Istanbul, Turkey) in September 2017 (no. 2017.145.IRGB2.051). The approval was prolonged for 2 more years and was expanded for urine sample collection in 2019. Tissue and urine samples were collected with the patients’ informed written consent following the guidelines of the Declaration of Helsinki.

Tumor and normal adjacent tissue (NAT) samples were collected during partial or radical nephrectomy in two hospitals in Istanbul (Koc University School of Medicine and Istanbul University Istanbul Faculty of Medicine). Samples were immediately stored at −80°C after collection. The age of the patients varied between 28 and 74 and reflected a typical gender representation for RCC with a 2.6:1 male-to-female-ratio. Urine samples were collected preoperatively from patients with ccRCC and a tumor-free control group in the Urology Department of the Koc University School of Medicine (Istanbul, Turkey). Samples were supplemented with complete EDTA-free protease inhibitor and immediately stored at −80°C. Detailed clinical information about the patient cohort is summarized in Supplementary Data S1.

Normal adjacent and tumor tissue samples from 13 patients were subjected to protein isolation using 8 mol/L urea, 1 mmol/L sodium orthovanadate, complete EDTA-free protease inhibitor mixture (0.5 tablet for 5 mL, Thermo Fisher Scientific), phosSTOP phosphatase inhibitor mixture (0.5 tablet for 5 mL, Roche), and 1% n-octylglucoside in 50 mmol/L ammonium bicarbonate (ABC, pH 8.0). For extensive protein extraction, 20–30 mg of tissue samples were homogenized with the Bullet Blender Tissue Homogenizer (Next Advance) using 0.5 mm zirconium oxide beads. The cells were further disrupted by passing the sample through a fine syringe needle (26G) 10 times on ice. The protein concentration was measured using the BCA assay (Pierce, Thermo Fisher Scientific), followed by reduction of disulfide bonds with 10 mmol/L dithiothreitol for 1 hour at 56°C and cysteine alkylation with 20 mmol/L iodoacetamide for 45 minutes in the dark at room temperature. The urea concentration was reduced to 1 mol/L with 50 mmol/L ABC and trypsin was added at 1:50 enzyme:protein ratio. The protein digest took place at 37°C overnight and was quenched by acidification using 10% formic acid (FA). The protein samples were desalted on SepPak C18 cartridges (Waters) using 0.1% FA for washing (RP A solution) and 0.1% FA in 80% acetonitrile (ACN) for elution (RP B solution).

Peptide samples were labeled at their primary amines with stable dimethyl isotopes according to the protocol described by Boersema and colleagues (23). For light labeling 4% formaldehyde solution (CH₂O) and 0.6 mol/L cyanoborohydride solution (NaBH₃CN) were mixed in 50 mmol/L sodium phosphate buffer (pH 7.5), while for heavy labeling 4% ¹³CD₂-labeled formaldehyde solution (¹³CD₂O) and 0.6 mol/L cyanoborodeuteride solution (NaBD₃CN) were mixed. For 7 patients, the NAT samples were light labeled while the tumor samples were heavy labeled. For the remaining 6 patients, the labeling was swapped to eliminate bias from the label choice. A total of 600–700 μg of peptides were loaded on 1cc SepPak C18 cartridges and flushed with the prepared light and heavy labeling solution, respectively. The labeled peptides were washed twice with RP A solution and eluted with RP B solution. The peptide mixtures were then fractionated on-column by SCX chromatography. For this purpose, the dried peptide samples were reconstituted in SCX loading buffer (7 mmol/L KH₂PO₄, 30% ACN, pH 2.65) and light and heavy labeled peptides were mixed at 1:1 ratio based on their median peptide intensities. The SCX resin packed cartridge (HyperSep, Thermo Fisher Scientific) was prepared by sequential addition of following solutions: methanol, elution buffer F9, HPLC-grade water, equilibration buffer (50 mmol/L K₂HPO₄, 500 mmol/L NaCl, pH 7.5), HPLC-grade water and loading buffer. The bound peptides were eluted by stepwise addition of following elution buffers which consisted of the loading buffer with increasing KCl concentrations (F1: 30 mmol/L; F2: 45 mmol/L; F3: 60 mmol/L; F4: 75 mmol/L; F5: 90 mmol/L; F6: 105 mmol/L; F7: 120 mmol/L; F8: 150 mmol/L; and F9: 500 mmol/L). The obtained 10 fractions per patient (flow-through included) were extensively desalted using 1cc SepPak columns.

The collected SCX fractions were reconstituted in 5% FA and 5% ACN (MS analysis solution) and analyzed in duplicate with 120-minute linear gradients on an UltiMate 3000 RSLCnano reversed phase chromatographic platform (Thermo Fisher Scientific) coupled to a Q Exactive hybrid quadrupole-Orbitrap mass spectrometer (Thermo Fisher Scientific). Approximately 300 ng of each fraction was loaded onto an in-house packed 100 μmol/L i.d. × 17 cm C18 column (Reprosil-Gold C18, 5 μmol/L, 200Å, Maisch) and run with a flow rate of 300 nL/minute. The chromatographic separation of the peptides started at 4% of solution B (0.1% FA in ACN) and gradually increased to 25% in 69 minutes. The gradient continued from 25% to 40% of solution B in the next 20 minutes. Peptides in the mass range 400–1,500 m/z and with a positive polarity were allowed for detection in data-dependent acquisition (DDA) mode. For the MS1 spectra acquisition, the resolution was set to 70,000, the automatic gain control (AGC) target to 1e6 and the maximum injection time to 32 ms. The top 15 most intense peptides per cycle were selected for fragmentation in the higher energy collisional dissociation cell with a normalized collision energy (NCE) of 26. MS2 spectra acquisition was conducted at a resolution of 17,500, an AGC target of 1e6, a maximum injection time of 85 ms and a fixed first mass of 120 m/z. Furthermore, the isolation window was set to 2.0 m/z, the dynamic exclusion was set to 35 seconds and the charge exclusion was set as unassigned, 1, 6–8, >8.

For the targeted proteomics approach, unlabeled digested tissue samples of patients were analyzed in parallel reaction monitoring (PRM) mode. Samples were reconstituted in MS analysis solution and run with a 90-minute linear gradient on a similar system as described above except that a Q Exactive HF mass spectrometer was engaged. Approximately 600 ng were loaded onto an in-house packed 100 μmol/L i.d. × 25 cm C18 column (Reprosil-Gold C18, 3 μmol/L, 200 Å, Maisch) and run with a flow rate of 300 nL/minute. The gradient started at 4% of solution B, increased to 25% in 49 minutes and continued from 25% to 45% of solution B in the next 15 minutes. Peptides with a positive polarity were targeted. MS2 spectra acquisition was performed at a resolution of 45,000, an AGC target of 2e5, a maximum injection time of 100 ms and a fixed first mass of 100 m/z. Moreover, the isolation window was set to 1.2 m/z and the NCE to 28.

Peptide identification and quantification was done with Proteome Discoverer (PD; v1.4, Thermo Fisher Scientific) using the SequestHT search engine. Peptide spectral matches were searched against a Swissprot database containing 21,039 entries for Homo sapiens retrieved from Uniprot in March 2016. Trypsin was selected as hydrolytic enzyme with a maximum number of allowed missed cleavages of 2. Regarding peptide identification, a mass tolerance of ±20 ppm for precursor masses and ±0.05 Da for fragment ions were selected. For quantitation, the 2-plex dimethyl heavy/light (H/L) method with a mass precision of 2 ppm for precursor measurements was used. Light and heavy dimethylation of peptide N termini and of lysine residues, as well as methionine oxidation, were set as dynamic modifications. Cysteine carbamidomethylation was set as fixed modification. The false discovery rates (FDR) for peptide and protein identifications were set to 1%. Only peptides with medium or high identification confidence, with a sequence length between 7 and 25 and with a peptide rank of minimum 1 were allowed.

Quantification ratios of samples with labeling swap were converted to H/L format. H/L ratios for biological replicates were average values from two technical replicates. The obtained data were filtered for proteins quantified in at least 1 out of the 13 patients (5,799 proteins), which are hereafter denoted as “quantified proteins.” The H/L values were sample-wise normalized by the sample median and log₂ transformed.

First, the pH value of urine samples was adjusted to pH 8.0 using 500 mmol/L ABC supplemented with EDTA-free protease inhibitor and then, the biospecimen was filtered through an Amicon Ultra-0.5 mL Centrifugal Filter Device (10 kDa, Merck Millipore). Collection of the retentate (proteins in ∼100 μL 50 mmol/L ABC) was performed by flipping the filter device upside-down and centrifugation at 1,000 × g for 2 minutes. The protein concentration was determined using the BCA assay.

For immunoblot analysis, 40 mg of normal tissue and 70–80 mg of tumor tissue were lysed in 0.1% Triton-X in 1× PBS and EDTA-free protease inhibitor. In total, 40 μg of tissue samples and 20 μg of urine samples were run on 8% to 10% SDS gels and transferred to nitrocellulose membranes for 3 to 4 hours at room temperature. Primary antibody incubation was done overnight at 4°C with following dilutions: anti-PLOD2 mouse at 1:500–1:1,000 (MAB4445, R&D Systems), anti-FERMT3 mouse at 1:500–1:1,500 (ab173416, Abcam), anti-SPARC mouse at 1:1,000 (33–5500, Invitrogen), anti-SIRPα rabbit at 1:1,500–1:3,500 (13379, Cell Signaling Technology) and anti-ACTB mouse at 1:5,000 (ab6276, Abcam). As secondary antibodies horseradish peroxidase (HRP)-conjugated anti-mouse IgG at 1:2,000 (7076S, Cell Signaling Technology) and HRP-conjugated anti-rabbit IgG at 1:2,000 (7074S, Cell Signaling Technology) were used. Protein expression was visualized using Clarity Western ECL Substrate (Bio-Rad). Densitometric analysis was performed with ImageJ (v1.46r), and visualization and statistical analysis (Student t test P < 0.05) with GraphPad PRISM v8, respectively.

Selected unlabeled peptides from digest experiments were targeted by PRM as quantifiable surrogates for the proteins of interest. Candidate peptides were selected as follows: uniqueness of peptides was ensured by BLASTp analysis, peptides with length between 7 and 18aa were preferred but peptides with unstable charge state across samples, with ragged ends, with missed cleavages and/or with posttranslational modifications, respectively, were avoided. When no alternative peptides were available for a protein, exceptions to these exclusion criteria were made. No exception was made for the charge state criteria. MS/MS spectra were inspected extensively and peaks with a mass window <10 ppm were integrated by manually setting the borders in Skyline (v19.1). The spectral library consisted of 60,833 spectra from DDA-based data acquisition of tumor and normal digest samples of the discovery cohort. Doubly and triply charged peptides in the range of 375–1,500 m/z with a library ion match tolerance of 0.5 m/z were accepted. Only spectra with a high correlation with the spectral library (dotp value >0.8) were considered. For each peptide, the summed peak area of its six fragment ions was calculated and then normalized by the summed peak area of the β-actin peptide GYSFTTTAER in the sample. Calculated tumor/normal peak area ratios were transformed to log₂.

IHC analyses were performed with the BenchMark Ultra autostaining platform (Ventana Medical Systems). Briefly, the tissue paraffin sections were cut to 3 μm thickness onto charged slides. Sections were deparaffinized and rehydrated through an alcohol series. Tissue sections were incubated with primary antibodies at 37°C with following dilutions: anti-SPARC mouse at 1:800 (33-5500, Invitrogen), anti-SIRPα rabbit at 1:25 (13379, Cell Signaling Technology), anti-FERMT3 mouse at 1:2,500 (ab173416, Abcam) and anti-PLOD2 mouse at 1:250 (MAB4445, R&D Systems). The UltraView Detection kit (760–500, Ventana Medical Systems) was used for detection. The reaction product was visualized with 3, 3′-diaminobenzidine chromogen and counterstained with hematoxylin. All stained slides were evaluated in a blinded fashion by a single pathologist. The intensity was scored as mild (score 1: any positivity that could be seen at high magnification), moderate (score 2: any staining between score 1 and 3), and marked (score 3: any staining that could be seen at low magnification). The prevalence of the staining was evaluated as score 0 (<5%), score 1 (5%–25%), score 2 (25%–75%), and score 3 (>75%).

Significantly dysregulated proteins between tumor and normal tissues were determined using the one-sample Wilcoxon signed-rank test. For this purpose, the Python function scipy.stats.wilcoxon (v0.14.0) was extended by the feature to omit missing quantification values and was then applied genewise on the data. Benjamini—Hochberg (BH) adjusted P < 0.05 were considered statistically significant. For the determination of statistically significant cluster-discriminative proteins, the two-sample Student t test was applied. P values were calculated genewise using the ttest_ind function of the scipy.stats module by omitting missing values (nan_policy = “omit”). P < 0.05 were considered statistically significant.

Significantly regulated proteins were annotated by their role in cancer according to the COSMIC Cancer Gene Census database (v91). For Gene Ontology (GO) cellular component and for pathway annotations, over-representation analysis was applied using the WebGestalt platform (v2019; ref. 24) with a BH-adjusted FDR of 0.02 and 0.05, respectively. Pathway annotations were combined analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, Hallmark50, and Wikipathway_cancer databases. For GO biological process annotation of the significantly regulated proteins, the PANTHER (v14) platform (25) was used with an FDR setting of 0.05.

Optimal subnetworks, which represent the list of significant proteins best, were reconstructed by integrating a reference human interactome with our proteomic data using the Forest module of Omics Integrator (v0.3.1; ref. 26) and the visualization tool Cytoscape (v3.8.0; ref. 27). The iRefWeb_2013 interactome file was used as global network reference. Nodes were annotated by GO cellular compartment using QuickGO (28), which was limited to Uniprot assignments. For multiple allocations, priority was given as follows: “extracellular”, “membrane”, “nucleus”, “cytoplasm”, “other”, or “unknown”. In addition, networks were annotated with GO biological process or Reactome pathways, respectively, using the Cytoscape plug-in ClueGO (v2.5.7; ref. 29). For cluster-discriminative proteins, the network was complemented with associated inhibitory FDA-approved drugs retrieved from the Drugbank database (v5.1.7, access August 2020; ref. 30). Chemical elements were excluded as drugs.

Unsupervised hierarchical clustering was applied sample-wise on the global proteome data. A distance matrix was prepared by calculating the square root of the sum of squared differences between samples (Euclidean distances) with an in-house generated Python script. The subsequent clustering was done using the linkage function of the scipy.hierarchy module with the weight method “ward” and the distance metric “Euclidean”. Furthermore, k-means clustering of the 5,799 quantified proteins was applied using the tool ConsensusCluster with 300 iterations (31). Principal component analysis (PCA) was employed for the first two components using the PCA function of the sklearn.decomposition module.

Two independent expression datasets of ccRCC were used for comparison with our study. Fragments per kilo base per million mapped reads (FPKM) normalized RNA-seq data of The Cancer Genome Atlas (TCGA)-kidney renal clear cell carcinoma (KIRC) cohort consisting of 539 tumor and 72 normal tissue samples was retrieved from the TCGAPortal (http://tcgaportal.org/download.html, access April 2020). Genes with more than 70% missing expression values were filtered out and tumor counts were genewise normalized to the upper quartile of respective normal counts and log₂ transformed. Moreover, the Clinical Proteomic Tumor Analysis Consortium (CPTAC) ccRCC proteome data was retrieved from the CPTAC portal (https://cptac-data-portal.georgetown.edu, access November 2019). Quality control samples and the reported contaminated sample (C3N-00314) were excluded. log₂ tumor/normal ratios were calculated for the 80 tumors with adjacent normal sample by eliminating the relation to the TMT reference sample. Pearson and Spearman correlations of median log₂ tumor/normal ratios were calculated using the respective Python scipy.stats functions.

For evaluation of the prognostic impact of the tumor clusters, the follow-up data of the TCGA-KIRC cohort for 534 patients were used (Firehouse Legacy, downloaded from cBioPortal, access August 2020). Survival analysis was performed using the KaplanMeierFitter and multivariate_logrank_test functions of the Python library lifelines (v0.25.0). Furthermore, statistical significance of the frequency of clinicopathologic factors in the clusters was determined using the fisher_exact function of the scipy.stats library. To determine the impact of clinicopathologic factors on survival, univariate and multivariate Cox regression analyses were performed using the CoxPHFitter function of the lifelines library.

To verify the observed clustering signature, TCGA and CPTAC datasets were subjected to two-sample Student t test after clustering (P < 0.05). Cluster assignment was then conducted on the basis of highest overlap of the significantly different proteins between the clusters of the validation datasets with those between the clusters of our proteomics study. To identify the key discriminatory elements of the clustering signature (marker genes), a cut-off value of ± 2 for the t-test difference and a cutoff >2 for the −log₁₀P value were set.

In this study, a comprehensive quantitative proteome analysis of ccRCC has been conducted to address the scarcity of reliable biomarker panels and to uncover underlying molecular phenotypes. To this end, dimethylation-based quantitative proteomics analysis with extensive homogenization of biospecimens has been performed with 13 ccRCC patient’s frozen tumor and normal adjacent tissue samples (Fig. 1A). The discovery cohort included patients of different age, grade, stage, and gender (Supplementary Data S1). We identified a total of 10,160 unique proteins from all samples with a FDR of 0.01. Protein identification and quantification were robust across patients with a coefficient of variation of approximately 6.7% for the identifications and 8.5% for the quantifications, respectively (Fig. 1B). On average, 4,453 proteins were identified and 3,118 proteins were quantified in the paired samples, respectively. Abundance levels of the ccRCC proteome encompassed almost 6 orders of magnitude in dynamic range and correlated with diverse cellular compartments (Fig. 1C). The most abundant proteins were associated with angiogenesis such as blood coagulation and platelet degranulation, and with translational processes such as chaperone activity and protein maturation in the endoplasmic reticulum. Notably, the low-abundance range was characterized by proteins predominantly localized to the mitochondrial matrix and mitochondrial inner membrane. In general, the distribution of quantified proteins followed a near normal pattern (Fig. 1D, top histogram) over more than 14 orders of magnitude with a noticeable skewness toward positive fold change values suggesting more pronunciation of overexpressed processes in ccRCC tumors compared with normal tissues (Fig. 1D and E).

In total, 5,799 of the 10,160 identified proteins were quantified by dimethyl labeling (Fig. 2A). We determined 955 proteins as statistically differentially expressed between tumors and NATs (BH-adjusted P < 0.05) with 409 significantly upregulated and 546 significantly downregulated proteins. Among these were 44 genes associated with mutations that are causally implicated in cancer, such as oncogenes, tumor suppressor genes and fusion genes, according to the COSMIC Cancer Gene Census database (Supplementary Fig. S1A and S1B). Significantly elevated oncogenes were EGFR, H3F3A, AKT1, MAPK1, SRC, CALR, and XPO1 (Fig. 2B; Supplementary Data S2). The activation of the EGFR-MAPK1 axis has previously been reported to enhance cell proliferation in RCC (22, 32). Also, our data implicated that the EGFR-PI3K-AKT1 and AKT1 related mTOR signaling pathways might be activated in ccRCC, which are both associated with tumor growth and survival (33). Another significantly upregulated oncogene is the tyrosine kinase SRC, which has previously been reported to be recruited by EGFR and to affect its downstream players RAS and AKT1 in RCC (34). Despite its categorization as a tumor suppressor, fatty acid synthase was highly upregulated in the tumor tissues relative to NATs, which has previously been associated with tumor aggressiveness and poor patient outcome in RCC (35). We inferred that several oncogenes possibly mediate tumor proliferation and survival in ccRCC tumors, emphasizing the need for further investigation especially of the EGFR-related molecular mechanisms in respect of alternative directed therapeutic intervention.

As many cancers, RCC undergoes substantial metabolic reprogramming by increasing glycolytic processes and diminishing mitochondrial ATP production to meet its excessive energy requirements for tumor growth and progression (35, 36). Our findings for the significantly altered biological processes were in agreement with the reported metabolic pathobiology of RCC. The significantly upregulated proteins were primarily associated with protein maturation and transport, PTM events such as glycosylation and phosphorylation of macromolecules, glycolysis-related processes, and immune response (Fig. 2C). In contrast, the significantly downregulated proteins were mainly involved in mitochondrial events such as oxidative phosphorylation, electron transport chain, proton transport, and tricarboxylic acid (TCA) cycle (Fig. 2D).

To unveil the most prominent cross-talks between significantly upregulated proteins, we constructed a protein-protein interaction subnetwork using Omics Integrator (26), which confirmed that the highly elevated oncogenes EGFR and SRC were central players in ccRCC tumorigenesis (Fig. 2E; Supplementary Fig. S2). While most of the hubs were associated with intracellular transport of macromolecules (SRC, RPS4X, SUMO1, and YWHAZ) and organelle organization (UBQLN4; Supplementary Fig. S2), EGFR was associated with immune system–related processes and molecule secretion (Fig. 2E). SUMOylation modifications have previously been reported to be elevated in a subgroup of low-grade ccRCC tumors, which were particularly associated with angiogenesis and absence of infiltrating immune cells (22). Furthermore, the regulatory protein YWHAZ (14–3-3ζ) has previously been noted to play a promoting role in primary as well as metastatic ccRCC (7, 8, 11).

Given that secretion of macromolecules is essential for sculpting the tumor microenvironment to promote tumor growth and metastasis, and that vesicle-mediated transport and secretion were one of the significantly enriched biological processes in the ccRCC tumors (Fig. 2), we sought to identify and characterize novel secreted biomarkers. To this end, we selected four biomarkers among the top 70 highly upregulated proteins predicted to be located in the “extracellular space.” To the best of our knowledge, the selected biomarkers, namely, PLOD2, FERMT3, SPARC, and SIRPα, have previously not been suggested or validated as ccRCC biomarkers in proteomics-based studies. According to the CGC annotation, SIRPα has a tumor suppressor role in cancer (Fig. 2B), however, in this study, SIRPα was with a median log₂ tumor/normal ratio of 2.15 clearly an upregulated protein in the ccRCC tumor tissues.

Immunoblot analysis confirmed the elevated expression of all four biomarkers in the tumors compared with the adjacent normal tissues (Fig. 3A). This was observed in grade 2, grade 3, as well as grade 4 patients. Next, we attempted to validate the biomarkers by targeted proteomics in digested unlabeled tissue samples. PRM allowed tracking of selected peptides, which served as surrogates for the proteins of interest, and whose summed product ion peak areas were used to calculate tumor/normal ratios (Fig. 3B). Besides the four selected biomarkers, six further putatively secreted proteins of interest were monitored via PRM (Fig. 3C and D). As expected, the selected 20 targets were more abundant in tumor samples compared with NATs and had comparable tumor/normal ratios across grades. We observed a significant positive Pearson correlation of 0.84 (P < 0.001; Supplementary Fig. S1C) between the discovery cohort (Fig. 3C) and an independent cohort consisting of 11 additional patients with ccRCC (Fig. 3D) of different grades, which supported the suitability of these proteins as potential biomarkers.

Next, we tested the impact of intratumoral heterogeneity on the expression of the biomarkers in random parts of tumor and adjacent normal tissues with morphologic diversity from 3 selected patients (Fig. 3E; Supplementary Fig. S3). We could not observe any remarkable local distortion of abundance for the biomarkers as evidenced by a comparable expression detected in the random sections, which further qualifies them as highly expressed and robust potential biomarkers for ccRCC.

We next investigated the expression of the four biomarkers by IHC (Fig. 4A). Although transitional regions illustrated that the expression of PLOD2 was more intense in tumor sections (Supplementary Fig. S1D), a similar prevalence (score 3, >75%) was detected in both tumor and NATs (Fig. 4A). In contrast, SIRPα, FERMT3 and SPARC had a low prevalence score of 1 (5%–25% spreading) and a low intensity score of 1 in the normal tissues, while the expression was largely distributed and markedly elevated in the tumor tissues with a prevalence score between 2 (25%–75%) and 3 (>75%), and an intensity score of 3. Furthermore, for SIRPα, we determined a membranous localization in the tumor cells. Interestingly, for FERMT3 we observed almost no staining in large areas of the tumor tissues (Supplementary Fig. S1E), but high prevalence and intense staining in inflammatory regions (Fig. 4A), coinciding with the previously reported role of FERMT3 in leukocyte transmigration, platelet degranulation, and inflammation (37, 38). We were intrigued by an intense staining (score 3) for SPARC in the intratumoral stroma, which suggested that SPARC might be a secreted protein and could possibly be detected in body fluids of patients with ccRCC.

To this end, we assessed the secretion of the four suggested proteins into urine, which is a readily available biospecimen and can be collected noninvasively (Fig. 4B and C). We detected three of the four biomarkers in the urine samples of patients with ccRCC of different grades, as well as the nontumoral control group. While FERMT3 and SIRPα did not display a significant expression difference, SPARC was significantly more abundant in urine samples collected from the patients with cancer in comparison with the control group (P < 0.05; Fig. 4C), which supports our proteomics approach to find novel secreted biomarkers for ccRCC in urine.

Our comprehensive proteomic analysis and rigorous comparison between normal and tumor tissues provided us with new biomarkers for ccRCC. Next, we focused on the interpatient comparison in an attempt to identify a novel biomarker panel for the classification of patients based on their individual molecular expression profiles. To this end, we performed unsupervised hierarchical clustering to unveil tumor clusters. Indeed, the tumors grouped into two distinct clusters based on the expression profile of the 5,799 quantified proteins (Supplementary Fig. S4A). Furthermore, we challenged the robustness of the clustering by applying k-means clustering on the data, which recapitulated the same classification of the patients (Supplementary Fig. S4B). Increasing the number of the k-value (degree of freedom) from 2 to 4 did not remarkably change the patient clustering, which confirms that the observed two main clusters are stable. Next, we employed PCA, which showed concordance with the results of the hierarchical clustering and further confirmed the two main clusters (Supplementary Fig. S4C). Furthermore, statistical analysis revealed that a total of 599 proteins were differentially expressed between clusters (Student t test P < 0.05), with 354 significantly higher abundant in cluster 1 (cluster 1-specific) and 245 significantly higher abundant in cluster 2 (cluster 2-specific), respectively (Fig. 5A). Cluster 1-specific proteins were primarily involved in processes associated with the immune system, such as regulation of complement cascade and innate immune response (Fig. 5B), indicating that the observed enhanced immune cell infiltration in RCC (Fig. 2C and E) was primarily abundant in this subset of patients. Moreover, platelet activation and hemostasis implied neovascularization (Fig. 5B), and therefore a concomitant angiogenic profile of these tumors. Also noteworthy is the enrichment of components associated with the regulation of the insulin-like growth factor (IGF), which implied a primary dependence on this growth factor for tumor progression. In contrast, cluster 2 tumors displayed an enrichment of metabolic processes related to lipid, fatty acid, and amino acid degradation (Fig. 5C). Moreover, proteins associated with mitochondrial processes such as oxidative phosphorylation and TCA cycle were higher abundant in cluster 2 tumors.

To obtain a predictive biomarker panel sufficient for cluster discrimination, we further reduced the 599 cluster-significant proteins to a subset of representative proteins by applying cut-off values on the statistics (−log₁₀P > 2, t-test difference ± 2). The resulting 73 discriminative proteins, hereafter denoted as “marker genes,” comprised of 54 cluster 1-specific and 19 cluster 2-specific proteins, respectively (Fig. 5D). Besides their expected role in angiogenesis and immune response–related processes (Fig. 5B), marker genes of cluster 1 were also associated with TGFβ-related metastasis, apoptosis and HIF1α signaling (Fig. 5D). This further implies a progressive and malignant clinical behavior of cluster 1 tumors. Contrarily, besides their expected primary involvement in metabolic processes (Fig. 5C), marker genes of cluster 2 were also associated with mTOR signaling and adipocyte development (Fig. 5D), both of which have previously been described to promote RCC progression (2, 39).

For the validation of the proposed clustering profile, two independent larger datasets were used. First, the TCGA-KIRC RNASeq dataset (40), consisting of 539 tumors (normalized to 72 normal tissue counts), was utilized. In general, more than 5,500 mRNA-protein pairs displayed a significant positive correlation (0.38–0.43, P < 0.0001; Supplementary Fig. S5A) similar to previous observations in proteotranscriptomics analysis of ccRCC (7, 22), warranting the suitability of this transcriptomics data for validating our proteome-derived clustering signature. Furthermore, 565 cluster-significant proteins had an mRNA counterpart, which also demonstrated significant positive correlation (P < 0.0001), with a Spearman correlation of 0.6 and a Pearson correlation of 0.42, respectively (Fig. 6A). Surprisingly, proteome-transcriptome regulation discrepancy was higher for cluster 1-specific pairs (Spearman 0.49, Pearson 0.38) compared with cluster 2-specific pairs (Spearman 0.74, Pearson 0.69). This may suggest that metabolic processes in cluster 2 tumors are less prone to discordant regulation between the transcriptome and proteome level than immune- and angiogenesis-related processes in cluster 1 tumors.

Unsupervised hierarchical clustering of TCGA data based on the expression of these 565 mRNA/protein pairs classified the tumors into three clusters (TCGA-clusters), with two main clusters and 1 outlier cluster (n = 19), which was excluded from the validation (Fig. 6B). For cluster assignment, the significantly differentially expressed genes between TCGA clusters were compared with those from the proteome-clusters of this study and assigned on the basis of their overlap. TCGA cluster 1 was assigned with more than 84% concordance with proteome-cluster1, and TCGA-cluster 2 was assigned with more than 62% overlap with proteome-cluster2 (Supplementary Fig. S5B).

We further validated the clustering signature using the CPTAC data for ccRCC (22), which only included patients with tumor and NATs (n = 80). Of the determined cluster-significant proteins in our study, 547 were detected in the CPTAC proteome (Supplementary Fig. S5C), which displayed significantly high concordance in regulation (correlations >0.82, P < 0.0001). Unsupervised hierarchical clustering of the CPTAC data based on the expression of these common cluster-significant proteins also classified the tumors into two main clusters (CPTAC-clusters; Supplementary Fig. S6A). The subsequent cluster assignment resulted in an overlap of approximately 99% with cluster 1 proteins and 94% with cluster 2 proteins of this study, respectively (Supplementary Fig. S6B).

Next, we examined the prognostic difference between the clusters using the TCGA-KIRC follow-up data. Strikingly, Kaplan–Meier analysis of the overall survival (OS) revealed highly significant discrepancy in clinical outcome between the clusters (log-rank test P < 0.001; Fig. 6C; Supplementary Table S1). Patients of cluster 1 had an approximately 2.8 [95% confidence interval CI, (2.03–3.74)] times higher risk of death (HR) compared with patients of cluster 2 (Fig. 6C). The HR value was even higher for the recurrence-free survival (RFS) of this cohort with more than four times higher risk of death [95% CI, (2.59–6.71)] for patients of cluster 1 (Supplementary Fig. S6C; Supplementary Table S2). Furthermore, the 5-year survival rate of the patients differed noticeably between clusters (Fig. 6C). While approximately 43% [95% CI, (34.4–51.1)] of the patients in cluster 1 survived 5 years after diagnosis, more than 75% [95% CI, (69.6–80.8)] of the patients in cluster 2 were alive. This significant difference in survival between the patient clusters was also recapitulated by the 73 marker genes (Supplementary Fig. S6D; Fig. 5D).

The relatively higher mortality rate of cluster 1 patients was also supported by the evaluation of the clinicopathologic features (Supplementary Table S3). Remarkably, cluster 1 patients demonstrated significantly higher correlation with unfavorable clinical factors such as the presence of lymph node or distant metastasis, higher stages and grades, and recurrent disease (P < 0.005). Moreover, to find out whether the proposed patient clustering was among the confounding factors having a significant impact on the survival of patients with ccRCC, univariate and multivariate Cox regression analyses were performed with the follow-up data of the TCGA cohort (Supplementary Table S4; Table 1). Indeed, besides the age at diagnosis (P < 0.01), tumor stage (P < 0.005), tumor grade (P < 0.05) and the presence of distant metastasis (P < 0.0001; also the presence of lymph node metastasis in univariate analysis, P < 0.0001), the cluster membership was determined as a statistically significant predictor of survival in both analyses (P < 0.0001). In consideration of all parameters, the assignment to cluster 1 tumors indicated a similarly high risk of mortality as the presence of distant metastasis (HR of 2.11 and 2.18, respectively; Table 1).

To elucidate the cross-talk between marker genes of both clusters, we created a protein–protein interaction network complemented with Steiner nodes (Fig. 7A, markers from cluster 1 are labelled yellow, markers from cluster 2 are labelled green). We observed that the most prominent hubs CAV1, F2R, FLOT2, CFTR, APOA1, CANX, GJA1, LGALS3, ELANE, and FN1 were highly interconnected (Fig. 7A, nodes with thick border). Considering that two of the hubs and most of their interactors were cluster 1-specific, we hypothesized that this underlying interconnectivity might be an attractive target for drug treatment of cluster 1 tumors. This prompted us to annotate the network with FDA-approved drugs with inhibitor function. A total of 13 of the 87 cluster-discriminative marker genes and Steiner nodes were targets for 39 FDA-approved drugs (Fig. 7B; Supplementary Data S3). The majority of the pharmaceutical options that are used in RCC treatment target the VEGF(R) (sunitinib, pazopanib, sorafenib, axitinib) or mTOR (everolimus, temsirolimus) pathways (22). Although these annotated drugs for the marker genes are generally applied for diseases other than cancer such as hypertension, cardiovascular conditions, pulmonary embolism, anxiety disorders, osteoarthritis or inflammatory conditions, their repurposed application in combination with the FDA-approved drugs for ccRCC might be worthwhile as a future treatment perspective. For example, two marker gene hubs are targetable with drugs, namely, ELANE with alpha-1-proteinase inhibitor and CFTR with ibuprofen and dexibuprofen. Also, as shown in our study, immune system–related processes play a major role in the highly lethal cluster 1 tumors. Considering that five druggable marker genes (CA1, C5, MME, ELANE and ACAA1) were involved in the innate immune system response (Fig. 7A), stressing the importance of developing novel immunotherapy combinations besides the approved therapeutics nivolumab and ipilimumab against immune checkpoints is crucial for future treatment options.

In this study, we performed quantitative proteomics analysis of frozen ccRCC tissues and NATs and identified a total of 10,160 proteins, of which 5,799 were quantified by dimethylation (Fig. 2A). Moreover, 955 of these proteins were determined as statistically significant between tumor and NATs. A remarkable undertaking for the comprehensive proteogenomic characterization of ccRCC was performed by the CPTAC, representing the deepest proteome coverage of ccRCC with a total of 11,355 proteins identified and 7,026 proteins quantified from 103 patients (22). Comparing both proteome datasets revealed high conformity, with 4,540 of the quantified proteins of this study (84.5%) being shared with the CPTAC study (Supplementary Fig. S7A) and displaying significant positive correlation (P < 0.05) as indicated by a Spearman and Pearson correlation of 0.68 for each (Supplementary Fig. S7B). A total of 382 proteins of the significantly regulated proteins identified in the CPTAC study, were also discovered in this study, validating our data (Supplementary Fig. S7C). In addition, we identified 573 unique, differentially expressed proteins (Supplementary Data S2), promising a distinct proteome characterization of ccRCC compared with the CPTAC approach.

Out of the detected 409 significantly upregulated proteins, we selected four biomarkers, namely, PLOD2, FERMT3, SPARC and SIRPα that are potentially secreted, and confirmed their significantly elevated expression in ccRCC tumor tissues compared with NAT samples (Fig. 3). None of the biomarkers were significantly abundant in the proposed tumor clusters (Fig. 5A), which supports that these proteins are unaffected not only by intratumor heterogeneity (Fig. 3E), but also by interpatient tumor heterogeneity, supporting their eligibility as potential biomarkers of ccRCC. Although all four putative biomarkers are in the list of significantly upregulated proteins in other ccRCC proteomics studies (9, 22, 41, 42), they have not been validated yet by other approaches orthogonal to proteomics.

PLOD2 (lysyl hydroxylase 2, LH2) modifies telopeptidyl lysine residues of pro-collagen α chains on the endoplasmic reticulum and is generally regulated by HIF1α, TGFβ, and in RCC also by the tumor-suppressive miRNAs miR-26a-5p and miR-26b-5p (43–45). Aside from PLOD2, our PRM analysis also showed a higher abundance for the other known PLODs, PLOD1 (LH1) and PLOD3 (LH3) in ccRCC tumors (Fig. 3C and D). Despite its previously reported possible secretion in lung cancer (45), we could not detect PLOD2 in urine samples (Fig. 4B and C), however, it is worthwhile to investigate its suitability in the clinical setting by using a larger cohort and by measuring the urinary level of its collagen degradation products.

FERMT3, also known as Kindlin-3 and URP2, functions in hemostasis and thrombosis and plays a pivotal role in integrin-mediated cell-to-cell cross-talk as well as cell-matrix junctions (46, 47). Although its role in cancer remains unclear, it has been reported to be highly expressed in different types of lymphoma and breast cancer, leading to increased tumor growth, angiogenesis, and metastasis (46, 47). Interestingly, FERMT3 has been reported to act as a tumor suppressor in RCC, breast tumors, and melanoma based on mRNA data (46). The ambiguous role of FERMT3 in cancer warrants the need for further investigations.

SIRPα is a transmembrane protein which interacts with the ubiquitously expressed transmembrane protein CD47 (“SIRPα-CD47 axis”; ref. 48). This cross-talk contributes to the escape of tumor cells from the immune system response (48). Our data confirmed that SIRPα is located to the membrane in ccRCC tumor cells (Fig. 4A), however, we could also detect SIRPα in urine samples of patients with ccRCC as well as control patients, implicating that SIRPα is also a secreted protein. Given its role in immune evasion and its high abundance in ccRCC tissues, the SIRPα-CD47 axis can be considered as another promising target for immunotherapy besides nivolumab against PD-1 or ipilimumab against CTLA-4 (49).

SPARC is a multi-faceted matricellular glycoprotein that is involved in tissue remodeling, morphogenesis, and bone mineralization (50). While it functions as a tumor suppressor in neuroblastomas, colorectal and ovarian cancer (50, 51), it has been reported to be elevated on the mRNA level in sarcomatoid RCC, as well as most of ccRCC tumors (70%; ref. 52). SPARC has been reported to be expressed only by tumor-associated stromal cells such as vascular endothelial cells and fibroblasts, but not by cancer cells themselves (52). In line with this, we observed for SPARC an intense staining in the intratumor stroma (Fig. 4A). Furthermore, we detected a significantly elevated abundance of SPARC in urine samples of patients with ccRCC proving our concept to find secreted biomarkers from our comprehensive proteomics analysis (Fig. 4B and C). However, more advanced and elaborated approaches in a larger independent cohort are required to validate SPARC as robust screening biomarker for urine tests.

Here, we also propose a biomarker panel for the appropriate classification of ccRCC tumors based on the underlying molecular phenotypes, since classification solely based on histologic or morphologic features of tumors and subsequent treatment decisions can lead to varying treatment response (20, 53). In this study, we identified two ccRCC tumor clusters that differed in their molecular features. Our clustering signature is distinct from previously suggested ones such as the widely investigated transcriptomics based ClearCode34 grouping system, which has only five proteins in common (Supplementary Fig. S7D; ref. 21). Both the ClearCode34 as well as our classification system showed similar HRs for the OS between the patient clusters (HR 2.4 and 2.2, respectively). However, regarding the RFS, our clustering signature showed a 4.2-fold difference in the risk of death between the clusters while the difference was 2.3-fold between the ClearCode34 clusters. Furthermore, our clustering signature distinguishes from the proteomics-based classification system ccRCC1–3 (22), which showed malignant features intermixed across the three identified clusters, while in this study we showed that highly malignant processes such as EMT, HIF1α signaling, innate immune response, and angiogenesis were condensed in cluster 1 tumors accounting for the relatively worse survival outcome (Fig. 5D; Fig. 6C), while cluster 2 tumors had a low-risk profile and a more favorable prognosis.

Advanced RCC is currently treated with monotherapies or combination therapies of immune checkpoint inhibitors and VEGF or mTOR inhibitors in first- and second-line strategies (22, 54). Hypothesizing from our data, the highly immunogenic and angiogenic cluster 1 tumors might be receptive to immunotherapy in combination with VEGF inhibitors (Supplementary Data S3). In recent years, first-line approaches combining immunotherapeutics such as avelumab or pembrolizumab with the VEGF inhibitor axitinib have proven to be promising alternatives for the treatment of progressed RCC (54). However, using immunotherapies as first-line treatment excludes them as second-line therapy due to poor response. Thus, monotherapies with cabozantinib or axitinib are the preferred choices for second-line treatment (54). In contrast to cluster 1 tumors, mTOR signaling was elevated in cluster 2 tumors, which might be targetable by the RCC-approved mTOR inhibitor everolimus (Supplementary Data S3).

Furthermore, we identified 13 cluster-discriminative marker genes targetable by repurposed FDA-approved drugs (Fig. 7) that can provide alternative treatment strategies to meet the need for novel anticancer therapeutics (55). One highly elevated EMT marker in cluster 1 tumors is the signaling protein S100A4, which is a target of the antipsychotic drug trifluoperazine (Fig. 7A). It has been shown that trifluoperazine can be repurposed as an anticancer drug and successfully prevent or minimize metastatic spread in different cancers by inhibiting S100A4 through protein oligomerization (56, 57). Furthermore, trifluoperazine has been shown together with doxorubicin to decrease multi-drug resistance in RCC cell lines (58). Targeting S100A4 could also be advantageous for reducing hypoxic signaling in cluster 1 tumors, because it was revealed as target of HIF1α (Supplementary Data S3). Another gene highly overexpressed in cluster 1 tumors upon hypoxic signaling is SLC2A1 (Supplementary Data S3), also known as glucose transporter 1 (GLUT1). This marker gene can be targeted by the intravenous anesthetic etomidate (Fig. 7A), which has been reported to reduce GLUT1-mediated glucose transport and associated tumor growth in lung cancer (59, 60). SLC2A1 has also been revealed to be related to angiogenic processes such as coagulation and platelet degranulation, as well as the marker CA1 (Supplementary Data S3), which participates in hypoxia-induced glucose transport, angiogenesis, and pH regulation (61). Multiple inhibitory FDA-approved sulfonamides such as acetazolamide, methazolamide, dichlorophenamide, dorzolamide, and brinzolamide can target CA1 (Fig. 7A), and have been reported to counteract metastatic dissemination in different RCC cell lines (62).

Taken together, our comprehensive proteomics analysis identified four novel biomarkers for the detection of ccRCC. Future directions should entail the investigation of the role of these biomarkers in ccRCC tumorigenesis and metastatic spread, as well as the application of SPARC as an analytical tool in clinical assays. We further uncovered two distinct molecular phenotypes of ccRCC tumors and identified a cluster-discriminative biomarker panel targetable by several repurposed FDA-approved drugs. Future directions warrant the investigation of tumor-oriented intervention by combining these novel therapeutics with current medical treatment options for RCC.

No disclosures were reported.

A. Senturk: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft. A.T. Sahin: Formal analysis, validation, investigation, visualization, methodology. A. Armutlu: Supervision, investigation, visualization. M.C. Kiremit: Resources, supervision, investigation. O. Acar: Resources, supervision, writing–review and editing. S. Erdem: Resources, supervision, writing–review and editing. S. Bagbudar: Resources, formal analysis. T. Esen: Supervision, writing–review and editing. N. Tuncbag: Conceptualization, supervision, project administration, writing–review and editing. N. Ozlu: Conceptualization, supervision, funding acquisition, project administration, writing–review and editing.

We gratefully acknowledge the KUPAM and KUTTAM OMICS facilities of Koc University, especially Busra Akarlar for technical support. We thank Ozgur Kurt and Gulsum Citak from Koc University School of Medicine for their essential help in pathologic assessment of tissue samples. We further thank Serkan Kir for valuable scientific comments on the results. The results shown here are in whole or part based upon data generated by TCGA Research Network: https://www.cancer.gov/tcga. Data used in this publication were generated by the Clinical Proteomic Tumor Analysis Consortium (NCI/NIH). The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE (63) partner repository with the dataset identifier PXD024696.