indicates AE + patient withdrawal.

All adjuvant TKI trials were placebo-controlled and aimed to show the DFS benefit, but only S-TRAC yielded a positive result, with a 1.2-year improvement in the DFS of the sunitinib arm. Tumor cell PD-L1 expression was not statistically significantly associated with DFS, whereas high tumor CD8+ T-cell density was predictive for longer DFS in the sunitinib arm of the S-TRAC trial [105]. The S-TRAC, PROTECT, and ATLAS trials included only ccRCC patients, and the majority (79% and 84%) of patients enrolled in the ASSURE and SORCE trials had ccRCC. Usually, the protocol-specified duration of adjuvant TKI therapy was 12 months. The ATLAS and the SORCE trials included cohorts with adjuvant TKI therapy up to 36 months, but the longer duration of TKI therapy did not lead to improved DFS. Adjuvant TKI therapy caused substantial toxicity (grade 3–4 adverse events 49–72%) and a significant proportion of the patients (23–49%) discontinued adjuvant TKI therapy because of intolerable toxicity or refused to continue study therapy (96–100). Currently, adjuvant TKI therapy is not recommended after complete resection of the primary tumor in the international RCC guidelines due to the substantial toxicity and the lack of OS benefit [4,94].

Immune checkpoint inhibitors (ICI) have replaced cytokines in the immune therapy of advanced RCC and are also being studied in randomized placebo-controlled prospective clinical trials in the adjuvant and neoadjuvant setting. IMmotion010 is evaluating 12-month adjuvant therapy with PD-L1 inhibitor atezolizumab, PROSPER neoadjuvant therapy (nivolumab two doses), followed by 9-month adjuvant therapy with PD-1 inhibitor nivolumab and CheckMate 914 6-month adjuvant therapy with the combination of CTLA-4 inhibitor ipilimumab and PD-1 inhibitor nivolumab in resected localized ccRCC patients, and RAMPART 12-month durvalumab adjuvant therapy and 12-month adjuvant CTLA-4 and PD-L1 inhibitor (tremelimumab and durvalumab) combination therapy. The first results of these trials are expected to be published in 2022–2024. The results from the KEYNOTE-564 trial evaluating 12-month adjuvant therapy with pembrolizumab in resected intermediate- or high-risk ccRCC patients showed a statistically significantly longer recurrence-free survival rate in the pembrolizumab arm compared to the placebo arm at 24 months (77.3% vs. 68.1%, HR for recurrence or death 0.68 (0.53–0.87)) (Table 2) [105]. As this was the first analysis, a longer follow-up will be needed to confirm the survival outcomes of the pembrolizumab adjuvant therapy. However, ICI may finally become a practice-changing adjuvant treatment option for RCC patients after complete resection of the primary tumor and lymph node or distant metastases.

#### **6. Discussion**

Numerous traditional histopathological factors and an increasing number of biomarkers have been identified to affect the postoperative prognosis of patients with localized ccRCC. The individual assessment of the risk for disease recurrence after radical or partial nephrectomy is important to tailor the intensity of postoperative follow-up imaging. Moreover, accurate risk assessment for disease recurrence is essential to select optimal patients for adjuvant therapy trials. However, there is no consensus regarding which is the best model or biomarker to choose to guide the clinical decision making. Limitations in the availability of biomarker analyses, time required to obtain the results, costs from the analyses, and, in particular, the lack of sufficient clinical validation still limit the use of prognostic biomarkers in clinical practice. Useful risk assessment tools for clinicians should be easy-to-use and include only a moderate amount of readily available risk factors (e.g., 3–5 traditional histopathological factors). Different clinicopathological features may be available in different centers. In the future, biomarkers, including those from plasma and urine (liquid biopsies), may supplement these prognostic algorithms.

Prognostic models with traditional histopathological and clinical factors should be easy-to-use and readily available. However, only the SORCE trial had incorporated a prognostic algorithm into the inclusion criteria of the trial. Adjuvant TKI trials underscored the fact that careful patient selection is required to avoid substantial toxicity and enrich higher-risk patients for adjuvant therapy. A meta-analysis of adjuvant TKI trials showed a DFS benefit in the high-risk (T3, Fuhrman grade 3–4; T4, or N1) population but not in the low-risk population (pooled HR for DFS 0.85 (0.75–0.97) and 0.98 (0.82–1.17), respectively) [106]. The optimal selection criteria for the high-risk localized ccRCC population remain to be defined. The results from the biomarker analyses of current neoadjuvant and adjuvant trials with ICI may shed more light on the issue.

#### **7. Conclusions**

Prognostic factors and validated prediction models help to evaluate the risk for disease recurrence after complete surgical resection of localized ccRCC. Better models to reduce follow-up imaging in low-risk patients and optimize the selection of patients for adjuvant trials are required. The combination of clinical and histopathological features with novel biomarkers may improve the prediction accuracy of prognostic models. The optimal set of prognostic factors and biomarkers to define high-risk patients for disease recurrence may be discovered in ongoing placebo-controlled randomized prospective clinical trials.

**Author Contributions:** Conceptualization, K.E.M., P.V. and P.M.J.; writing—original draft preparation, K.E.M., P.V. and P.M.J.; writing—review and editing, K.E.M., P.V. and P.M.J.; supervision, P.M.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research work was supported by the Finnish Cancer Unions and Turku University Hospital (project 13031).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **The Hypoxic Microenvironment Induces Stearoyl-CoA Desaturase-1 Overexpression and Lipidomic Profile Changes in Clear Cell Renal Cell Carcinoma**

**Juan Pablo Melana 1, Francesco Mignolli 2, Tania Stoyanoff 1, María V. Aguirre 1, María A. Balboa 3,4, Jesús Balsinde 3,4,\*,† and Juan Pablo Rodríguez 1,\*,†**


**Simple Summary:** Clear cell renal cell carcinoma (ccRCC) is characterized by a high rate of cell proliferation and an extensive accumulation of lipids. Uncontrolled cell growth usually generates areas of intratumoral hypoxia that define the tumor phenotype. In this work, we show that, under these microenvironmental conditions, stearoyl-CoA desaturase-1 is overexpressed. This enzyme induces changes in the cellular lipidomic profile, increasing the oleic acid levels, a metabolite that is essential for cell proliferation. This work supports the idea of considering stearoyl-CoA desaturase-1 as an exploitable therapeutic target in ccRCC.

**Abstract:** Clear cell renal cell carcinoma (ccRCC) is the most common histological subtype of renal cell carcinoma (RCC). It is characterized by a high cell proliferation and the ability to store lipids. Previous studies have demonstrated the overexpression of enzymes associated with lipid metabolism, including stearoyl-CoA desaturase-1 (SCD-1), which increases the concentration of unsaturated fatty acids in tumor cells. In this work, we studied the expression of SCD-1 in primary ccRCC tumors, as well as in cell lines, to determine its influence on the tumor lipid composition and its role in cell proliferation. The lipidomic analyses of patient tumors showed that oleic acid (18:1*n*-9) is one of the major fatty acids, and it is particularly abundant in the neutral lipid fraction of the tumor core. Using a ccRCC cell line model and in vitro-generated chemical hypoxia, we show that SCD-1 is highly upregulated (up to 200-fold), and this causes an increase in the cellular level of 18:1*n*-9, which, in turn, accumulates in the neutral lipid fraction. The pharmacological inhibition of SCD-1 blocks 18:1*n*-9 synthesis and compromises the proliferation. The addition of exogenous 18:1*n*-9 to the cells reverses the effects of SCD-1 inhibition on cell proliferation. These data reinforce the role of SCD-1 as a possible therapeutic target.

**Keywords:** kidney; hypoxia; tumor microenvironment; SCD-1; oleic acid

#### **1. Introduction**

Renal cell carcinoma (RCC) is the most common malignancy of the urinary system. Although the incidence of RCC has remained stable, the mortality rates have decreased by only 0.9% each year from 2007 to 2016 [1,2].

**Citation:** Melana, J.P.; Mignolli, F.; Stoyanoff, T.; Aguirre, M.V.; Balboa, M.A.; Balsinde, J.; Rodríguez, J.P. The Hypoxic Microenvironment Induces Stearoyl-CoA Desaturase-1 Overexpression and Lipidomic Profile Changes in Clear Cell Renal Cell Carcinoma. *Cancers* **2021**, *13*, 2962. https://doi.org/10.3390/ cancers13122962

Academic Editors: Claudia Manini and José I. López

Received: 12 April 2021 Accepted: 10 June 2021 Published: 13 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Clear cell renal cell carcinoma (ccRCC) represents the most common subtype (>80%) of RCC. The most striking phenotypic feature of ccRCC is its clear cell morphology, which has been linked to a high lipid and glycogen accumulation [3]. Neutral lipids such as triacylglycerol (TAG) and cholesterol esters (CE) are stored in prominent cytoplasmic lipid droplets (LD), which are critical for cell growth and maintenance of the cell membrane [4]. Although the presence of these droplets in ccRCC is critical for sustained tumorigenesis, their contribution to lipid homeostasis and tumor cell viability is not completely understood [5].

A ubiquitous metabolic event in cancer is the constitutive activation of the pathway for fatty acid biosynthesis. Saturated fatty acids (SFAs), monounsaturated fatty acids (MUFAs) and polyunsaturated fatty acids (PUFAs) are synthetized to sustain the growing demand for phospholipids (PLs) that are used for the assembly of new membranes, energy storage and cell signaling [6–8]. The activation of enzymes such as ATP-citrate lyase (ACL), acetyl-CoA carboxylase (ACC) and fatty acid synthase (FAS) leading to an increased synthesis of SFAs has been extensively studied [9–12]. SFAs later become MUFAs by the action of stearoyl-CoA desaturase-1 (SCD-1) [13]. SCD-1 is a Δ9-fatty acyl-CoA desaturase that catalyzes the insertion of a double bond in the cis-Δ9 position of several saturated fatty acyl-CoAs—mainly, palmitoyl-CoA and stearoyl-CoA—to produce palmitoleoyland oleoyl-CoA, respectively [14]. It has been reported that these unsaturated fatty acids affect several crucial biological functions of tumor cells, such as proliferation, signaling, invasiveness and apoptosis. It was shown that oleic acid (18:1*n*-9), one of the most prevalent free fatty acids (FFAs) in human plasma, increases the proliferation of human prostate, breast and renal cancer cells [15,16]. As noted above, it was suggested that SCD-1 could be a therapeutic target in oncology, since its pharmacological inhibition induces tumor cell apoptosis [14,17–21]. Controversially, it is known that certain fatty acids such as 18:1*n*-9 exert anticancer effects on many tumors, inhibiting cell proliferation and favoring apoptosis [22,23].

Solid tumors such as ccRCC often show hypoxic areas as a result of uncontrolled tumor growth, without a proper development of its associated vascular network [13]. Hypoxia inducible factors (HIF-1α and HIF-2α) are commonly stabilized key players connected to cell growth and metabolic reprogramming in ccRCC. Both factors modulate tumoral hypoxic responses through altering the cell energy metabolism, including the modification of glucose consumption [24] and the expression of a lipid metabolism-associated gene [13,25,26].

We have previously shown that SCD-1 expression correlates with other cellular markers of the tumor hypoxic microenvironment, such as EPO, EPO-R, VEGF and VEGF-R in ccRCC [27]. However, up to now, tumor hypoxia in ccRCC has not been associated with the induction of SCD-1 and its consequent modification of the tumor lipidomic profile.

In this work, we demonstrate that cellular hypoxia favors the induction of SCD-1, and this influences the cellular lipid phenotype. Furthermore, we observed that SCD-1 inhibition deprives cells of essential lipid metabolites for cell proliferation. These data provide evidence to consider SCD-1 as an exploitable therapeutic target in ccRCC.

#### **2. Materials and Methods**

#### *2.1. Patients and Sampling Procedures*

Samples of ccRCC (*n* = 12) were obtained from patients treated for radical nephrectomy in the Urology Unit of the Hospital Dr. José Ramón Vidal (Corrientes, Argentina) between 2015 and 2018. The normal distal tissues and ccRCC of the same affected kidney were surgically removed. Samples were aseptically transported to the laboratory and quickly processed. They were then fixed for histopathology and immunohistochemistry procedures.

The design and methods of this research were approved by the Bioethics Committee of the Medical Research Department at Dr. José Ramón Vidal Hospital in Corrientes, Argentina. Written informed consent was obtained from each donor. The researchers received the samples in an anonymous manner.

#### *2.2. Cell lines, Proliferation and Viability Assays*

Caki-1 and Caki-2 cell lines (derived from ccRCC), originally from the American Type Cell Culture Collection, were generously provided by Dr. Alfredo Martínez Ramírez (Centro de Investigaciones Biomédicas de La Rioja, Logroño, Spain) and Dr. Ricardo Sánchez Prieto (Universidad de Castilla-La Mancha, Albacete, Spain), respectively. The cells were cultured in McCoy's 5a medium modified (Thermo Fisher, Madrid, Spain) with 10% fetal bovine serum and 50 μg/mL of gentamicin (Invitrogen, Carlsbad, CA, USA) at 37 ◦C under humidified conditions with 5% CO2. Cell proliferation and viability were measured with a Neubauer chamber and also using the CellTiter 96® AQueous One Solution Cell Proliferation Assay Kit (Promega, Madison, WI, USA).

For in vitro experiments, cells were subcultured every 3 to 4 days after reaching 80–90% confluence. The cells were trypsinized, centrifuged and resuspended in the medium at a suitable density. Experiments utilizing exogenous 18:1*n*-9 were performed under serum-restrictive conditions (1%) [28].

#### *2.3. Hypoxic Microenvironment*

To achieve a hypoxic microenvironment similar to the tumor and the effects of HIF stabilization, Caki-2 cells were exposed to different nontoxic concentrations of CoCl2 in McCoy's 5a medium modified with 1% fetal bovine serum [29–35]. It has been previously shown that CoCl2 inhibits the hydroxylation of HIF-1α, thus stabilizing HIF-1α and achieving the desired hypoxic effect [36–39].

#### *2.4. Real-Time Quantitative PCR (RT-qPCR)*

*SCD-1*, *HIF-1A* and *HIF-2A* mRNA were determined by RT-qPCR. Total RNA was extracted using the TRIzol reagent method (Invitrogen) according to the manufacturer's protocols. The obtained total RNA was purified using Ambion® TURBO DNA-free™. First-strand cDNA was obtained by using the Moloney murine leukemia virus reverse transcriptase (Promega) from 1 μg of RNA. qPCRs were then performed using specific primers for *SCD-1* as follows: 5 -TTCCTACCTGCAAGTTCTACACC-3 (forward) and 5 -CCGAGCTTTGTAAGAGCGGT-3 (reverse) with a product of 116 bp. *HIF-1A*: 5 - TGCTGGGGCAATCAATGGAT-3 (forward) and 5 -CTACCACGTACTGCTGGCAA-3 (reverse) with a product of 590 bp. *HIF-2A*: 5 -TATAGTGACCCCGTCCACGT-3 (forward) and 5´-AGGGCAACACACACAGGAAA-3 (reverse) with a product of 572 bp. *B-ACTIN*: 5 -CATGTACGTTGCTATCCAGGC-3 (forward) and 5 -CTCCTTAATGTCACGCACGAT-3 (reverse) with a product of 250 bp was used as the housekeeping gene.

All primers were tested for specificity using the primer BLAST program available at the National Center for Biotechnology Information website (www.ncbi.nlm.nih.gov; accessed on 1 February 2020). Cycling conditions were: 1 cycle at 95 ◦C for 12 min, 40 cycles at 95 ◦C for 15 s, 60 ◦C for 20 s, 72 ◦C for 20 s and a final extension at 72 ◦C for 10 min.

#### *2.5. SCD-1 Inhibition Assays*

CAY 10566, a potent selective SCD-1 inhibitor, was purchased from Cayman Chemical (Ann Arbor, MI, USA), dissolved in DMSO and used (3 μM) according to the manufacturer's recommendations (noncytotoxic concentrations). Caki-2 cells were cultured in McCoy's 5a medium modified with 1% fetal bovine serum.

#### *2.6. Analyses of Fatty Acids by Gas Chromatography Coupled to Mass Spectrometry (GC/MS)*

Cellular lipids were extracted using the method of Bligh and Dyer [40]. After the addition of the appropriate standards, lipids were separated by thin-layer chromatography (TLC) using as the stationary phase silica gel 60 and a mobile phase consisting of

*n*-hexane/ethyl ether/acetic acid (70:30:1 *v*/*v*/*v*) [41]. Glycerolipids and glycerophospholipids were transmethylated with 500 μL of 0.5-M KOH in methanol for 30 min at 37 ◦C, and 500 μL of 0.5-M HCl was added to neutralize. For the transmethylation of cholesterol esters, the samples were resuspended in 400-μL methyl propionate and 600 μL of 0.84-M KOH in methanol for 1 h at 37 ◦C. Afterward, 50 mL and 1 mL of acetic acid and water, respectively, were added to neutralize. Analysis of the fatty acid methyl esters was carried out using an Agilent 7890A gas chromatograph coupled to an Agilent 5975C mass-selective detector operated in the electron impact mode (EI, 70 eV) equipped with an Agilent 7693 autosampler and an Agilent DB23 column (60-m length × 250-μm internal diameter × 0.15-μm film thickness) under the conditions described previously [42–45]. Data analysis was carried out with Agilent G1701EA MSD Productivity Chemstation software, revision E.02.00.

#### *2.7. Confocal Microscopy*

Caki-2 cells attached to coverslips were incubated for 24 h in McCoy 5a medium modified with different concentrations of CoCl2 and a positive control with 30-μM oleic acid. Cells were then washed with phosphate-buffered saline and incubated with BODIPY 493/503 staining solution (2 μg/mL) for 15 min at 37 ◦C. Cells were subsequently washed, fixed with 4% paraformaldehyde and washed again. The coverslips were mounted on slides with a DAPI reagent (1 μg/mL). Untreated cells were used as negative controls. Fluorescence was monitored by microscopy using a Bio-Rad confocal system Radiance 2100 laser scanner (Bio-Rad, Richmond, VA, USA). The images were analyzed with ImageJ software.

#### *2.8. Apoptosis Detection by Flow Cytometry*

The effect of SCD-1 inhibition on apoptosis was evaluated by flow cytometry. Based on the preliminary time–course data, the exposure time was set to 18 h, and apoptosis was analyzed by labeling with the annexin V-fluorescein isothiocyanate (FITC) apoptosis detection kit (BD Bioscience, San Jose, CA, USA), which recognizes phosphatidylserine exposure on the outer leaflet of the plasma membrane. After washing the cells, cell fluorescence was quantified by flow cytometry in FL1 (Gallios; Beckman Coulter, Barcelona, Spain). Data were analyzed with FlowJo software version 8.7. The propidium iodide (PI; Sigma-Aldrich, Madrid, Spain) uptake was analyzed by incubating cells with 50-μg/mL PI in PBS in the dark for 5 min. Fluorescence was quantified by flow cytometry in FL3. Data were analyzed with FlowJo version 8.7.

#### *2.9. Statistics*

Statistics were performed using GraphPad Prism 8.0 via an unpaired *t*-test or oneway analysis of variance (ANOVA), followed by Bonferroni's or Tukey's comparison tests. Differences were considered to be significant at *p* < 0.05.

#### **3. Results**

#### *3.1. The Lipidomic Profile of ccRCC Is Dependent on the Tumor Area Analyzed*

ccRCC tumors frequently show visible macroscopic differences with defined boundaries between the center and external areas. Therefore, we first performed a lipidomic analysis of fatty acids by GC/MS of two arbitrarily separated tumor sections: the core and periphery. Figure 1 shows that the fatty acid profile of cellular PLs did not show marked differences when the control was compared with the different tumor sections: core or periphery. Both had similar amounts and types of fatty acids, with the exception of oleic acid (18:1*n*-9) and arachidonic acid (20:4*n*-6), which were increased in the core. In contrast, the distribution of fatty acids in the TAG and CE fractions showed a higher amount of lipids in the core (Figure 1B,C). Palmitic acid (16:0), stearic acid (18:0) and, particularly, 18:1*n*-9 were greatly increased in the core compared to normal tissue or periphery.

**Figure 1.** Lipidomic profile of ccRCC. (**A**) The profile of major phospholipid fatty acids in healthy distal normal tissue (blue bars) or tumors (core and periphery: red and green, respectively) were determined by GC/MS after converting the fatty acid glyceryl esters into fatty acid methyl esters. (**B**,**C**) Profile of fatty acids present in neutral lipids (TAG and CE). Data are expressed as the means ± SEM (*n* = 12). \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001 and \*\*\*\* *p* < 0.0001, significantly different from the control.

Since interindividual genotypic variations create great variability in primary cell culture models derived from tumors [46], in the following series of experiments, we used the Caki-1 and Caki-2 cell lines as an in vitro model of ccRCC. To compare the lipidomic profile of the Caki-1/-2 cell lines, we first analyzed the total cellular fatty acid content. Similar to that observed in tumors, Figure 2 shows that the most abundant fatty acids in these cell lines were also 16:0, 18:0, 18:1*n*-9 and 20:4*n*-6. Although the two cell lines showed a similar fatty acid distribution, Caki-2 had a slightly higher amount of 18:1*n*-9. While both cell lines are ccRCC models, Caki-2 was established from a primary clear cell carcinoma of the kidney, and Caki-1 was derived from a skin metastasis. Consequently, we decided to use the Caki-2 cell line for the in vitro experiments.

**Figure 2.** Lipidomic fatty acids profile of the Caki-1 and Caki-2 cell lines. (**A**) Total cellular fatty acid content. (**B**) Distribution of SFAs, MUFAs and PUFAs. Data are expressed as the means ± SEM and are representative of three independent experiments.

#### *3.2. The Hypoxic Microenvironment Promotes SCD-1 Overexpression, Lipid Droplet Formation and Changes in the Cellular Fatty Acids Profile*

Differences in the lipid composition, depending on the area of the tumor analyzed, should be in line with the expression pattern of the enzymes involved in their cellular metabolic pathways. Likewise, enzyme induction is strictly linked to the tumor microenvironment. Analyzing renal tumors, we previously showed that there is a statistical association of some hypoxia markers (e.g., *HIF-1A*) with the expression of *SCD-1* [27]. Thus, to investigate whether this physiological condition is actually responsible, at least in part, for SCD-1 induction, the Caki-2 cells were treated with different CoCl2 concentrations for 24 h to generate chemical hypoxia in vitro [30,39].

We first determined the cytotoxicity after the treatment with CoCl2 for 24 h and cell proliferation rates (growth constants (k) and cell doubling times; Supplementary Materials Figure S1). We observed that concentrations below 300 μM did not affect the cellular viability [47], but higher concentrations, up to 400 μM, induced cell death (5–10%). To verify that this salt did indeed generate cellular hypoxia, we evaluated by RT-qPCR the expression of hypoxia markers such as *HIF-1A* and *HIF-2A* (Figure 3A). Under the same conditions, the expression of *SCD-1* mRNA significantly increased (up to 200-fold) with 300-μM CoCl2 (Figure 3B). We hypothesized that this increase in SCD-1 would be expected to lead to elevated amounts of intracellular 18:1*n*-9, which, in turn, would induce modifications in the lipid profile and/or cytoplasmic morphological changes. Hence we next evaluated the Caki-2 cells treated with 30-μM 18:1*n*-9 as a positive control for LD formation [48] with the cells treated with CoCl2 (0–300 μM) by confocal microscopy. Using BODIPY® staining, we detected increased LD biogenesis that was dependent on the hypoxic conditions (Figure 3C).

**Figure 3.** Chemical hypoxia promotes SCD-1 expression in vitro. (**A**) Caki-2 cells were exposed to different concentrations of CoCl2 (0–400 μM) for 24 h, and the development of the hypoxic microenvironment was tested with the expression of *HIF-1A* and *HIF-2A* by RT-qPCR. (**B**) Then, *SCD-1* overexpression was detected under the same experimental conditions by RT-qPCR. (**C**) The evaluation of LD formation was determined by confocal microscopy. Caki-2 cells were exposed to 18:1*n*-9 (30 μM) for 24 h as a positive control for LD formation. Magnification 400×. Data are expressed as the means ± SEM and are representative of three independent experiments. \*\* *p* < 0.01, \*\*\* *p* < 0.001 and \*\*\*\* *p* < 0.0001, significantly different from the control.

In line with these phenotypic modifications, we evaluated the hypoxia-induced lipidomic changes in Caki-2 cells using GC/MS. Notably, the saturated fatty acids (16:0 and 18:0) experienced a significant reduction (*p* < 0.001) in all treatments in the PL fraction (Figure 4A). Consistent with the induction of *SCD-1* and the appearance of cytoplasmic LD, we observed a significant increase in 18:1*n*-9 but only in the TAG and CE fractions under hypoxic conditions (*p* < 0.01) (Figure 4B,C).

**Figure 4.** Lipidomic profile of the Caki-2 cells under hypoxic conditions. The cells were treated with the indicated concentrations of CoCl2, and (**A**) the fatty acid contents in the phospholipids (PL), (**B**) triacylglycerol (TAG) and (**C**) cholesterol esters (CE) were determined by GC/MS. Data are expressed as the means ± SEM and are representative of three independent experiments. \*\* *p* < 0.01 and \*\*\*\* *p* < 0.0001, significantly different from the control.

#### *3.3. Oleic Acid Is Essential for ccRCC Cell Proliferation*

As additional evidence for the role of SCD-1 in the biosynthesis of MUFAs, the pharmacological inhibition of the enzyme was performed in Caki-2 cells using the selective inhibitor CAY 10566 at a nontoxic concentration (3 μM) for 24 h [49]. As shown in Figure 5A, marked lipid changes in the 18:0 and 18:1*n*-9 levels were detected. The total cellular fatty acid profile, considering all the lipid fractions simultaneously, showed a significant increase of 18:0, with a consequent decrease of 18:1*n*-9 (*p* <0.01).

In line with that demonstrated by other authors with different methodologies [50], we noted that prolonged exposure times (longer than 24 h) to the enzyme inhibitor (3 μM) induced drastic decreases in the cell viability (34.65% ± 2.97; *p*< 0.001).

**Figure 5.** SCD-1 pharmacological inhibition. (**A**) Caki-2 cells were treated with the enzyme inhibitor CAY 10566 (3 μM) for 24 h. (**B**) Experimental design used in SCD-1 inhibition (slanted arrows indicate the addition of 18:1*n*-9 at different times). (**C**) Cell viability was measured with the CellTiter 96® kit. In all cases, the cells were cultured with 1% fetal bovine serum. Data are expressed as the means ± SEM and are representative of three independent experiments. \* *p* < 0.05, \*\* *p* < 0.01 and \*\*\* *p* < 0.001, significantly different from the control.

If the reduced cellular levels of 18:1*n*-9 play a role in arresting the cell growth or triggering apoptosis, its addition to the cell culture would restore, at least in part, the cell proliferation. In order to check this hypothesis, 18:1*n*-9 (50 μM) was added at different times (Figure 5B,C) to Caki-2 cells with and without a treatment with CAY 10566 (3 μM). We observed that the addition of 18:1*n*-9, along with CAY 10566 or two hours after, preserved the cell viability; this effect was not observed if the fatty acid was added at later time points (Figure 5C).

In addition to evaluating the cell viability with CellTiter, we determined whether the inhibition of SCD-1 with CAY 10566 induced apoptosis in ccRCC. The Caki-2 cells were stained with annexin V-FITC and propidium iodide (50 μg/mL), as detailed in the Materials and Methods. Figure 6 shows that the cells treated with CAY 10566 manifested a small increase in apoptotic cell death, which was fully preventable if 18:1*n*-9 was present in the incubation media. Collectively, these data suggest that the decrease in cell viability that the cells experienced when SCD-1 was inhibited by CAY 10566 (Figure 5) was due to a reduced proliferation rate rather than drug-induced apoptotic cell death.

**Figure 6.** Analysis of the apoptotic markers in Caki-2 cells treated with a SCD-1 inhibitor and 18:1*n*-9. Apoptotic changes in the plasma membrane were detected by simultaneous staining with annexin V-fluorescein isothiocyanate (FITC) (FL1) and propidium iodide (PI; FL3) (**A**) Gating strategy in the control and treated cells. (**B**,**C**) Cells were untreated (control) or treated with 3-μM CAY 10566, respectively. (**D**–**F**) Afterward, 18:1*n*-9 (50 μM) was added to the cells at 2, 6 and 18 h, respectively.

#### **4. Discussion**

ccRCC tumors characteristically show a bright yellow color in the center as a result of its abundant lipid content and a variegated appearance in the boundary with hemorrhage, necrosis and/or fibrosis with a frequently well-circumscribed capsule or pseudo-capsule that separates the tumor from the adjacent tissues [51]. In order to address these differences in tumor heterogeneity, we developed a separate core and periphery analysis compared to the normal distal tissues, similarly to other ccRCC-focused studies [52,53]. In all the samples individually analyzed, the highest lipid content in the core was a constant pattern, despite the interindividual differences usually determined in the oncological analysis [54–56]. Consistent with these results, Saito et al., developing a broader study of untargeted lipidomics, determined that these tumors have large accumulations of CE and TAG among the other lipids [53].

These imbalances are obviously associated with the particular conditions that tumor cells show to adapt their metabolism to an uncontrolled and continuous growth [57]. Several enzymes within the fatty acid biosynthesis pathway have been found to be essential for cancer cell growth or survival and are currently tracked as possible targets for therapeutic development [13,58]. Among them, SCD-1 was demonstrated to be a key regulator of the MUFA/SFA balance in several cancer cells, and its blockade triggers apoptosis [8,59,60]. In particular, SCD-1 was shown to be overexpressed in ccRCC [1] and, therefore, proposed as a possible therapeutic target for future pharmacological actions [61]. However, in none of these previous studies was enzyme inhibition associated with the lipidomic cell profile and its hypoxic context.

We previously associated the expression of cellular hypoxia markers with SCD-1 overexpression in a large number of tumor samples (*n* = 24) [27]. In this study using Caki-2 cells exposed to in vitro chemical hypoxia [47], we effectively determined that the enzyme was highly overexpressed, mimicking the microenvironmental conditions of the tumor core. In this sense, although there is evidence that SCD-1 can be modulated by posttranscriptional mechanisms involving ubiquitin proteasome-dependent and -independent

pathways [62], the levels of mammalian SCD-1 appear to be principally determined by its rate of transcription [6].

Adapting to hypoxic stress is pivotal in tumor progression and determining tumor malignancy [25,63,64]. Since cellular hypoxia increases 18:1*n*-9 production, the Caki-2 cells showed increased LD biogenesis, and this was in line with the dose of CoCl2 used. Additionally, we used an 18:1*n*-9 overload of the same cells to simultaneously demonstrate that the 18:1*n*-9 increase is responsible for the increased LD production. The accumulation of these droplets is a sign of adaptation to stress and/or increased cell confluence [65]. The lipotoxicity that the increased synthesis of lipids brings about is neutralized, at least in part, with the production of these organelles [66,67]. LD biogenesis and degradation need to be discussed in the context of the synthesis and turnover of their major components: neutral lipids. Their synthesis is driven by the availability of precursors like TAG and CE [68,69]. Thus, we next performed a lipidomic analysis, under hypoxic conditions, using GC/MS to measure the levels of fatty acids in PLs and neutral lipids. In the PL fraction, large decreases were observed only in the SFAs (16:0 and 18:0), consequent with the increased expression of SCD-1. However, we did not observe changes at the 18:1*n*-9 level in this fraction. Conversely, the increases in 18:1*n*-9 were observed in neutral lipids, mainly in CE. Since these decreases in the SFAs are not quantitatively related to the newly formed 18:1*n*-9 by SCD-1 catalysis, most of 18:1*n*-9 must be redirected to mitochondrial β-oxidation [70] or, more likely (given the hypoxic context), released to the extracellular space. Thus, all types of cancer overexpressing this enzyme show high rates of cell proliferation, and this can only be done in cell contexts with high metabolic energy production [6,8,18,19,60,71]. On the other hand, Kamphorst et al. [72] found in other tumor lines (breast, lung and cervix) that hypoxia inhibits the catalytic activity of SCD-1 (since it uses oxygen as an electron acceptor) and shows that 18:1*n*-9 could be imported from the extracellular space.

The essential role of SCD-1 in cancer cell mitogenesis was unambiguously demonstrated by several works in which the suppression of SCD-1 by genetic and pharmacological means led to a slower rate of cell proliferation and decreased survival [61]. In this study, we observed that the viability of CAY 10566-treated Caki-2 cells was strongly correlated with the degree of inactivation of SCD-1, firmly establishing a positive relationship between the rate of MUFA synthesis and cell replication. Thus, at different times, we restored 18:1*n*-9 in the cell culture and observed that the cell viability improved, compared with CAY 10566 -treated Caki-2 cells. Taken together, the lipidomic profile and cell viability experiments allowed us to assume that changes in the 18:1*n*-9 levels are critical for these tumor cells. Simultaneously, it has been observed that the excess content of long-chain fatty acids, especially SFAs, triggers programmed cell death in a process known as lipid-mediated toxicity or lipoapoptosis [73]. Thus, these two effects (SFA increase and 18:1*n*-9 decrease, caused by SCD-1 inhibition) evidently synergize and explain the cytotoxicity observed by long-term pharmacological inhibition.

Finally, the findings described here support the concept that SCD-1 may be a potentially useful target for ccRCC treatments [61,74]. The specific design of small-molecule inhibitors of SCD-1 activity could be of great potential for possible therapeutic agents [75]. Likewise, the association of SCD-1 inhibitors with therapeutic agents that target signaling pathways and their receptors (i.e., tyrosine kinase-mediated cascades, such as pazopanib, sunitinib, axitinib and cabozantinib, or temsirolimus, which targets mTOR, among others) already in use in medical oncology could also be an attractive option for future implementation [2]. Despite these hypothetical considerations, establishing the value of SCD-1 inhibitors as a protective agent for the treatment of ccRCC will require more extensive experimental testing and careful preclinical validation.

#### **5. Conclusions**

In this work, we provided evidence supporting the hypothesis that the lipid composition of ccRCC depends on the hypoxic microenvironment prevailing in certain areas, such as the center of the tumor. Our results added to this concept by demonstrating that there is

a tumor microheterogeneity in terms of the fatty acid distribution in different lipid species such as PLs, TAG and CE. In line with the above, we showed that SCD-1 is particularly influenced by hypoxia, since it is overexpressed under these conditions and catalyzes the conversion of 18:0 into 18:1*n*-9, favoring tumor cell proliferation. In addition, we provided evidence to reinforce the idea that SCD-1 is a meaningful pharmacological target pondering the global hypoxic context of the tumor microenvironment.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/cancers13122962/s1: Figure S1: Cell proliferation rates (growth constants (k) and cell-doubling times) in all the conditions tested. Data are expressed as the means ± SEM and are representative of three independent experiments.

**Author Contributions:** Conceptualization: J.P.M., J.B. and J.P.R. Methodology: J.P.M., F.M. and T.S. Software: J.P.R. and F.M. Validation: J.B. and M.A.B. Formal analysis: J.P.R. and F.M. Investigation: J.P.M. Resources: J.B., M.V.A., J.P.R. and F.M. Data curation: M.V.A. and M.A.B. Writing—original draft preparation: J.P.M., J.B. and J.P.R. Writing—review and editing: J.P.M., J.P.R., M.A.B. and J.B. Visualization: J.P.M. and J.P.R. Supervision: J.P.R., M.V.A., M.A.B. and J.B. Project administration: J.P.R., M.V.A. and J.B. Funding acquisition: J.P.R., M.A.B. and J.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Spanish Ministry of Science and Innovation, grant number PID2019-105989RB-I00, and Universidad Nacional del Nordeste de Argentina, grant number PI 18-I009. CIBERDEM is an initiative of Instituto de Salud Carlos III.

**Institutional Review Board Statement:** This study was conducted according to the guidelines of the Declaration of Helsinki and was approved by the Department of Medical Research and Bioethics Committee of the José Ramón Vidal Hospital in Corrientes, Argentina (date of approval, June 30 2015).

**Informed Consent Statement:** Informed consent was obtained from all the subjects involved in the study.

**Data Availability Statement:** The datasets used during the current study are available from the corresponding author upon reasonable request.

**Acknowledgments:** The authors thank Juan Todaro for the excellent technical advice.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of the data; in the writing of the manuscript or in the decision to publish the results.

#### **References**


## *Article* **PBRM1 Immunohistochemical Expression Profile Correlates with Histomorphological Features and Endothelial Expression of Tumor Vasculature for Clear Cell Renal Cell Carcinoma**

**Kazuho Saiga 1,†, Chisato Ohe 1,\*,†, Takashi Yoshida 2, Haruyuki Ohsugi 2, Junichi Ikeda 1,2, Naho Atsumi 1, Yuri Noda 1, Yoshiki Yasukochi 3, Koichiro Higasa 3, Hisanori Taniguchi 2, Hidefumi Kinoshita <sup>2</sup> and Koji Tsuta <sup>1</sup>**


**Simple Summary:** The PBRM1 protein, whose gene is the most frequently mutated one in clear cell renal cell carcinoma (ccRCC) following *von Hippel-Lindau*, has been proposed as a potential biomarker for ccRCC. However, the association of the PBRM1 immunohistochemical expression with histomorphological features of ccRCC and the endothelial expression of tumor vasculature, which is an important role of the tumor microenvironment related to treatment response, is little known. Recently, our research team has established a vascularity-based architectural classification of ccRCC correlated with angiogenesis and immune gene expression signatures, which could provide prognostic information and function as a surrogate for treatment selection. In the present study, we found the PBRM1 expression was correlated with the architectural patterns. Furthermore, we demonstrated that endothelial expression tended to be lost in cases with low PBRM1 expression. This correlation implied the orchestrated expression of PBRM1, raising the possibility that the cancer cells and their microenvironment interact in ccRCC.

**Abstract:** Loss of the polybromo-1 (PBRM1) protein has been expected as a possible biomarker for clear cell renal cell carcinoma (ccRCC). There is little knowledge about how PBRM1 immunohistochemical expression correlates with the histomorphological features of ccRCC and the endothelial expression of tumor vasculature. The present study evaluates the association of architectural patterns with the PBRM1 expression of cancer cells using a cohort of 425 patients with nonmetastatic ccRCC. Furthermore, we separately assessed the PBRM1 expression of the endothelial cells and evaluated the correlation between the expression of cancer cells and endothelial cells. PBRM1 loss in cancer cells was observed in 148 (34.8%) patients. In the correlation analysis between architectural patterns and PBRM1 expression, macrocyst/microcystic, tubular/acinar, and compact/small nested were positively correlated with PBRM1 expression, whereas alveolar/large nested, thick trabecular/insular, papillary/pseudopapillary, solid sheets, and sarcomatoid/rhabdoid were negatively correlated with PBRM1 expression. PBRM1 expression in vascular endothelial cells correlated with the expression of cancer cells (correlation coefficient = 0.834, *p* < 0.001). PBRM1 loss in both cancer and endothelial cells was associated with a lower recurrence-free survival rate (*p* < 0.001). Our PBRM1 expression profile indicated that PBRM1 expression in both cancer and endothelial cells may be regulated in an orchestrated manner.

**Keywords:** clear cell renal cell carcinoma; histomorphological features; PBRM1; immunohistochemistry; architectural patterns; endothelial cells

**Citation:** Saiga, K.; Ohe, C.; Yoshida, T.; Ohsugi, H.; Ikeda, J.; Atsumi, N.; Noda, Y.; Yasukochi, Y.; Higasa, K.; Taniguchi, H.; et al. PBRM1 Immunohistochemical Expression Profile Correlates with Histomorphological Features and Endothelial Expression of Tumor Vasculature for Clear Cell Renal Cell Carcinoma. *Cancers* **2022**, *14*, 1062. https://doi.org/10.3390/ cancers14041062

Academic Editors: Claudia Manini and José I. López

Received: 17 January 2022 Accepted: 16 February 2022 Published: 20 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Clear cell renal cell carcinoma (ccRCC), the most frequently diagnosed histologic subtype of adult RCC [1], is associated with a hyperangiogenic state due to the overproduction of vascular endothelial growth factor (VEGF) by loss of *von Hippel-Lindau (VHL)* gene function [2]. In addition to targeted therapy for these angiogenesis pathways such as VEGF receptor—tyrosine kinase inhibitors (TKIs) [3], novel systemic immunotherapy agents have improved patient survival in metastatic RCC [4,5]. However, predictive biomarkers for both the prognostic and therapeutic implications of RCC remain lacking in a clinical setting [6].

Recent genomic advances using exome sequencing revealed that the *PBRM 1* gene encoding the protein polybromo-1, which is a subunit of the SWI/SNF chromatin remodeling complex, is a second major ccRCC cancer gene, following the *VHL* gene [7,8]. Several studies have shown that the loss of PBRM1 protein has been confirmed as a possible biomarker for ccRCC, which is associated with adverse pathological factors and poor patient outcomes [9,10]. Subsequently, our research team presented a novel scoring system to predict recurrence after radical surgery using standard pathologic factors incorporating immunohistochemical (IHC) expression of PBRM1 [11]. Furthermore, because *PBRM1* is considered not only a key driver gene of ccRCC but also a key regulator of tumor cellautonomous immune response in ccRCC, the influence of PBRM1 loss on the response to immune checkpoint inhibitors (ICIs) has been investigated [7,12,13].

Recently, we first demonstrated that histological phenotypes, such as clear or eosinophilic types, were significantly correlated with survival outcomes and response to TKIs and ICIs in patients with ccRCC, which could be applied as a predictive marker for treatment selection [14]. Additionally, we established the vascularity-based architectural classification of ccRCC in accordance with nine architectural patterns, which corresponded to both angiogenesis and immune gene expression signatures [15]. Although the prognostic and therapeutic significance for architectural patterns of ccRCC has been shown [16,17], there is little knowledge on how genomics and subsequent protein expressions are reflected in histomorphological features [18].

To evaluate the association of the PBRM1 expression with histomorphological features, we semiquantitatively re-evaluated the expression by using the PBRM1-stained slides used in our previous study [11]. In addition, we noticed that the expression in vascular endothelial cells, which has been used as one of the internal positive controls in some studies [10,11], tended to decrease or disappear in the PBRM1 loss cases. However, there is little evidence regarding PBRM1 expression of the tumor vasculature, which plays an important role in the tumor microenvironment [19]. In the present study, we aimed to evaluate whether the histomorphological features of ccRCC correlate with the PBRM1 expression of cancer cells. Furthermore, we separately evaluated the PBRM1 expression of the vascular endothelial cells and examined the PBRM1 expression profiles of cancer cells and endothelial cells.

#### **2. Materials and Methods**

#### *2.1. Case Selection*

This study was performed under the institutional review board's approval at Kansai Medical University Hospital (No. 2018109 and No. 2020222). As in our previous report [15], data for 436 patients who underwent extirpative surgery for nonmetastatic ccRCC were identified from the institutional database between 2006 and 2017. Of these, 11 patients were excluded from this study due to an insufficient supply of pathological materials for immunohistochemistry. Thus, 425 cases with nonmetastatic ccRCC (cT1-4N0-1M0) were retrospectively analyzed. Our institutional database of RCC contains pathological findings, which were re-evaluated by a genitourinary pathologist (C.O.) based on the 2016 World Health Organization (WHO) classification [20] and the 2017 TNM staging system [21] as previously described [11,14,15]. All ccRCCs were histologically diagnosed when the carcinoma contained typical ccRCC histology and/or showed diffuse membranous positivity

of carbonic anhydrase IX by immunohistochemistry [20]. Pathological prognostic factors, including pathological TNM stage, WHO/International Society of Urological Pathology (WHO/ISUP) grade, and necrosis, were collected [22].

#### *2.2. Evaluation of Histomorphological Features*

All histomorphological features were evaluated by C.O., blinded to clinical outcomes, using whole-tissue sections of H&E-stained slides. Histological phenotype, based on cytoplasmic features, such as clear, mixed, or eosinophilic, and vascularity-based architectural classification, based on nine architectural patterns, such as compact/small nested, macrocyst/microcystic, tubular/acinar, alveolar/large nested, thick trabecular/insular, papillary/pseudopapillary, solid sheets, and sarcomatoid and rhabdoid, were determined at the highest-grade area as previously described [15].

#### *2.3. Tissue Microarray (TMA) Construction and Immunohistochemistry of PBRM1*

As previously described [11,23,24], TMA was constructed from duplicate 2 mm cores of representative tumor locations (including the highest-grade area) in each case. The morphological patterns of each core were also assessed based on the nine architectural patterns included in the vascularity-based architectural classification [15]. A primary antibody against PBRM1 (rabbit polyclonal, dilution 1:200; Atlas Antibodies AB, Bromma, Sweden) was used according to the manufacturer's protocols of the Ventana Discovery Ultra Autostainer (Roche Diagnostics, Indianapolis, IN, USA). PBRM1 was visualized with OptiView and an amplification kit (Ventana Medical System, Tucson, AZ, USA). The same PBRM1-stained slides from our previous study [11] were used in the present study. The nuclear expression of cancer cells was semiquantitatively assessed, referring to the internal positive controls (inflammatory cells or stromal fibroblasts), using the H-score. The score was determined by multiplying the staining intensity (0, none; 1, weak; 2, moderate; and 3, strong) and the percentage of positive cells (range: 0–300). The final scores (average H-score for the two cores) were determined as previously described [23]: H-score ≤ 20 was considered for PBRM1 loss, and H-score > 20 was considered for PBRM1 retention in cancer cells. An IHC evaluation was performed by two pathologists (K.S. and C.O.), and discordant cases were resolved by consensus. Next, we separately evaluated the nuclear expression of endothelial cells within the tumor area and scored them as follows: 0, none; 1, focal weak; 2, diffuse weak; or 3, diffuse strong. The scores of endothelial cells were finally stratified as PBRM1 loss (score: 0–1) and PBRM1 retention (score: 2–3). The representative PBRM1 expressions of cancer cells and endothelial cells are presented in Figure 1.

**Figure 1.** Representative PBRM1 expressions of cancer cells and endothelial cells. The staining intensity of cancer cells is assessed as follows: 0, none (internal control shows positive staining); 1, weak; 2, moderate; 3, strong. The score of endothelial cells is separately assessed as follows: 0, none; 1, focal weak; 2, diffuse weak; or 3, diffuse strong. The negative and positive expressions of endothelial cells are indicated by yellow and red arrows, respectively. Scale bar: 20 μm.

#### *2.4. Statistical Analysis*

Statistical analyses were performed using EZR version 1.54 (Saitama Medical Center, Jichi, Japan) [25]. A two-sided *p* < 0.05 was considered statistically significant. A Chi-squared test for categorical variables was used to evaluate the statistical significance among two or more groups. The t-statistic in linear regression analysis and one-way ANOVA analysis were used to evaluate the statistical significance among the architectural patterns. Interobserver agreement was statistically assessed using kappa statistics. Correlations between the two variables were evaluated using Spearman's rank correlation test. Recurrence-free survival (RFS; recurrence was calculated on imaging from the date of surgery to the date of recurrence) was assessed using the Kaplan–Meier method with the log-rank test.

#### **3. Results**

#### *3.1. Patients' Characteristics and PBRM1 Expression in Cancer Cells*

The median age of the patients was 65 years (IQR, 56–73 years). The male to female ratio was 2.8:1 (312 males and 113 females). The rate of TNM stage III or IV, WHO/ISUP grade 3 or 4, and the presence of necrosis was 24.0% (102/425), 32.3% (137/425), and 15.3% (65/425), respectively. Of the 425 patients, 57 (13.4%) experienced a recurrence of ccRCC during a median follow-up of 62.6 months (IQR, 33.8–94.0 months).

Cases with PBRM1 loss and PBRM1 retention were observed in 148 (34.8%) and 277 (65.2%) patients, respectively. The interobserver variability showed good agreement between the two pathologists (kappa = 0.84). The PBRM1 expression of clinicopathological factors is shown in Table 1.


**Table 1.** PBRM1 expression in cancer cells with clinicopathological factors in 425 cases with nonmetastatic ccRCC.

#### *3.2. Association of PBRM1 Expression in Cancer Cells with Clinicopathological Factors*

Loss of PBRM1 expression was significantly associated with worsened pathological prognostic factors, such as TNM stage, WHO/ISUP grade, and the presence of necrosis (all *p* < 0.001; Figure 2A). Regarding the association of PBRM1 expression with histomorphological features, cases with PBRM1 loss were significantly observed in the eosinophilic type, which is related to high gene expression signature scores of effector T-cells, immune checkpoint molecules, and epithelial and mesenchymal transitions [14], among other histologic phenotypes. Similarly, cases with PBRM1 loss were significantly observed in category 3, which is associated with a low gene signature of angiogenesis and high gene signatures of effector T-cell and immune checkpoint [15], among vascularity-based architectural categories (both *p* < 0.001; Figure 2B).

**Figure 2.** Association of PBRM1 expression in cancer cells with pathological factors. (**A**) Percentage of cases of PBRM1 expression and conventional pathological prognostic factors; (**B**) Percentage of cases

of PBRM1 expression and histological phenotype and vascularity-based architectural classification.

#### *3.3. Association of PBRM1 Expression in Cancer Cells with Architectural Patterns*

Regarding the association of PBRM1 expression with architectural patterns in the highest-grade area, tumors with PBRM1 loss were observed in 50/177 (28.2%) of compact/small nested, 1/36 (2.8%) in macrocyst/microcystic, 7/63 (11.1%) in tubular/acinar, 20/47 (42.6%) in alveolar/large nested, 37/55 (67.3%) in thick trabecular/insular, 10/20 (50%) in papillary/pseudopapillary, 8/9 (88.9%) in solid sheet, and 15/18 (83.3%) in sarcomatoid/rhabdoid patterns (Table 2). Representative images of PBRM1 expression in each architectural pattern are shown in Figure 3.


**Table 2.** PBRM1 expression in cancer cells with histomorphological features in 425 cases with nonmetastatic ccRCC.

**Figure 3.** Representative images of each architectural pattern and PBRM1 immunohistochemical expression. Compact/small nested, macrocyst/microcystic, and tubular/acinar patterns are highly associated with PBRM1 retention, whereas the other patterns are highly associated with PBRM1 loss. Scale bar: 20 μm.

To evaluate the correlation between architectural patterns and PBRM1 expression (Hscore), multiple linear regression analysis was performed (Figure 4). Macrocyst/microcystic (t statistic = 7.734, *p* < 0.001), tubular/acinar (t statistic = 4.228, *p* < 0.001), and compact/small nested (t statistic = 1.95, *p* = 0.0519) were positively correlated with the PBRM1 expression although compact/small nested was not statistically significant. On the other hand, thick trabecular/insular (t statistic = −5.98, *p* < 0.001), sarcomatoid/rhabdoid (t statistic = −3.829, *p* < 0.001), solid sheets (t statistic = −2.965, *p* = 0.0032), alveolar/large nested (t statistic = −2.935, *p* = 0.0035), and papillary/pseudopapillary (t statistic = −2.016, *p* = 0.0444) were negatively correlated with the PBRM1 expression.

**Figure 4.** Association of architectural patterns with PBRM1 expression in cancer cells; correlation analysis between architectural patterns in the highest-grade area and PBRM1 expression (*n* = 425). \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001 using multiple linear regression analysis.

Of 403 cases where two cores were assessed for PBRM1 expression (22 out of 425 cases were missing one core), 77 (19.1%) showed heterogeneity of PBRM1 expression (H-score ≤ 20 vs. >20) between cores. Therefore, we examined whether PBRM1 expression was correlated with the architectural patterns of the corresponding area by assessing a total of 828 cores. It was revealed that PBRM1 expression was correlated with the architectural patterns among all of the evaluated cores. Notably, this association between PBRM1 expression and the architectural patterns assessed in the highest-grade area, namely, macrocyst/microcystic, tubular/acinar, and compact/small nest, had significantly higher PBRM1 expressions (H-score) compared to the other patterns (*p* < 0.001, *p* < 0.001, and *p* < 0.05, respectively) (Figure 5).

**Figure 5.** Association of architectural patterns with PBRM1 expression in cancer cells based on H-score in each TMA core (*n* = 828). The histogram shows the mean ± standard error of the mean H-score of PBRM1 expression in cancer cells. One-way analysis of variance with the Tukey test was used for statistical analysis (N.S. means not statistically significant: \* *p* < 0.05, \*\*\* *p* < 0.001).

#### *3.4. Association between Cancer Cells and Endothelial Cells*

3.4.1. Correlation between PBRM1 Expression in Cancer Cells and Endothelial Cells

A positive correlation between PBRM1 expression in cancer cells and endothelial cells was confirmed (correlation coefficient = 0.834, *p* < 0.001; Figure 6A).

**Figure 6.** Association between cancer cells and endothelial cells. (**A**) Correlation between PBRM1 expression in cancer cells and endothelial cells. Correlations between the two variables were evaluated using Spearman's rank correlation test. (**B**,**C**) Kaplan–Meier curve of recurrence-free survival (RFS) stratified by PBRM1 expression. (**B**) PBRM1 expression of cancer cells. (**C**) PBRM1 expression of endothelial cells.

#### 3.4.2. Prognostic Significance of PBRM1 Expression in Cancer Cells and Endothelial Cells

Survival curve analysis showed that the 5-year RFS rate was significantly lower in patients with PBRM1 loss than in those with PBRM1 retained in cancer cells (71.1% versus 96.1%, *p* < 0.001; Figure 6B). Similarly, the 5-year RFS rate was significantly lower in patients with PBRM1 loss than in those with PBRM1 retained in endothelial cells (72.5 versus 95.6%, *p* < 0.001; Figure 6C).

#### **4. Discussion**

Typical histological features of ccRCC consist of neoplastic cells with clear cytoplasm and a vascular network of small and thin-walled blood vessels, activated by hypoxiainducible factors following *VHL* inactivation [20]. Although the most common architectural pattern of ccRCC is compact/small nested with an extensive vascular network, the morphologic intratumoral heterogeneity of ccRCC has been recognized [15–17]. Recent findings have shown that *VHL* mono-driver tumors are characterized by low-grade and indolent behavior with minimum intratumoral heterogeneity [26]. In contrast, tumors characterized by high-grade and aggressive behavior include multiple clonal drivers that exhibit truncal aberrations of ccRCC epigenetic-related genes: the SWI/SNF chromatin remodeling complex gene *PBRM1*, histone deubiquitinate gene *BAP1,* and histone methyltransferase gene *SETD2* [8,26]. Högner et al. also showed that the combined loss of PBRM1 and VHL may contribute to tumor aggressiveness [27]. However, little is known about the ways these genetic abnormalities impact the histomorphological features of ccRCC.

In the current study, we provided several insights into the PBRM1 IHC expression profile of ccRCC. First, we revealed the association of PBRM1 expression with histological phenotype based on cytoplasmic features [14] and vascularity-based architectural classification [15] (Figure 2B), both of which stratify patient prognosis. For histological phenotype, the eosinophilic type was significantly correlated with PBRM1 loss, followed by mixed type, whereas for vascularity-based architectural classification, category 3 was significantly enriched in the PBRM1 loss group, followed by category 2. These results indicated that PBRM1 loss was correlated with novel poor prognostic factors based on histomorphological features. Consistent with the previous reports [28–31], we showed the adverse prognostic factors of ccRCC, such as high TNM stage and WHO/ISUP grade or presence of necrosis, were significantly associated with PBRM1 loss (Figure 2A). While a study of localized RCC using TMA failed to show the prognostic role of PBRM1 loss after adjusting for the significant prognostic clinicopathological parameters [32], multivariable models of our prior study showed that PBRM1-negativity is an independent prognostic factor for RFS [11]. Thus, our findings suggest that the additional epigenetic change increases the aggressiveness of ccRCC and results in a poor prognosis.

Second, we showed architectural patterns based on a vascularity-based architectural classification [15], assessed in the highest-grade area in 425 nonmetastatic ccRCC, correlated with the PBRM1 expression profile. Macrocyst/microcystic, tubular/acinar, and compact/small nested patterns characterized by enrichment of the vascular network (corresponded to category 1) were positively correlated with PBRM1 expression, whereas alveolar/large nested, thick trabecular/insular, and papillary/pseudopapillary patterns characterized by the widely spaced-out vascular network (corresponded to category 2), or solid sheets and sarcomatoid/rhabdoid patterns characterized by scattered vascularity without a vascular network (corresponded to category 3) were negatively correlated with the PBRM1 expression (Figure 4). These results indicate that PBRM1 expression patterns differ among the architectural patterns of ccRCC with or without an extensive vascular network.

Third, in the evaluation of 828 cores considering intratumor heterogeneity, we also demonstrated architectural patterns in macrocyst/microcystic, tubular/acinar, and compact/small nested associated with significantly higher PBRM1 expression (H-score) compared to the other patterns (Figure 5), which suggested that PBRM1 expression profile correlated well with the ccRCC architectural patterns, even with intratumoral heterogeneity. Although intratumoral heterogeneity of ccRCC has been reported based on DNA sequencing and chromosome aberration analysis [33], we showed that loss of PBRM1 protein reflects morphologic heterogeneity and aggressive architectural patterns of ccRCC.

The role of PBRM1 protein expression for clinical decisions is not only being a biomarker of prognostic prediction but also providing information on molecular mechanisms and potential therapeutic targets. In the present study, we showed the prognostic predictive ability of PBRM1 loss in nonmetastatic ccRCC, while Cai et al. also showed that PBRM1 could improve the predictive accuracy for survival outcomes of metastatic RCC patients treated with tyrosine kinase inhibitors (TKIs) [34]. Recently, the effectiveness of systemic therapies (TKIs vs. ICIs) in patients with the *PBRM1* mutation status of ccRCC has also been investigated [12,13,35–37]. Although some studies have shown that patients with PBRM1 loss in ccRCC experience increased clinical benefit from ICIs [12,35], data on the effect of PBRM1 loss regarding immune responsiveness are inconsistent [13,36,37]. According to our previous study, category 3 of the vascularity-based architectural classification, which is related to loss of PBRM1 expression, was significantly associated with an inflamed and excluded immunophenotype in the localized ccRCC cohort and significantly enriched in effector-T cell and immune checkpoint gene signatures in the TCGA-KIRC cohort [15]. We have also shown that in ccRCC, including eosinophilic features related to loss of PBRM1 expression, significant clinical benefit was observed in the ICI therapy group compared to the TKI therapy group (*p* = 0.035) [14].

Contrary to our findings, however, some studies showed that *PBRM1* mutations were associated with increased angiogenesis, decreased immune infiltrates, and poor response to ICIs [13,37]. While these controversial findings have yet to be resolved, the *PBRM1* mutation does not directly determine the loss of the corresponding protein or function [38]. Because some discrepancies between *PBRM1* mutation and PBRM1 IHC expression have been reported, a comprehensive investigation, including *PBRM1* mutation, PBRM1 expression, and histomorphological features, should be conducted. Recently, Lin et al. evaluated the influence of PBRM1 loss for treatment response, focusing on the "immunogenic" tumor microenvironment [13]. However, the "non-immunogenic" tumor

microenvironment, including endothelial cells, is also an important factor for appropriate treatment strategies because combined therapies of TKIs and ICIs have been applied for metastatic ccRCC [19,39]. Nevertheless, there are a few studies focusing on the expression of PBRM1 in endothelial cells of ccRCC.

To the best of our knowledge, we are the first to have demonstrated that the PBRM1 IHC expression of endothelial cells is correlated with the expression of cancer cells, which suggests that the vascular endothelial cells may also be genetically or immunohistochemically abnormal (Figure 6). Although we should consider the possibility of a marked reduction in the protein expression due to insufficient or unequal fixation [40], positive expression of internal control such as inflammatory cells or stromal fibroblasts was confirmed in the present study (Figure 1). Angiogenesis also plays a central role in ccRCC tumorigenesis and progression, regulating the immune landscape through abnormal tumor vessel formation [39]. Our observation showed that the tumor vasculature among the vascularity-based architectural pattern of category 1 vs. categories 2 and 3 was different. The specific mechanism underlying the association of decreased PBRM1 expression with the architectural patterns without a vascular network is still unclear, but the interaction of cancer cells and endothelial cells may be suggested. In the current treatment strategies, including angiogenic therapy, the understanding of the epigenetic abnormality between cancer cells and endothelial cells should be considered. Further investigation by single-cell analysis is required to determine the mechanism of the interaction between cancer cells and endothelial cells in the tumor microenvironment.

Our current work has some limitations. The PBRM1 expression was evaluated using only TMA, including the highest-grade area. Even considering intratumoral heterogeneity, however, we showed that the PBRM1 expression was correlated with the architectural patterns. Next, we semiquantitatively assessed PBRM1 IHC expression in cancer cells using an H-score. Furthermore, we could not validate the association of architectural patterns with *PBRM1* mutation status. Despite these limitations, we comprehensively showed the association of the PBRM1 expression profile with clinicopathological factors, including detailed histomorphological features.

#### **5. Conclusions**

We demonstrated that PBRM1 expression of cancer cells correlated with histomorphological features of ccRCC and correlated with the expression of vascular endothelial cells. Our PBRM1 expression profile indicated that PBRM1 expression in both cancer and endothelial cells may be regulated in an orchestrated manner.

**Author Contributions:** Conceptualization, K.S., C.O., T.Y. and K.T.; Data curation, T.Y., J.I. and H.O.; Formal analysis, C.O., T.Y. and J.I.; Investigation, K.S. and C.O.; Methodology, C.O., T.Y., N.A. and Y.N.; Resources, H.O., Y.Y. and H.T.; Supervision, K.H., H.K. and K.T.; Visualization, C.O.; Writing—original draft, K.S. and C.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Japan Society for the Promotion of Science KAKENHI fund (Grants No. 19K16875 to C.O. and Grants No. 20K16457 to H.O.).

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of Kansai Medical University Hospital (No. 2018109 and No. 2020222).

**Informed Consent Statement:** Informed consent was obtained in the form of opt-out on the website of Kansai Medical University Hospital. No one expressed a refusal.

**Data Availability Statement:** The data are available upon reasonable request by contacting the corresponding author.

**Acknowledgments:** We are grateful to Ryosuke Yamaka for his technical assistance in the construction of tissue microarray.

**Conflicts of Interest:** C.O. received research funding from Chugai Pharmaceutical Co. Ltd. outside the submitted work. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. The remaining authors declare no conflicts of interest.

#### **References**


## *Review* **Radiogenomics in Clear Cell Renal Cell Carcinoma: A Review of the Current Status and Future Directions**

**Sari Khaleel 1, Andrew Katims 2, Shivaram Cumarasamy 2, Shoshana Rosenzweig 2, Kyrollis Attalla 2, A Ari Hakimi <sup>1</sup> and Reza Mehrazin 2,\***


**Simple Summary:** Clear renal cell carcinoma (ccRCC) is the most common type of renal cancer. As with other malignancies, knowledge of the genetic makeup of ccRCC tumors may provide insights for tumor management and outcomes. However, this normally requires obtaining tissue specimens from the tumor by invasive interventions—surgery or biopsy. Radiogenomics is a field that aims to non-invasively predict the genetic makeup of the tumor based on the tumor's appearance on conventional imaging, such as CT scans. To achieve this, radiogenomics uses complex machine learning (artificial intelligence) algorithms to process imaging data and build predictive models that can infer a tumor's genetic makeup and clinical outcomes from its features on conventional imaging. In this article, we searched scientific literature databases for radiogenomic studies in ccRCC, offering a review and critical analysis of these studies. More research and validation are needed before applying radiogenomics in clinical practice.

**Abstract:** Radiogenomics is a field of translational radiology that aims to associate a disease's radiologic phenotype with its underlying genotype, thus offering a novel class of non-invasive biomarkers with diagnostic, prognostic, and therapeutic potential. We herein review current radiogenomics literature in clear cell renal cell carcinoma (ccRCC), the most common renal malignancy. A literature review was performed by querying PubMed, Medline, Cochrane Library, Google Scholar, and Web of Science databases, identifying all relevant articles using the following search terms: "radiogenomics", "renal cell carcinoma", and "clear cell renal cell carcinoma". Articles included were limited to the English language and published between 2009–2021. Of 141 retrieved articles, 16 fit our inclusion criteria. Most studies used computed tomography (CT) images from open-source and institutional databases to extract radiomic features that were then modeled against common genomic mutations in ccRCC using a variety of machine learning algorithms. In more recent studies, we noted a shift towards the prediction of transcriptomic and/or epigenetic disease profiles, as well as downstream clinical outcomes. Radiogenomics offers a platform for the development of non-invasive biomarkers for ccRCC, with promising results in small-scale retrospective studies. However, more research is needed to identify and validate robust radiogenomic biomarkers before integration into clinical practice.

**Keywords:** radiogenomics; translational; clear cell renal cell carcinoma

#### **1. Introduction**

Renal cell carcinoma (RCC) is the most common malignant kidney tumor, accounting for approximately 85% of cases [1]. Clear cell carcinoma (ccRCC) is the most common histologic RCC subtype, particularly in advanced RCC (approximately 60–70%, and 90%, respectively) [2]. With increased use of computed tomography (CT) and magnetic resonanceguided imaging (MRI), the incidence of RCC is rising in developed countries, usually at

**Citation:** Khaleel, S.; Katims, A.; Cumarasamy, S.; Rosenzweig, S.; Attalla, K.; Hakimi, A.A.; Mehrazin, R. Radiogenomics in Clear Cell Renal Cell Carcinoma: A Review of the Current Status and Future Directions. *Cancers* **2022**, *14*, 2085. https:// doi.org/10.3390/cancers14092085

Academic Editors: Claudia Manini and José I. López

Received: 15 October 2021 Accepted: 21 November 2021 Published: 22 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the clinically localized stage [3]. Despite the advancements in cross-sectional imaging technology, their ability to differentiate RCC subtypes and their underlying molecular profiles remain limited [4,5].

One approach to improve the diagnostic ability of conventional imaging has been the adoption of advanced computational and statistical methods to process high throughput radiologic features extracted from conventional imaging, giving rise to the field of radiomics [6]. In parallel, our understanding of the genomic profiles of cancers and their potential as diagnostic, prognostic, and therapeutic biomarkers has been advanced by the application of complex computational and statistical methods to analyze high-throughput next-generation sequencing data, allowing for complex genomic, transcriptomic, and epigenomic analyses of tumor specimens. Such analyses in the field of RCC have revealed that in addition to histologic variance, RCC is a genetically diverse disease, with distinct molecular genomic and transcriptomic profiles that correlate with clinical outcomes such as recurrence, progression, and response to systemic therapies [7–10].

Despite the above advances in molecular and radiologic profiling of RCC in general and ccRCC in particular, the current prognostic models remain based on clinical, pathologic, and laboratory characteristics, with the pathologic stage heavily influencing cancer-specific survival [11–14]. The reliance of these models on pathologic staging makes them inherently invasive, requiring tissue diagnosis based on surgical extirpation or tissue biopsy, with no standardized non-invasive or pre-treatment biomarkers that can be used to classify RCC or predict tumor behavior. This limitation applies to genomic profiling tools, as well, as they also require tissue extraction for their analyses, along with complex and cost-prohibitive translational infrastructures that currently limit their applicability in clinical practice.

Radiogenomics is a novel field that circumvents the above challenges by utilizing computational machine learning algorithms to correlate radiomic features of disease (radiologic phenotype) with its underlying molecular profile (genotype), thereby offering a platform for the development of non-invasive biomarkers to aid in treatment decisions and disease [15,16]. Of note, while the term "radiogenomics" has been used interchangeably with "radiomics" in literature to describe the study of radiologic features of predictive treatment outcomes, "radiogenomics" is more commonly used to describe the study of the molecular changes underlying the radiologic phenotype of a disease process, including genetic mutations, gene expression, and methylation (epigenetic) changes [15–19].

Here, we present an in-depth review of the current state of radiogenomics in ccRCC, and examine the variety of innovative computational models that have been developed in this field to infer the molecular profile of ccRCC from its radiologic phenotype, concluding with a discussion of the field's current limitations and future directions.

#### **2. Methods**

A literature review was performed by querying the PubMed, Medline, Cochrane Library, Google Scholar, and Web of Science databases. We attempted to identify all articles pertaining to radiogenomics and ccRCC. The search terms included "radiogenomics and ... " one of the following MeSH search terms: "renal cell carcinoma", "clear cell renal cell carcinoma", "kidney cancer", or "renal cancer". Titles and abstracts of the articles retrieved from the above search were then screened for relevance. Inclusion criteria were (1) publication in the English language, (2) publication between 1/2009 and 9/2021, (3) and article topic pertaining to radiogenomics of ccRCC. Exclusion criteria included (1) publication before 1/2009, (2) not published in the English language, (3) study topics not pertaining to ccRCC or radiogenomics, (4) duplicate articles, and (5) non-primary literature, e.g., abstracts, review articles, and letters to the editor, which were excluded after being reviewed to identify any missed primary studies.

#### **3. Results**

Overall, 141 articles were identified in our initial search, of which 16 fit our inclusion criteria described above. Radiogenomic features related to mutational status were the most commonly described and targeted features for modeling (eight articles), followed by gene expression (five articles) and epigenetic features (one article). Only two articles developed clinical prognostic models utilizing radiogenomic data. Most articles focused on multiphasic, contrast-enhanced CT scan as the modality of choice, with two paper(s) discussing MRI features. A PRISMA flow chart of our search with inclusion and exclusion criteria can be seen in Figure 1. A list of the included articles along with a summary of their methodology and targeted predictive outcomes can be found in Table 1.

**Figure 1.** PRISMA flow diagram of article selection criteria.





**Table**

**1.**

*Cont.*





PFS with IO between the two groups

for angiogenesis

immune infiltration.

 or

#### *3.1. Key Genetic Mutations in ccRCC*

Key gene mutations identified in ccRCC include *VHL*, *PBRM1*, *BAP1*, *SETD2*, and *KDM5C*; most of which are located on the short arm of chromosome 3 [35]. Key genetic mutations and their radiogenomic characteristics as well as prognostic value are discussed below, and are summarized in Table 2.


**Table 2.** Summary of the top 5 most common gene mutations in ccRCC.

#### 3.1.1. *VHL*

*VHL* gene alteration is the most common mutation in solid ccRCC, with very high frequency (>90%) of biallelic inactivation due to deletion, mutation, or loss of heterozygosity [36,37]. As normal *VHL* protein complexes with other proteins to degrade hypoxiainducible factor (HIF), *VHL* loss or mutation results in constitutive activation of HIF, promoting cell growth and neo-angiogenesis through the VEGF pathway [38]. Despite its prevalence in ccRCC, the presence of a *VHL* mutation in patients with ccRCC has no prognostic value [10,36,39,40].

#### 3.1.2. *PBRM1*

*PBRM1* is the second most commonly mutated tumor suppressor gene in ccRCC (40–50%), and is often co-deleted with *VHL*. This gene encodes for a nucleosome remodeling complex which limits DNA accessibility to RNA polymerase and transcription factors [35,41]. The prognostic value of *PBRM1* mutation is unclear, with a recent meta-analysis suggesting that mutation and/or loss of in *PBRM1* is a poor prognostic factor in localized disease and a good prognostic factor in advanced disease [42,43]. Other analyses suggest that *PBRM1* mutation status may be predictive of response to immune checkpoint inhibitors [44,45]. *PBRM1* mutations are most associated with solid ccRCC on imaging [20,21].

#### 3.1.3. *BAP1*

*BAP1* gene, present on the short arm of chromosome 3, is mutated in 10–15% of ccRCC, and is typically mutually exclusive of *PBRM1* mutation [35,46]. This tumor suppressor gene encodes a ubiquitin carboxyl-terminal hydrolase that regulates with downstream

targets involved in cell breakdown and replication, with *BAP1* inactivation resulting in uncontrolled cell proliferation [41,47]. *BAP1* mutation has been associated with more aggressive disease and lower overall survival in ccRCC, with coagulative necrosis and high Furman grade on tumor pathology [48,49].

Typical radiologic features associated with *BAP1* mutation include renal vein invasion, ill-defined tumor margins, and intratumor calcifications. Of note, *BAP1* mutations were absent in multicystic ccRCC [20,21].

#### 3.1.4. *SETD2*

As with *BAP1*, *SETD2* is a tumor suppressor gene located on the short arm of chromosome 3, and is mutated in approximately 10–15% of ccRCC [35]. *SETD2* loss has been associated with poor prognosis in nonmetastatic ccRCC [48]. Radiomic analyses note *SETD2* mutation to be absent in multicystic ccRCC, with no consistent CT imaging findings predictive of *SETD2* mutation in solid ccRCC [20,21].

#### 3.1.5. *KDM5C*

*KDM5C* is mutated in approximately 6–7% of ccRCC [35]. The prognostic value of *KDM5C* remains debated, with one series noting an association with prolonged survival in metastatic ccRCC [50]. Tumors with *KDM5C* mutation were consistently associated with renal vein invasion on CT and absent in multicystic ccRCC [20,21].

#### *3.2. Overview of Radiogenomics Workflow*

As mentioned earlier, radiomics refers to the extraction and analysis of quantitative imaging features from cross-sectional imaging modalities, while radiogenomics refers to the study of the translational phenotype underlying these imaging features [51]. A typical radiogenomic workflow is shown in Figure 2. First, the region of interest (ROI), being the tumor and/or specific tumor sub-region(s), is "segmented", i.e., outlined in all slices of the imaging study using manual or semi-automated segmentation software, generating a 3D rendering of the ROI. Next, specialized software is used to extract hundreds to thousands of radiomic features from the ROI "agnostically", with no knowledge of its clinical context or molecular profile, such as malignant/benign status, RCC subtype, or mutational profile. Extracted features may include first-order statistics of voxel intensity and distribution, as well as higher level metrics of tumor shape, texture, and 2D/3D features, extracted from one or more phases of the imaging study. Next, machine learning (ML) algorithms are used to process these raw features to identify the subset of features that are predictive of an outcome of interest, which in radiogenomics would include specific gene mutation, gene expression profile, or clinical outcome [52]. The radiogenomic model constructed from this subset of features is usually "trained" using one dataset, followed by external cross-validation in an independent dataset.

**Figure 2.** Flowchart showing typical radiogenomic workflow. Using cross-sectional images, a region of interest (ROI) that contains either the whole tumor or subregions within the tumor can be identified and outlined using manual in process called segmentation, using semi-automated, or automated segmentation software. Some segmentation software, such as 3D Slicer (shown above) allow for further ROI rendering in 3D dimensions. Quantitative radiomic features are extracted from ROI using separate or built-in radiomic feature extraction modules. Finally, this data is integrated with corresponding tumor molecular profile, as well as patient clinical data. These data are then processed using machine learning algorithms to develop diagnostic, predictive, or prognostic models for outcomes of interest.

#### *3.3. Mutational Radiogenomic Biomarkers*

In this section, we review articles that develop radiogenomic models to predict tumor gene mutational profile in ccRCC, which mostly focused on the previously discussed *PBRM1* and *BAP1* mutations.

Chen et al. (2018) presented a radiogenomic predictive model to predict multiple ccRCC gene mutations (*VHL*, *PBRM1*, and *BAP1*) using quantitative CT features. To achieve this, they developed a new multi-classifier multi-objective (MCMO) model to train their model against multiple objectives (ccRCC mutations of interest) rather than a single objective. After training their model using an institutional cohort of 57 patients, it was validated using The Cancer Genome Atlas's Kidney-Renal Cell Carcinoma (TCGA-KIRC) cohort. Their model achieved prediction accuracy of 0.81, 0.78, and 0.90 for *VHL*, *PBRM1*, and *BAP1* genes, respectively, with AUC ≥ 0.86, sensitivity ≥ 0.75, and specificity ≥ 0.80 [22].

Focusing on *PBRM1* mutation, which is a likely good prognostic factor in advanced ccRCC [42,53], Kocack et al. (2019) developed two predictive radiogenomic models using artificial neural network algorithm (ANN) and RF algorithms to differentiate ccRCC tumors by *PBRM1* mutations status in the TCGA-KIRC cohort (45 patients; 29 without and 16 with *PBRM1* mutation). Their ANN model demonstrated an accuracy of 88.2% (AUC = 0.925) compared to 95.0% (AUC = 0.987) with the RF algorithm. However, they did not directly evaluate their model's correlation with clinical outcomes [24].

The same group (Kocak et al., 2020) then developed an RF-based radiogenomic model for the prediction of *BAP1* mutation status, which carries poor prognostic implications in ccRCC [20,21,25], in a subset of 65 patients from TCGA-KIRC (13 with and 52 without *BAP1* mutation). This model correctly classified *BAP1* mutation status in 84.6% of cases (AUC = 0.897) [20,21,25]. The same algorithm (RF) and dataset (TCGA-KIRC) were used by Feng et al. (2020) to also predict *BAP1* mutation status, but using different segmentation and radiomic feature extraction platforms; the model accurately classified 83% (AUC = 0.77) of *BAP1* mutation status with a sensitivity of 0.72, specificity of 0.87, and precision of 0.65 [26]. Finally, targeting *BAP1* status as well, Ghosh et al. (2015) developed four imaging phase-specific *BAP1* classifiers for the non-contrast, cortico-medullary, nephrographic, and

excretory phases of CT studies from the TCGA-KIRC cohort (78 patients). Interestingly, their model utilized 3D feature extraction to evaluate intra-tumoral heterogeneity (ITH), which they hypothesized reflected *BAP1* mutational status [27]. In contrast, none of the previously discussed studies considered ITH in their model design, despite its known prevalence and influence on clinicopathological and molecular assessment of ccRCC, as one tumor area's molecular profile may be different from another's, with downstream implications for the extracted radiomic features in these models [54,55].

#### *3.4. Beyond Gene Mutations: Transcriptomic and Epigenetic Radiogenomic Biomarkers*

As discussed earlier, the clinical relevance of some of the most common mutations in ccRCC remains unclear, particularly given the low prevalence of some of these mutations, limiting their potential as clinical biomarkers [20,21,35–39,41–50,56]. In contrast, transcriptional (gene-expression based) signatures have been shown to be better tools for classifying ccRCC into clinically-relevant molecular subtypes [57,58]. Such subgrouping classifications include clear cell type A (ccA) and clear cell type B (ccB) described by Brannon et al. using microarray data. Using these tumor classifications, they noted a prognostic difference between the two groups; ccA was significantly associated with better survival compared to ccB [8]. As transcriptomic research shifted from microarray to next-generation RNA-sequencing (RNA-Seq) technology, Brooks et al. developed a 34-gene expression signature, ClearCode34, for the classification of localized ccRCC tumors to ccA and ccB categories using RNA-Seq data [57]. Another attempt at transcriptomic profiling of ccRCC performed by the TCGA group using unsupervised clustering of RNA-Seq data identified four subgroups, m1–m4. Supervised clustering of these subgroups against ccA/ccB subgrouping noted cluster m1 to correspond to ccA, m2–m3 to correspond to cluster ccB, and m4 to correspond to the 15% of tumors that did not align with either ccA or ccB. As with ccA, the m1 subgroup had a survival advantage over m2–m4, sharing some of its genes with the *PBRM1* mutation group and functions within the chromatin remodeling process. In contrast, the m3 subgroup harbored mutations of *PTEN* and *CDKN2A*, while patients within the m4 subgroup exhibited a higher frequency of *BAP1* mutations [35]. Furthermore, these subtypes were associated with distinct radiomic features; the m1 subgroup had well-defined tumor margins (vs. ill-defined, OR = 2.104; CI 1.024–4.322), while the m3 subtype was less frequently associated with well-defined tumor margins (OR = 0.421; CI 0.212–0.834) and had more collecting system invasion (OR = 2.164; CI 1.090–4.294) and renal vein invasion (OR 2.120; CI 1.078–4.168). There were no significant CT findings with the m2 or m4 subgroups [7].

In this section, we explore radiogenomic models that correlate radiologic tumor "phenotype" to its underlying transcriptomic and epigenetic molecular profile, rather than genetic mutational profile. In addition to their ability to reflect variations in individual tumor gene expression and hypermethylation patterns, these models are potentially more applicable to clinical practice than radiogenomic models that predict only genomic mutational profile, given the ability of their targeted molecular expression profiles to better reflect survival and therapeutic outcomes [29,30].

#### 3.4.1. Transcriptomic Radiogenomic Biomarkers

Using the aforementioned transcriptomic ccA/ccB ccRCC subtype, Yin et al. developed a model utilizing radiomic features extracted from MRI/PET data to classify ccRCC into ccA or ccB subtypes, using sparse partial least squares discriminant analysis (SPLS-DA) to build two predictive models—one with the radiomic features alone, and one incorporating clinical characteristics, mRNA, microvascular density, and molecular subtype information. The correct classification rate was 87% vs. 95.6% using the radiomic signature alone vs. the combined signature, respectively [29]. However, the study utilized a small cohort (23 specimens from eight primary ccRCC patients), and PET imaging is not usually used for evaluation nor surveillance of localized ccRCC.

#### 3.4.2. Epigenetic Radiogenomic Biomarkers

At the epigenetic level, DNA methylation, particularly the runt-related transcription factor 3 (RUNX3) gene, has been correlated with overall survival [30]. Cen et al. (2019) evaluated the correlation between RUNX3 methylation levels and certain imaging features on CT in ccRCC. Among somatic CT findings, margin status (ill vs. well-defined; OR 2.685; CI 1.057–6.820) and intratumoral vascularity (present or absent; OR 3.286; CI 1.367–7.898) were significant independent predictors of high RUNX3 methylation levels on multivariate logistic regression [30].

#### *3.5. Beyond Predicting Molecular Profile: Radiogenomic Models as Clinical Biomarkers*

While the above reviewed studies present impressive analyses and methods for inferring tumor biology using radiomic features, the clinical relevance of their proposed features and models remains unproven without direct assessment of their ability to predict clinical outcomes. In this section, we review a few notable radiogenomic studies that go beyond correlating only radiomic and molecular features to also demonstrating a direct correlation between their radiogenomic biomarkers with clinical outcomes for ccRCC.

Focusing on radiologic features predictive of survival outcomes, Huang et al. performed radiogenomic analysis of CT imaging for ccRCC cases with corresponding RNA expression data in the TCGA-KIRC cohort. LASSO-COX regression was used to identify prognostic radiomic features and prognostic gene signatures. An RF algorithm was then used to combine prognostic and radiomic features into a radiogenomic prognostic model. The radiogenomic model outperformed the radiomic features-only model at predicting overall survival at 1, 3, and 5 years (average AUCs for 1-, 3-, and 5-year survival of 0.814 vs. 0.837, 0.74 vs. 0.806, and 0.689 vs. 0.751, respectively) [31].

In another study, Jamshidi et al. constructed a radiogenomic risk score (RSS) using a cohort of patients who underwent nephrectomy with corresponding micro-array-derived gene expression data. Following CT imaging feature extraction, multivariate regression was used to identify features most predictive of variation in supervised principal component (SPC) gene expression analysis. These features were used to constitute their RSS, which was validated in a separate patient cohort (70 for validation of the signature's correlation with micro-array results, 77 for correlation of signature with disease-free survival). The RRS exhibited a statistically significant correlation with micro-array SPC variation (R = 0.57, *p* < 0.001, classification accuracy 70.1%, *p* < 0.001) and disease-specific survival (log-rank *p* < 0.001), accounting for stage, grade, and performance status (multivariate Cox model *p* < 0.05, log-rank *p* < 0.001) [32]. In a separate study, the RRS was validated in a cohort of 41 mRCC patients undergoing cytoreductive nephrectomy (CRN) and pre-surgical bevacizumab, noting that it was able to stratify radiological progression-free survival (rPFS) in this cohort; patients with a low RSS vs. high RSS had longer rPFS (25 months vs. 6 months; *p* = 0.005) and OS (37 months vs. 25 months; *p* = 0.03) [33].

Focusing on micro-RNA (miRNA) expression in RCC, Marigliano et al. (2019) evaluated the correlation between a variety of radiomic features extracted from a cohort of 20 ccRCC patients, and their expression levels of selected microRNAs. Specifically, they examined the correlation of these features with miR-21-5p, miR-210-3p, miR-185-5p, miR-221-3p, and miR-145-5p, which had been shown to correlate with clinical outcomes in ccRCC [59]. They found no significant correlation between their extracted features and expression of any of their evaluated miRNAs [28].

While the molecular profiling of tumors using transcriptomic and epigenetic signatures offers more clinically meaningful biomarkers than genomic mutational signatures, it overlooks the critical role of the tumor's stromal and immune background, collectively referred to as the tumor microenvironment (TME), in the prognosis and therapeutic response of ccRCC. This role has been increasingly recognized with the rise of immunotherapy (IO) regimens, which target the immune component of the TME as monotherapy or in combination with TKI agents, which target the angiogenic component of the TME, as well [60–63].

In this regard, Udayakumar et al. (2021) utilized dynamic contrast-enhanced MRI (DCE-MRI) imaging to identify areas of high and low colocalized enhancement within tumor regions of 49 ccRCC patients undergoing DCE-MRI prior to nephrectomy, followed by targeted sampling and RNA-sequencing of nephrectomy specimen regions corresponding to these areas. They found enhancement-high tumors to exhibit upregulated angiogenesisrelated TME gene expression signatures, while enhancement-low areas exhibited higher levels of immune (T-cell infiltration) TME signatures, confirmed by immunohistochemical analysis. They then validated their model's ability to predict response to TKI or immunotherapy (IO) treatments in a cohort of 19 patients with metastatic ccRCC, noting better PFS with TKI in the enhancement-high compared to enhancement-low tumor groups (adjusted *p* < 0.0001), but no significant difference in PFS with IO between the two groups [34].

#### **4. Discussion**

In this review, we provided an overview of radiogenomic studies in ccRCC, the most common subtype of RCC, and renal malignancies in general. While the majority of studies focused on developing models for the prediction of tumor gene mutational profiles, we noted a shift towards the prediction of gene expression patterns and epigenetic changes within the tumor as well as the tumor microenvironment, which provide better insights into tumor biology and potential therapeutic response than isolated gene mutation profiles. A minority of the reviewed models were also shown to be predictive of relevant clinical outcomes, such as cancer-specific survival and response to systemic therapy in advanced ccRCC. Such models may complement the management of localized renal tumors to confirm whether the tumor exhibits high- or low-risk features that may warrant more aggressive management vs. surveillance, and in advanced ccRCC to determine the optimal systemic treatment regimens based on radiogenomic assessment of the tumor and its microenvironment.

However, the clinical applicability of these models remains limited by several factors. First, all the predictive models presented by the reviewed studies were developed using relatively small cohorts, mostly utilizing the same publicly available cohort (TCGA-KIRC), potentially overfitting their models to this cohort, with only a few performing external validation in independent cohorts. Second, the quality of CT studies is dependent on a variety of technical factors, such as the CT scanner, acquisition mode, and voxel reconstruction algorithms, thereby affecting the quality of extracted radiomic data. Third, the extracted radiomic features come from segmented tumor images, which are usually manually or semi-automatically delineated by a human user—a process that is inherently subjective and liable to inter-observer variability. Fourth, there are no standardized protocols or software tools for radiomic feature extraction, with the concern that the hundreds to thousands of radiomic features extracted by one software package are often redundant and difficult to replicate by other software packages [64], thus limiting the external validity of the models developed from these features. The Image Biomarker Standardization Initiative is a recent attempt at addressing this issue, establishing a standardized set of unique radiomics features [65], although compliance with this initiative has yet to be seen in radiogenomics publications. This lack of a unified radiomic feature extraction protocol or terminology limits our ability to compare the subsets of predictive radiomic features across different models, which consequently limits the ability to identify any consistent radiomic features across different models. Furthermore, it hinders attempts to identify the biologic processes that may underlie changes in these radiomic features. Fifth, most of the models did not consider intra-tumoral heterogeneity, despite its known influence on clinicopathological and molecular assessment of ccRCC, with different tumor regions expressing different pathologic phenotypes and molecular profiles, with implications for therapeutic response. Therefore, a radiogenomic model that was trained to treat the entire tumor region as a single homogenous entity may not accurately predict a tumor's molecular profile or its correlated clinical outcomes. Finally, while the ultimate measure of any biomarker is to show reliable and independent correlation with clinical outcomes, complementing standard-of-care biomarkers and predictive tools, most of these studies focused on developing models to predict molecular profiles without directly demonstrating clinical relevance as an independent biomarker of key prognostic and therapeutic outcomes, or in combination with established predictive models and nomograms. These are critical limitations that must be addressed for radiogenomics to be reliably used as a tool in clinical practice.

Despite these limitations, the above studies demonstrate the potential of radiogenomics as a non-invasive biomarker of tumor biology, utilizing complex computational tools to identify radiologic tumor features that correlate with genomic, transcriptomic, and/or epigenetic features of the tumor, and their downstream clinical implications.

#### **5. Conclusions**

The field of radiogenomics is a potentially promising tool in constructing personalized cancer care, offering a novel non-invasive translational biomarker that can be used for molecular profiling of clear cell renal carcinoma. However, this field remains relatively immature, and all the reviewed studies in the field rely on retrospective analyses, with no large-scale prospective trials, a critical requirement for the implementation of this technology in clinical practice.

**Author Contributions:** Conceptualization, A.K. and R.M; methodology, A.K., S.C. and S.K.; writing original draft preparation, S.K. and A.K.; writing—review and editing, S.K., A.K., K.A., A.A.H. and R.M.; visualization, S.K. and S.R.; supervision, A.A.H. and R.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Prognostic Gene Expression-Based Signature in Clear-Cell Renal Cell Carcinoma**

**Fiorella L. Roldán 1,2,†, Laura Izquierdo 1,2,†, Mercedes Ingelmo-Torres 1,2, Juan José Lozano 3, Raquel Carrasco 1,2, Alexandra Cuñado 1, Oscar Reig 4, Lourdes Mengual 1,2,5,\*,‡ and Antonio Alcaraz 1,2,‡**


**Simple Summary:** In this study, we identified molecular markers for disease progression from ccRCC tissue samples. Using the selected biomarkers and clinical data from the TCGA cohort, we developed a gene expression-based signature which enhances the prognostic prediction of clinicopathological variables and could help to provide personalized disease management.

**Abstract:** The inaccuracy of the current prognostic algorithms and the potential changes in the therapeutic management of localized ccRCC demands the development of an improved prognostic model for these patients. To this end, we analyzed whole-transcriptome profiling of 26 tissue samples from progressive and non-progressive ccRCCs using Illumina Hi-seq 4000. Differentially expressed genes (DEG) were intersected with the RNA-sequencing data from the TCGA. The overlapping genes were used for further analysis. A total of 132 genes were found to be prognosis-related genes. LASSO regression enabled the development of the best prognostic six-gene panel. Cox regression analyses were performed to identify independent clinical prognostic parameters to construct a combined nomogram which includes the expression of *CERCAM*, *MIA2*, *HS6ST2*, *ONECUT2, SOX12, TMEM132A,* pT stage, tumor size and ISUP grade. A risk score generated using this model effectively stratified patients at higher risk of disease progression (HR 10.79; *p* < 0.001) and cancer-specific death (HR 19.27; *p* < 0.001). It correlated with the clinicopathological variables, enabling us to discriminate a subset of patients at higher risk of progression within the Stage, Size, Grade and Necrosis score (SSIGN) risk groups, pT and ISUP grade. In summary, a gene expression-based prognostic signature was successfully developed providing a more precise assessment of the individual risk of progression.

**Keywords:** gene expression; clear-cell renal cell carcinoma; disease progression; prognostic factors; biomarkers; RNA sequencing

#### **1. Introduction**

Renal cell carcinoma (RCC) ranks third among the urological cancers with the highest incidence. Over 431,000 new cases and more than 170,000 RCC-related deaths were reported worldwide last year [1]. Clear-cell RCC (ccRCC) is the most common histological subtype and has the worst prognosis among all RCCs [2]. Currently, most of the newly diagnosed ccRCC cases are organ-confined tumors; however, after curative treatment, up to 30–40%

**Citation:** Roldán, F.L.; Izquierdo, L.; Ingelmo-Torres, M.; Lozano, J.J.; Carrasco, R.; Cuñado, A.; Reig, O.; Mengual, L.; Alcaraz, A. Prognostic Gene Expression-Based Signature in Clear-Cell Renal Cell Carcinoma. *Cancers* **2022**, *14*, 3754. https:// doi.org/10.3390/cancers14153754

Academic Editors: Jose I. Lopez and Claudia Manini

Received: 20 May 2022 Accepted: 22 July 2022 Published: 1 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

will develop tumor metastases. Unfortunately, metastatic patients have a very poor fiveyear survival rate, varying between 0–20% [3,4].

According to the European Urological Guidelines, the standard management for all localized ccRCCs receiving surgery is limited to radiographic surveillance. No adjuvant treatment is approved for patients with a higher risk of progression [2]. Moreover, all surveillance recommendations are only based on clinical parameters, even though these have been proven insufficient to accurately predict disease progression, either to select ccRCC patients for adjuvant treatments or guide disease management [5,6].

Gene expression profiling has been used extensively in cancer research and has led to the discovery of new molecular prognostic markers and potential therapeutic targets. Several genetic models have been proposed in ccRCC [7–10]; however, none of those classifiers have been widely accepted nor implemented in routine clinical practice. Biomarker research for ccRCC still faces multiple challenges, mainly due to tumor heterogeneity and lack of validation studies. In addition, the use of high-throughput assays and the identification of a significant number of markers in a relatively small number of patients increase the complexity of data analysis [3]. The currently validated gene signatures comprise a large number of biomarkers, hindering their applicability and reproducibility [8,9]. Therefore, in this study we sought to develop a novel and high-performing gene expression-based signature using data generated from our cohort and The Cancer Genome Atlas (TCGA) cohort, to provide a more accurate assessment of the individual risk of progression for patients with localized ccRCC.

#### **2. Materials and Methods**

#### *2.1. Patients, Datasets Sources and Study Design*

This study was split into a three-stage approach: an initial molecular profiling, a selection and verification of prognosis-related genes and a signature development phase (Figure 1). The initial molecular profiling phase included a total of 26 localized ccRCCs who underwent partial or radical nephrectomy between 2001 and 2010 in our center (Hospital Clinic of Barcelona, Barcelona, Spain). These 26 cases consisted of 13 progressive and 13 nonprogressive patients and met the following criteria: no neoadjuvant or adjuvant treatment, no prior or concomitant malignancies or a history of inherited von Hippel-Lindau disease, all patients had thoracoabdominal CT scan staging within two month before surgery to ensure organ-confined disease. Tumors were considered progressive when local relapse or distant metastasis developed during follow-up. All patients were followed up postoperatively according to the European Urology guidelines. Any progressive patient within two months of surgery was excluded from the study. Non-progressive patients had a minimum follow-up of 10 years to ensure their status as appropriate controls. Tissue samples were obtained under institutional review board-approved protocols (HBC/2016/0333).

The selection of prognosis-related genes and signature development phases was carried out using The Cancer Genome Atlas (TCGA) dataset. Level 3 RNAseq expression data and the corresponding clinical data from TCGA ccRCC samples were obtained from the portal (https: //firebrowse.org (accessed on 8 October 2021) [11] (Supplementary Material Table S1). Survival data were obtained from portal (https://www.sciencedirect.com/science/article/pii/S0 092867418302290?via%3Dihub#app2 (accessed on 8 October 2021) [12]. After selecting samples matching our selection criteria and excluding patients without survival status or missing clinical data, a total of 356 ccRCC samples from TCGA, 68 progressive and 288 nonprogressive, were selected and the gene expression of 20,532 genes was downloaded.

**Figure 1.** Flowchart of the whole study. Abbreviations: ccRCC, clear-cell renal cell carcinoma; DEGs, differentially expressed genes.

#### *2.2. Tissue Specimens and RNA Isolation*

Formalin-fixed paraffin-embedded (FFPE) tissue blocks were reviewed. The tumor area was macro-dissected from slides (total thickness 80 μm) and RNA was isolated from FFPE specimens using the kit RecoverAll™ Total Nucleic Acid Isolation for FFPE (Ambion, Inc. Austin, TX, USA), following manufacturers' instructions. RNA was quantified by spectrophotometric analysis at 260 nm (NanoDrop Technologies, Wilmington, DE, USA) and RNA integrity was assessed using Agilent 2100 Bioanalyzer System.

#### *2.3. Molecular Profiling by RNA Sequencing*

Library preparation and sequencing method: Following rRNA removal (Ribo-Zero® rRNA Removal Kit, Illumina), RNA from 26 selected ccRCC samples was processed for library preparation using the TruSeq® RNA Access Library Preparation Kit (Illumina, San Diego, CA, USA) that allows generating libraries starting from degraded RNA. Briefly, cDNA strands were synthetized from input RNA in order to be adaptor-tagged, labeled and amplified. cDNA was then pooled and enriched by a double step of probes hybridization. The enriched targets were captured by streptavidin labeled beads, cleaned up and amplified to obtain the final multiplexed libraries. The libraries were then sequenced on an Illumina HiSeq® 4000 platform (Illumina®).

Read alignment and differential gene expression analysis: Paired-end RNA-Seq FASTQ files were trimmed from a 3 end to a fixed length based on the Phred quality score (trimmed if score fell below 20, with a minimum read length of 25) [13]. Trimmed RNA-seq reads were aligned to the GRCh38 reference genome with STAR [14] and gene counts were determined using quantMode GeneCounts. Trimmed reads were then aligned using STAR. We used limma-voom transformation and cyclic-loess to normalize the non-biological variability. An assessment of differential expression between groups was evaluated using moderated t-statistics [15].

Significant DEGs between progressive and non-progressive patients were identified based on an adjusted *p*-value of <0.05 and a fold change (FC) ≥ ±2. The heatmap and statistical analyses were performed using the R statistical package (v3.3.2). Gene set enrichment analysis (GSEA) was performed using GSEA2-2.2.0 software for testing specific gene sets based on Gene Ontology (GO) Biological Processes [16]. The "EnrichmentMap" plug-in of Cytoscape was used to create an enrichment map of the GSEA results, depicting the overlap among pathways, with similar biological processes grouped together as subnetworks [17–19]. A conservative overlap coefficient (0.5) was used to build the enrichment map. The "AutoAnnotate" plug-in identified clusters in an automated manner, visually annotating them with a summary label [20]. RNAseq files and clinical information were deposited in the Gene Expression Omnibus (GEO) with accession number GSE175648.

#### *2.4. Selection and Verification of Prognosis-Related Genes*

The DEGs identified in the previous phase were intersected with the TCGA gene expression dataset and overlapping genes were used for further analysis. The raw counts of RNA-sequencing data from the TCGA cohort were normalized using log2-based transformation. This normalized expression was used to build a multigene signature panel. Firstly, we performed univariate Cox regression analysis (considering 637 genes and 3 clinical variables) to identify the potential prognosis-related variables. Then, LASSO (Least Absolute Shrinkage and Selection Operator) regression was applied to construct the gene-based signature for predicting tumor progression in ccRCC using the cvr.glmnet function from ipflasso R package using ten-fold cross-validation and repeated it ten successive runs to increase reliability and robustness [21]. In the machine learning procedure, we fixed the three clinical variables. For all statistical analysis, a *p*-value < 0.05 was considered significant. All statistical analyses were performed using SPSS 19.0 (Statistical Product and Service Solutions; IBM Corporation, Armonk, NY, USA) and R version 3.4.2 (R Foundation for Statistical Computing, Vienna, Austria).

#### *2.5. Development of a Gene Expression-Based Signature*

Based on the expression of each gene discovered and the three clinical variables, each patient's risk score (RS) was calculated according to the risk score model. The risk score model was then used to evaluate the ccRCC prognosis according to the general form RS = exp Σ*βixis*, where *i* = 1, *k* index variables, *βi* represents the coefficient for each variable estimated from the Cox regression model, and *xis* the corresponding value for each variable in a given patient. RS was subjected to a Receiver Operating Characteristics (ROC) curve analysis to choose the most appropriate threshold for predicting tumor progression. Thereafter, Kaplan–Meier curves were generated using the selected cut-off point and compared according to the log-rank test.

The endpoints were disease progression, defined as any local relapse or distant metastasis demonstrated by radiological imaging, and cancer-specific survival. We investigated the role of the gene panel alone, clinicopathological variables alone, and a combined model including gene expression and clinicopathological variables as potential predictors of disease progression and cancer-specific survival.

#### *2.6. Pathway Enrichment Analysis*

Ingenuity Pathway Analysis (IPA) software was used to identify interactions and networks between the prognostic markers included in our gene signature, possible altered canonical pathways, regulators, diseases and functions based on direct/indirect and experimental targets.

#### **3. Results**

#### *3.1. Clinical Features of the Cohort*

The clinicopathological characteristics of patients divided by study phase are summarized in Table 1. The median (range) follow-up of the cohort was 45.6 (2.1–135.8) months. During the follow-up period, a total of 68 patients (19.1%) developed tumor progression and a total of 36 patients died of ccRCC. The median time to tumor progression and cancer-related death was 19.9 (2.1–125.5) and 48.5 (2.5–151.2) months, respectively.

**Table 1.** Demographic and pathological characteristics of enrolled patients.


\* Stage, Size, Grade and Necrosis (SSIGN) score [22].

#### *3.2. Molecular Profiling of ccRCC Samples*

Overall, we identified 1380 transcripts that were differentially expressed (*p* < 0.05) between progressive and non-progressive ccRCC samples. Of these, 639 were proteincoding genes; 217 were downregulated and 422 upregulated in progressive compared with non-progressive cases. A heat map based on the most DEGs between the two groups of ccRCC patients is shown in Figure 2A. Gene set enrichment analysis (GSEA) identified several enriched biological processes, such as the dependent toll-like receptor signaling pathway, metabolic process and immune response regulating cell surface receptor signaling pathway (Figure 2B). The full list of GO biological processes is available in Supplementary Materials (Table S2). To aid interpretation of these enriched pathways, we used enrichment maps to create a network-based representation of our results. The most prominent cluster of significantly enriched pathways recapitulated changes in the Catabolism Biological Process and Toll Signaling Pathway (Figure 2C).

**Figure 2.** DEGs in the discovery phase and gene-set enrichment analysis. (**A**) Heat map displaying the 50 most DEGs between progressive and non-progressive localized ccRCC patients. Red pixels correspond to upregulated genes, whereas green pixels indicate downregulated genes. (**B**) GSEA shows positive correlation of DEGs in biological processes involved in tumor progression. (**C**) Enrichment map where nodes represent gene sets (pathways) and edges (blue lines) denote overlapping genes between 2 pathways. Node size denotes gene set size. Predicted pathways are grouped as circles, where shades in red correspond to up-regulated gene-sets and shades in light blue correspond to down-regulated gene-sets. Highly redundant gene sets are grouped together as clusters. Abbreviations: DEGs, differentially expressed genes. GSEA, gene set enrichment analysis.

#### *3.3. Identification of Prognosis-Related Genes in an External Data Set*

To validate the 639 genes identified as DEGs in the previous study phase, these genes were intersected with the 20,532 genes from the TCGA cohort (Figure 1). As a result, we obtained 637 overlapping genes. Of those, univariate Cox regression analysis identified 132 prognosis-related DEGs and three clinicopathological variables (pT stage, tumor size and ISUP grade).

LASSO regression analysis was used to select the best combination of genes significantly associated with disease progression and to build a six-gene signature. The expressions of five of these genes: *CERCAM*, *HS6ST2*, *ONECUT2*, *SOX12* and *TMEM132A* were upregulated, while *MIA2* expression was downregulated in progressive compared with non-progressive cases. According to Ingenuity Pathway Analysis (IPA), these six validated genes were enriched in cancer, organismal injury, abnormalities, cell-to-cell signaling and interactions, cell-mediated immune response, cellular development, cellular growth and proliferation, carbohydrate metabolism, and angiogenesis, among others. Significant IPA canonical pathways are depicted in Supplementary Material Table S3. The network generated shows that there were no direct interactions between the six prognostic genes (Figure S1).

Gene expression values for each selected gene were used for Cox regression analysis. High expression of *CERCAM*, *HS6ST2*, *ONECUT2*, *SOX12* and *TMEM132A* and low expression of *MIA2* related to poor outcomes for progression-free survival (Figure S2) and cancer-specific survival (Table 2). Moreover, the clinicopathological variables pT stage, tumor size and ISUP were also found as to be prognostic factors for both survival endpoints.

**Table 2.** Univariate Cox regression analysis of statistically significant genetic and clinical variables in the validation set (TCGA cohort).


#### *3.4. Development of a Prognostic Signature*

The risk score (RS) for disease progression was calculated for each patient according to a mathematical algorithm containing the six-gene expression values; pT stage, tumor size and ISUP grade (Supplementary Material Table S4). An ROC analysis of this combined gene expression–clinicopathological model was performed and allowed the selection of a threshold of 0.789 (sensitivity 90% and specificity 60%) and 0.799 (sensitivity 94% and specificity 55%) to categorize patients into high- and low-risk groups for tumor progression and cancer-related death, respectively. The Kaplan–Meier curve of the combined generated gene expression-based model was able to discriminate two groups with significantly different probabilities of tumor progression (hazard ratio (HR) 10.79; 95%, *p* < 0.001) and cancer-specific survival (HR 19.27; 95%, *p* < 0.001) (Figure 3). Notably, the performance of the combined gene expression-based model (Area under the Curve [AUC] 0.824) was higher than that of clinicopathological variables alone (AUC 0.766) or gene expression data alone (AUC 0.753) (Figure S3).

**Figure 3.** Kaplan–Meier curves of the combined gene expression-based model for (**A**) disease progression-free survival and (**B**) cancer-specific survival for TCGA cohort.

#### *3.5. Correlation Analysis of the RS with Clinical Characteristics for Disease Progression*

Given the clinical significance of the RS in ccRCC, we sought to investigate the potential correlation between RS and clinical features. The Mann–Whitney test revealed that higher RSs correlated with a higher risk group within the SSIGN model, higher pT stage and higher ISUP grade (Figure 4). Furthermore, the Kaplan–Meier curve indicated that our established RS was capable of identifying ccRCC patients at the highest risk of progression within the groups stratified by SSIGN, pT stage and ISUP grade (all *p* < 0.05; Figure S4).

**Figure 4.** Box plots for the correlation analysis of RS with clinical characteristics for disease progression. (**A**) SSIGN risk groups, (**B**) pT stage and (**C**) ISUP grade.

#### **4. Discussion**

Currently, clinicopathological variables are the most valuable tool for predicting disease outcomes in ccRCC. However, due to the highly variable behavior of ccRCC, the prediction of tumor progression is still an important clinical challenge. For decades, surgery has remained the only treatment approved for localized ccRCC [2]. At present, disease management of these patients at high risk of progression is changing and new adjuvant treatments are being considered, opening a door to more personalized medicine in localized renal carcinoma [23,24].

Molecular profiling helps our understanding of the molecular mechanism underlying ccRCC and affords great potential to identify new biomarkers of clinical utility [25,26]. However, intratumor heterogeneity and the methodology for sample processing, readout and expression normalization have been a strong challenge in the development of a robust gene signature. Next-generation sequencing is the most advanced technique for gene expression profiling [27]; here, we used this technology to analyze the entire transcriptome profiling of tissue samples from progressive and non-progressive ccRCC patients. We then established a gene signature for predicting disease progression based on the gene expression and clinical data obtained from the TCGA cohort. We improved the gene selection and accuracy of the model by using LASSO regression analysis; this allowed us to include all DEGs found in the training set and avoid the preselection and validation of only a subset of these DEGs [28].

This study demonstrated that our gene expression-based signature was able to identify localized ccRCC patients with high and low risk of disease progression in the whole cohort and within the SSIGN risk groups. It properly correlated with clinical parameters and was proven to enhance the predictive value of the current clinicopathological variables. Furthermore, it was also predictive of cancer-specific survival. Therefore, our signature may constitute an important step forward in treatment decisions for ccRCC patients.

Remarkably, the developed gene-based panel demonstrated a greater value for prognostic prediction (HR 10.79, *p* < 0.001; AUC = 0.824) compared with similar, previously described models. Dai et al. proposed a four-gene signature with an HR of 3.1, *p* < 0.001, whereas Zhao et al. described a 15-gene model with an AUC of 0.737 [29]. Unfortunately, the different designs and methodologies of several other studies thwart any performance comparisons with their proposed genetic models. Thus, Brook et al. [9] assessed the performance of ClearCode34, which classifies ccRCC into ccA and ccB subtypes, ccB presented tumor relapse more frequently (HR 2.1; *p* = 0.001), whereas Rini et al. [8] generated a 16-gene signature associated with tumor recurrence with an HR per 25-unit increase in score of 3.37 (*p* < 0.001). Likewise, other authors have built molecular signatures aiming to predict overall survival (OS), thus making the classifiers' performance non-comparable [30–32].

Biologically, the genes from our panel are unrelated to each other and many of them have been shown to have either prognostic or biologic relevance in tumor metastasis development. According to previous reports, some of our selected genes were consistent with previously discovered biomarkers; therefore, we have further validated their value for ccRCC progression. Briefly, *CERCAM* (cerebral endothelial cell adhesion molecule) is an adhesion molecule found to be an unfavorable prognostic marker in several tumors [33–35]. Its overexpression promotes cell viability, proliferation and invasion, it is involved in the PI3K/AKT pathway [36] and is part of an immune prognostic signature for colon and rectal cancer [35]. *HS6ST2* (heparan sulfate D-glucosaminyl 6-*O*-sulfotransferase-2) is a glycolysis-related gene and has also been associated with poor disease outcomes in numerous malignancies [10,37]. Interestingly, our group previously validated this gene as an independent prognosis biomarker in intermediate/high-risk ccRCC and found it to be associated with angiogenesis, epithelia–mesenchymal transition (EMT), and indirectly related to the PD-1, PDL-1 cancer immunotherapy pathway [37–39]. *MIA2* (melanoma inhibitory activity 2) has been found in several malignancies and can act as a tumor suppressor [40] or as a proto-oncogene depending on the receptor-related signaling differences [41]. We found *MIA2* to be downregulated in progressive ccRCC. This is congruent with the human protein atlas findings, where high expression was a favorable prognostic factor in renal cancer [33,42]. *ONECUT2* (One cut domain family member 2) is a transcription factor able to activate oncogenic pathways and lineage-specific genes; hence, it is involved in EMT, angiogenesis, neural differentiation, proliferation, extracellular matrix organization, cell

locomotion and migration, among others. Overexpression of *ONECUT2* has been described in several tumors and is related to poor prognosis [43–46]. *SOX12* (Sex-determining region Y-box12) is a transcription factor, its upregulation promotes tumor progression and it is involved in EMT, apoptosis and cell proliferation [47–49]. It functions as an oncogeneregulating Wnt/B-catenin signaling to promote the growth of multiple myeloma cells [50]. As for *TMEM132A*, few studies were found in the literature, so further investigations are required to establish its role in tumor development.

Our study has multiple strengths. The first advantage of our gene expression-based signature is that it contains a low number of genes, making its clinical application easier. Our model did match genes from previous models and some of them exceeded the mere field of ccRCC, highlighting the prognostic power of the selected genes. The fact that they are involved in different carcinogenic mechanisms confers an advantage to our signature compared with others that only target one single pathway [10,51]. The high-throughput technology used to analyze the samples and the statistical methodology makes our gene model a reliable tool for predicting disease progression in ccRCC and adds important prognostic information to the clinicopathological parameters. However, we acknowledge that this study has several limitations. First, the retrospective design and the relatively small sample might have influenced our findings. Second, the definition of CSS in the TCGA cohort should be taken with caution since it was estimated [12]. Finally, despite the good performance of our six-gene model, further validations in larger cohorts are required.

#### **5. Conclusions**

A gene expression-based prognostic signature to predict disease progression in ccRCC was successfully developed; it could discriminate two groups with different probabilities of tumor recurrence. In addition, our model was also useful in predicting cancer-specific survival. The combination of genetic and clinical information enhanced the current risk stratification of the localized ccRCC patients. Refining prognostic algorithms could help to improve the disease management and follow-up of ccRCC patients.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cancers14153754/s1, Figure S1. Network of gene-gene interactions for the six genes from the gene expression-based signature. Figure S2. Kaplan–Meier curves of each gene of the model evaluating progression free survival (data from the TCGA cohort). Figure S3. ROC curves showing the performance of the gene expression-based signature. The figure depicts that the combined gene expression-based model outperforms clinicopathological variables alone or gene expression data alone. Figure S4. Kaplan–Meier curves of the combined gene expression-based model to evaluate progression-free survival within (A) SSIGN risk groups (low, intermediate and high), (B) pT stage groups (pT1/2 and pT3/4) and (C) ISUP grade groups (ISUP 1/2 and ISUP 3/4) (data from the TCGA cohort). Table S1. Clinical database of TCGA cohort used in this study. Table S2. List of GO biological processes from 1380 differentially expressed transcripts between progressive and non-progressive ccRCC samples. Table S3. Significant IPA canonical pathways from the six genes of the signature. Table S4. Risk Score for disease progression.

**Author Contributions:** Conceptualization, A.A., L.M., F.L.R. and L.I.; methodology O.R. and F.L.R.; software and validation, F.L.R. and J.J.L.; formal analysis, F.L.R., L.I., L.M. and A.A.; investigation, F.L.R., L.I. and L.M.; resources A.C., M.I.-T. and R.C.; data curation J.J.L. and F.L.R.; writing—original draft preparation, F.L.R. and L.I.; writing—review and editing, all authors contributed to reviewing and editing the writing; visualization, F.L.R. and L.M.; supervision, L.M., L.I. and A.A.; project administration, A.C., M.I.-T. and R.C.; funding acquisition, F.L.R. and A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by an Emili Letang grant from the Hospital Clínic de Barcelona to FLR/2019.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of Hospital clinic of Barcelona (protocol code HBC/2016/0333 and date of approval was 24 January 2019).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The data presented in this study are available in the article and supplementary material. RNAseq files and clinical information were deposited in the Gene Expression Omnibus (GEO) with accession number GSE175648. Further details can be obtained on request to the corresponding author.

**Acknowledgments:** We are indebted to the IDIBAPS Biobank for sample and fecha procurement. We thank Azucena Salas and Lluís Revilla for IPA software support. This work was developed at the building Centre de Recerca Biomèdica Cellex, Barcelona. We thank Helena Kruyer the English revision of the manuscript.

**Conflicts of Interest:** The authors have no conflict of interest to disclosure.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Cancers* Editorial Office E-mail: cancers@mdpi.com www.mdpi.com/journal/cancers

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-5518-8