Purpose: The currently used prognostic models for patients with nonmetastatic clear cell renal cell carcinoma (ccRCC) are based on clinicopathologic features and might be improved by adding molecular markers. Epigenetic alterations occur frequently in ccRCC and are promising biomarkers. The aim of this study is to identify prognostic promoter methylation markers for ccRCC.

Experimental Design: We integrated data generated by massive parallel sequencing of methyl-binding domain enriched DNA and microarray-based RNA expression profiling of 5-aza-2′-deoxycytidine–treated ccRCC cell lines to comprehensively characterize the ccRCC methylome. A selection of the identified methylation markers was evaluated in two independent series of primary ccRCC (n = 150 and n = 185) by methylation-specific PCR. Kaplan–Meier curves and log-rank tests were used to estimate cause-specific survival. HRs and corresponding 95% confidence intervals (CI) were assessed using Cox proportional hazard models. To assess the predictive capacity and fit of models combining several methylation markers, HarrellC statistic and the Akaike Information Criterion were used.

Results: We identified four methylation markers, that is, GREM1, NEURL, LAD1, and NEFH, that individually predicted prognosis of patients with ccRCC. The four markers combined were associated with poorer survival in two independent patient series (HR, 3.64; 95% CI, 1.02–13.00 and HR, 7.54; 95% CI, 2.68–21.19). These findings were confirmed in a third series of ccRCC cases from The Cancer Genome Atlas (HR, 3.60; 95% CI, 2.02–6.40).

Conclusions: A four-gene promoter methylation marker panel consisting of GREM1, NEURL, LAD1, and NEFH predicts outcome of patients with ccRCC and might be used to improve current prognostic models. Clin Cancer Res; 23(8); 2006–18. ©2016 AACR.

Translational Relevance

Of the patients diagnosed with localized clear cell renal cell carcinoma (ccRCC), about 30% relapse despite curative surgery. Currently, the best predictors of patient outcome are models that combine clinical parameters, such as performance status, with pathological features such as tumor–node–metastasis (TNM) stage, tumor size, and Fuhrman grade. However, with increased knowledge of the biology of ccRCC development and progression, these models might not fully represent tumor biology and therefore fall short in adequately predicting outcome of the individual patient. We show that promoter methylation of a panel of four genes, that is, GREM1, NEURL, LAD1, and NEFH, was associated with significantly worse survival in two independent patient series. These findings were validated in a third series of ccRCC cases from The Cancer Genome Atlas. Integration of this four-marker panel with current prognostic models might improve their predictive capability and help to accurately determine ccRCC prognosis.

Clear cell renal cell carcinoma (ccRCC) is the most common adult renal neoplasm and its incidence is still increasing (1). Approximately 30% of patients with ccRCC have metastatic disease at time of diagnosis, whereas in another 30%, recurrence develops after complete resection of the primary tumor (2). Once metastatic, the survival of patients with ccRCC is poor with a 5-year survival rate of 12% (3), and current pharmaceutical interventions have limited success. Several prognostic models for ccRCC have been developed, such as the University of California Los Angeles (UCLA) Integrated Staging System (UISS) and the Stage Size Grade Necrosis (SSIGN) Risk Score (4–7). These models use clinical parameters and pathologic features for recurrence prognostication. However, with growing insight in the molecular biology of ccRCC, these prognostic models fall short in adequately predicting outcome of the individual patient, and molecular markers that reflect the biologic behavior of ccRCC might add to their prognostic value. A key event in the development of ccRCCs is loss of function of the von Hippel-Lindau (VHL) gene (8), although (epi)genetic alterations of the VHL gene do not appear to be directly associated with prognosis (9, 10). Whole-genome and -exome sequencing approaches identified novel mutations in genes encoding chromatin-modifying proteins such as PBRM1 (11) and BAP1 (12, 13), and mutations in BAP1 are associated with poorer prognosis (14, 15). The identification of these mutations in a subset of ccRCCs points to downstream epigenetic alterations in these tumors.

Promoter methylation–induced gene silencing is an important aspect of RCC biology (16, 17) and represents promising biomarkers for the early detection of cancer and for the prediction of prognosis and response to therapy. Developments in approaches to analyze the DNA methylome, based on microarray and next-generation sequencing (NGS) technology (18), have facilitated the identification of cancer genes inactivated by DNA methylation.

In the present study, we integrated methyl-binding domain (MBD) affinity–based massive parallel sequencing data and global transcript expression microarray data obtained by pharmacologic inhibition of DNA methylation, to provide a comprehensive analysis of the ccRCC DNA methylome and to identify prognostic promoter methylation markers, which were evaluated in 3 independent ccRCC patient series.

Study populations

DNA from formalin-fixed, paraffin-embedded primary tumor samples from 2 independent well-characterized ccRCC patient series was used. The first population was a hospital-based series (n = 150) derived from the archives of the Department of Histopathology (University Hospital of Leuven, Leuven, Belgium) and Department of Pathology (Maastricht University Medical Center, Maastricht, the Netherlands). The second population (n = 185) was obtained from the Netherlands Cohort Study on Diet and Cancer (NLCS; ref. 19), a prospective, population-based cohort study in which tumor tissue was collected from 51 pathology laboratories throughout the Netherlands. All patients in both the hospital- and population-based series underwent nephrectomy without any neoadjuvant therapy. Median follow-up in the hospital-based series was 64 months (range, 1–153 months) and in the population-based series was 79 months (range, 0–218 months). Tissue collection, DNA isolation, and patient characteristics (see also Tables 1 and 2) for both study populations have been described in detail elsewhere (20, 21). In addition, histologically normal renal tissue samples (formalin-fixed, paraffin-embedded) of 20 non-RCC patients were collected from the archives of the Department of Pathology (Maastricht University Medical Center).

Table 1.

GREM1, LAD1, NEURL, and NEFH promoter methylation associated with patient and tumor characteristics, including gene panels—hospital-based series

All ccRCCMethylated GREM1Methylated LAD1Methylated NEFHMethylated NEURLMethylated GREM1 and LAD1Methylated GREM1 and NEURL and LAD1Methylated GREM1 and LAD1 and NEURL and NEFH
Characteristicsn = 150 (100%)n = 28 (20%)n = 36 (27%)n = 48 (37%)n = 62 (44%)n = 14 (11%)n = 10 (8%)n = 7 (6%)
Sex 
 Male, n (%) 99 (66) 22 (79) 22 (61) 36 (75) 45 (79) 10 (71) 9 (90) 6 (86) 
 Female, n (%) 51 (34) 6 (21) 14 (39) 12 (25) 17 (21) 4 (29) 1 (10) 1 (14) 
P (U vs. M)  0.11 0.63 0.14 0.27 0.91 0.53 0.21 
Age, mean (±SD), y 60.1 (±12.1) 60.3 (±10.4) 61.3 (±13.1) 63.4 (±10.4) 61.6 (±12.3) 60.7 (±12.8) 61.9 (±14.1) 63.3 (±11.2) 
P (U vs. M)  0.96 0.44 0.01 0.20 0.79 0.34 0.29 
Tumor size, mean (±SD), cm 6.36 (±3.45) 7.72 (±3.17) 8.03 (±3.86) 6.88 (±3.00) 6.96 (±3.56) 8.11 (±3.44) 8.72 (±3.80) 7.92 (±3.72) 
P (U vs. M)  0.02 0.0007 0.30 0.04 0.002 0.002 0.007 
Fuhrman grade, n (%) 
 1 15 (10) 1 (4) 1 (3) 4 (9) 5 (8) 0 (0) 0 (0) 0 (0) 
 2 72 (48) 11 (39) 10 (28) 17 (35) 25 (40) 4 (29) 3 (30) 2 (28.5) 
 3 46 (31) 9 (32) 14 (39) 17 (35) 22 (36) 4 (28) 3 (30) 2 (28.5) 
 4 17 (11) 7 (25) 11 (30) 10 (21) 10 (16) 6 (43) 4 (40) 3 (43) 
P (U vs. M)  0.04 <0.001 0.10 0.19 0.002 0.021 0.04 
Tumor stage, n (%) 
 1 8 (5) 0 (0) 0 (0) 1 (2) 2 (3) 0 (0) 0 (0) 0 (0) 
 2 89 (60) 12 (43) 16 (45) 26 (54) 31 (50) 5 (36) 3 (30) 3 (43) 
 3 30 (20) 7 (25) 12 (33) 10 (21) 17 (28) 4 (28) 4 (40) 2 (28.5) 
 4 22 (15) 9 (32) 8 (22) 11 (23) 12 (19) 5 (36) 3 (30) 2 (28.5) 
P (U vs. M)  0.013 0.016 0.118 0.063 0.025 0.007 0.075 
ECP %, mean (±SD) 1.4 (±2.17) 3.0 (±2.97) 2.5 (±3.35) 1.9 (±2.93) 2.1 (±2.92) 4.0 (±4.21) 4.5 (±4.8) 4.4 (±5.84) 
P (U vs. M)  0.0001 0.0011 0.1295 0.0039 <0.0001 <0.0001 0.0001 
TCP %, mean (±SD) 6.7 (±9.25) 12.0 (±11.90) 10.6 (±12.91) 7.8 (±10.99) 9.9 (±12.43) 15.8 (±13.6) 16.4 (±15.1) 14.3 (±17.88) 
P (U vs. M)  0.001 0.002 0.372 0.0006 0.0001 <0.0001 0.0008 
MVD per mm2, mean (±SD) 197.5 (±91.0) 152.5 (±92.2) 152.8 (±81.9) 178.1 (±83.6) 180.0 (±84.9) 114.0 (±51.5) 93.3 (±26.4) 96.1 (±28.5) 
P (U vs. M)  0.001 0.0008 0.132 0.053 0.0003 <0.0001 0.0001 
OS, mean (±SD), mo 63.5 (±36.8) 53.8 (±41.1) 52.6 (±32.9) 51.9 (±29.4) 57.4 (±39.2) 51.3 (±43.4) 49.8 (±41.0) 45.2 (±35.1) 
P (U vs. M)  0.11 0.05 0.02 0.08 0.16 0.13 0.20 
MFS, mean (±SD), mo 51.8 (±41.0) 35.8 (±41.8) 35.0 (±33.3) 38.5 (±32.3) 40.9 (±39.9) 32.1 (±39.3) 25.3 (±29.0) 31.1 (±32.9) 
P (U vs. M)  0.02 0.007 0.01 0.004 0.018 0.007 0.019 
HR (sex/age)  2.3 (CI, 1.3–4.2) 2.3 (CI, 1.2–4.3) 1.5 (CI, 0.8–2.8) 2.5 (CI, 1.4–4.5) 3.5 (CI, 1.6–7.9) 3.6 (CI, 1.4–9.6) 3.6 (CI, 1.0–13.0) 
HR (multi)  0.9 (CI, 0.5–1.7) 0.8 (CI, 0.4–1.6) 0.6 (CI, 0.3–1.2) 1.4 (CI, 0.7–2.7) 0.9 (CI, 0.3–2.4) 0.7 (CI, 0.2–2.4) 0.3 (CI, 0.1–1.9) 
All ccRCCMethylated GREM1Methylated LAD1Methylated NEFHMethylated NEURLMethylated GREM1 and LAD1Methylated GREM1 and NEURL and LAD1Methylated GREM1 and LAD1 and NEURL and NEFH
Characteristicsn = 150 (100%)n = 28 (20%)n = 36 (27%)n = 48 (37%)n = 62 (44%)n = 14 (11%)n = 10 (8%)n = 7 (6%)
Sex 
 Male, n (%) 99 (66) 22 (79) 22 (61) 36 (75) 45 (79) 10 (71) 9 (90) 6 (86) 
 Female, n (%) 51 (34) 6 (21) 14 (39) 12 (25) 17 (21) 4 (29) 1 (10) 1 (14) 
P (U vs. M)  0.11 0.63 0.14 0.27 0.91 0.53 0.21 
Age, mean (±SD), y 60.1 (±12.1) 60.3 (±10.4) 61.3 (±13.1) 63.4 (±10.4) 61.6 (±12.3) 60.7 (±12.8) 61.9 (±14.1) 63.3 (±11.2) 
P (U vs. M)  0.96 0.44 0.01 0.20 0.79 0.34 0.29 
Tumor size, mean (±SD), cm 6.36 (±3.45) 7.72 (±3.17) 8.03 (±3.86) 6.88 (±3.00) 6.96 (±3.56) 8.11 (±3.44) 8.72 (±3.80) 7.92 (±3.72) 
P (U vs. M)  0.02 0.0007 0.30 0.04 0.002 0.002 0.007 
Fuhrman grade, n (%) 
 1 15 (10) 1 (4) 1 (3) 4 (9) 5 (8) 0 (0) 0 (0) 0 (0) 
 2 72 (48) 11 (39) 10 (28) 17 (35) 25 (40) 4 (29) 3 (30) 2 (28.5) 
 3 46 (31) 9 (32) 14 (39) 17 (35) 22 (36) 4 (28) 3 (30) 2 (28.5) 
 4 17 (11) 7 (25) 11 (30) 10 (21) 10 (16) 6 (43) 4 (40) 3 (43) 
P (U vs. M)  0.04 <0.001 0.10 0.19 0.002 0.021 0.04 
Tumor stage, n (%) 
 1 8 (5) 0 (0) 0 (0) 1 (2) 2 (3) 0 (0) 0 (0) 0 (0) 
 2 89 (60) 12 (43) 16 (45) 26 (54) 31 (50) 5 (36) 3 (30) 3 (43) 
 3 30 (20) 7 (25) 12 (33) 10 (21) 17 (28) 4 (28) 4 (40) 2 (28.5) 
 4 22 (15) 9 (32) 8 (22) 11 (23) 12 (19) 5 (36) 3 (30) 2 (28.5) 
P (U vs. M)  0.013 0.016 0.118 0.063 0.025 0.007 0.075 
ECP %, mean (±SD) 1.4 (±2.17) 3.0 (±2.97) 2.5 (±3.35) 1.9 (±2.93) 2.1 (±2.92) 4.0 (±4.21) 4.5 (±4.8) 4.4 (±5.84) 
P (U vs. M)  0.0001 0.0011 0.1295 0.0039 <0.0001 <0.0001 0.0001 
TCP %, mean (±SD) 6.7 (±9.25) 12.0 (±11.90) 10.6 (±12.91) 7.8 (±10.99) 9.9 (±12.43) 15.8 (±13.6) 16.4 (±15.1) 14.3 (±17.88) 
P (U vs. M)  0.001 0.002 0.372 0.0006 0.0001 <0.0001 0.0008 
MVD per mm2, mean (±SD) 197.5 (±91.0) 152.5 (±92.2) 152.8 (±81.9) 178.1 (±83.6) 180.0 (±84.9) 114.0 (±51.5) 93.3 (±26.4) 96.1 (±28.5) 
P (U vs. M)  0.001 0.0008 0.132 0.053 0.0003 <0.0001 0.0001 
OS, mean (±SD), mo 63.5 (±36.8) 53.8 (±41.1) 52.6 (±32.9) 51.9 (±29.4) 57.4 (±39.2) 51.3 (±43.4) 49.8 (±41.0) 45.2 (±35.1) 
P (U vs. M)  0.11 0.05 0.02 0.08 0.16 0.13 0.20 
MFS, mean (±SD), mo 51.8 (±41.0) 35.8 (±41.8) 35.0 (±33.3) 38.5 (±32.3) 40.9 (±39.9) 32.1 (±39.3) 25.3 (±29.0) 31.1 (±32.9) 
P (U vs. M)  0.02 0.007 0.01 0.004 0.018 0.007 0.019 
HR (sex/age)  2.3 (CI, 1.3–4.2) 2.3 (CI, 1.2–4.3) 1.5 (CI, 0.8–2.8) 2.5 (CI, 1.4–4.5) 3.5 (CI, 1.6–7.9) 3.6 (CI, 1.4–9.6) 3.6 (CI, 1.0–13.0) 
HR (multi)  0.9 (CI, 0.5–1.7) 0.8 (CI, 0.4–1.6) 0.6 (CI, 0.3–1.2) 1.4 (CI, 0.7–2.7) 0.9 (CI, 0.3–2.4) 0.7 (CI, 0.2–2.4) 0.3 (CI, 0.1–1.9) 

NOTE: Tumor stage based on International Union against Cancer (UICC) TNM Classification of Malignant Tumors, 4th, fully revised edition, 1987.

Abbreviations: HR (multi): HR (RCC-related death) corrected for age, sex, tumor size, stage, and grade; HR (sex/age): HR (RCC-related death) corrected for sex and age; MFS, metastasis-free survival; OS, overall survival; P, P value for comparison of methylated (M) and unmethylated (U) cases.

Table 2.

GREM1, LAD1, NEURL, and NEFH promoter methylation associated with patient and tumor characteristics, including gene panels—population-based series

All ccRCCMethylated GREM1Methylated LAD1Methylated NEFHMethylated NEURLMethylated GREM1 and LAD1Methylated NEURL and LAD1 and NEFHMethylated GREM1 and LAD1 and NEURL and NEFH
Characteristicsn = 185 (100%)n = 68 (40%)n = 60 (35%)n = 74 (52%)n = 88 (47%)n = 31 (20%)n = 26 (20%)n = 22 (18%)
Sex 
 Male, n (%) 108 (59) 44 (65) 42 (70) 50 (68) 46 (52) 22 (71) 17 (65) 15 (68) 
 Female, n (%) 74 (41) 24 (35) 18 (30) 24 (32) 42 (48) 9 (29) 9 (35) 7 (32) 
P (U vs. M)  0.32 0.07 0.17 0.05 0.06 0.11 0.16 
Age, mean (±SD), y 67.4 (±4.7 ) 67.2 (±5.2) 68.6 (±4.7) 68.0 (±4.4) 67.0 (±4.6) 67.4 (±5.3) 67.8 (±5.1) 67.2 (±4.5) 
P (U vs. M)  0.49 0.02 0.78 0.58 0.29 0.33 0.83 
Tumor size, mean (±SD), cm 7.9 (±3.9) 7.5 (±3.17) 7.9 (±3.40) 7.3 (±3.01) 6.9 (±3.30) 7.4 (±2.8) 8.4 (±3.3) 8.1 (±2.78) 
P (U vs. M)  0.16 0.01 0.87 0.08 0.01 0.02 0.08 
Fuhrman grade, n (%) 
 1 42 (23) 10 (15) 9 (15) 11 (15) 30 (34) 3 (10) 2 (8) 2 (9) 
 2 66 (36) 22 (32) 15 (25) 20 (27) 35 (40) 5 (16) 3 (11) 2 (9) 
 3 47 (26) 23 (34) 23 (38) 29 (39) 17 (19) 15 (48) 12 (46) 11 (50) 
 4 27 (15) 13 (19) 13 (22) 14 (19) 6 (7) 8 (26) 9 (35) 7 (32) 
P (U vs. M)  0.03 0.03 0.01 <0.001 0.007 0.003 0.007 
Tumor stage, n (%) 
 1 4 (2) 0 (0) 0 (0) 0 (0) 3 (3) 0 (0) 0 (0) 0 (0) 
 2 81 (45) 27 (40) 19 (33) 26 (36) 47 (54) 11 (35) 6 (23) 5 (23) 
 3 63 (35) 23 (34) 26 (45) 27 (37) 26 (30) 10 (32) 10 (38) 7 (32) 
 4 32 (18) 18 (26) 13 (22) 20 (27) 11 (13) 10 (32) 10 (38) 10 (45) 
P (U vs. M)  0.08 0.035 0.009 0.008 0.006 0.001 0.001 
OS, mean (±SD), mo 82.9 (±64.1) 63.8 (±60.7) 51.3 (±50.4) 60.5 (±57.3) 64.8 (±58.6) 52.6 (±52.2) 46.2 (±56.1) 47.1 (±56.3) 
P (U vs. M)  0.001 <0.0001 <0.0001 0.001 <0.0001 <0.0001 <0.0001 
HR (sex/age)  2.3 (CI, 1.5–3.5) 2.3 (CI, 1.5–3.6) 2.8 (CI, 1.7–4.6) 1.9 (CI, 1.2–2.9) 3.9 (CI, 2.2–7.1) 6.7 (CI, 2.8–15.9) 7.5 (CI, 2.7–21.2) 
HR (multi)  2.3 (CI, 1.5–3.6) 1.8 (CI, 1.1–3.1) 2.6 (CI, 1.5–4.6) 1.5 (CI, 1.0–2.4) 3.4 (CI, 1.7–6.6) 5.7 (CI, 2.0–15.7) 6.6 (CI, 2.1–20.5) 
All ccRCCMethylated GREM1Methylated LAD1Methylated NEFHMethylated NEURLMethylated GREM1 and LAD1Methylated NEURL and LAD1 and NEFHMethylated GREM1 and LAD1 and NEURL and NEFH
Characteristicsn = 185 (100%)n = 68 (40%)n = 60 (35%)n = 74 (52%)n = 88 (47%)n = 31 (20%)n = 26 (20%)n = 22 (18%)
Sex 
 Male, n (%) 108 (59) 44 (65) 42 (70) 50 (68) 46 (52) 22 (71) 17 (65) 15 (68) 
 Female, n (%) 74 (41) 24 (35) 18 (30) 24 (32) 42 (48) 9 (29) 9 (35) 7 (32) 
P (U vs. M)  0.32 0.07 0.17 0.05 0.06 0.11 0.16 
Age, mean (±SD), y 67.4 (±4.7 ) 67.2 (±5.2) 68.6 (±4.7) 68.0 (±4.4) 67.0 (±4.6) 67.4 (±5.3) 67.8 (±5.1) 67.2 (±4.5) 
P (U vs. M)  0.49 0.02 0.78 0.58 0.29 0.33 0.83 
Tumor size, mean (±SD), cm 7.9 (±3.9) 7.5 (±3.17) 7.9 (±3.40) 7.3 (±3.01) 6.9 (±3.30) 7.4 (±2.8) 8.4 (±3.3) 8.1 (±2.78) 
P (U vs. M)  0.16 0.01 0.87 0.08 0.01 0.02 0.08 
Fuhrman grade, n (%) 
 1 42 (23) 10 (15) 9 (15) 11 (15) 30 (34) 3 (10) 2 (8) 2 (9) 
 2 66 (36) 22 (32) 15 (25) 20 (27) 35 (40) 5 (16) 3 (11) 2 (9) 
 3 47 (26) 23 (34) 23 (38) 29 (39) 17 (19) 15 (48) 12 (46) 11 (50) 
 4 27 (15) 13 (19) 13 (22) 14 (19) 6 (7) 8 (26) 9 (35) 7 (32) 
P (U vs. M)  0.03 0.03 0.01 <0.001 0.007 0.003 0.007 
Tumor stage, n (%) 
 1 4 (2) 0 (0) 0 (0) 0 (0) 3 (3) 0 (0) 0 (0) 0 (0) 
 2 81 (45) 27 (40) 19 (33) 26 (36) 47 (54) 11 (35) 6 (23) 5 (23) 
 3 63 (35) 23 (34) 26 (45) 27 (37) 26 (30) 10 (32) 10 (38) 7 (32) 
 4 32 (18) 18 (26) 13 (22) 20 (27) 11 (13) 10 (32) 10 (38) 10 (45) 
P (U vs. M)  0.08 0.035 0.009 0.008 0.006 0.001 0.001 
OS, mean (±SD), mo 82.9 (±64.1) 63.8 (±60.7) 51.3 (±50.4) 60.5 (±57.3) 64.8 (±58.6) 52.6 (±52.2) 46.2 (±56.1) 47.1 (±56.3) 
P (U vs. M)  0.001 <0.0001 <0.0001 0.001 <0.0001 <0.0001 <0.0001 
HR (sex/age)  2.3 (CI, 1.5–3.5) 2.3 (CI, 1.5–3.6) 2.8 (CI, 1.7–4.6) 1.9 (CI, 1.2–2.9) 3.9 (CI, 2.2–7.1) 6.7 (CI, 2.8–15.9) 7.5 (CI, 2.7–21.2) 
HR (multi)  2.3 (CI, 1.5–3.6) 1.8 (CI, 1.1–3.1) 2.6 (CI, 1.5–4.6) 1.5 (CI, 1.0–2.4) 3.4 (CI, 1.7–6.6) 5.7 (CI, 2.0–15.7) 6.6 (CI, 2.1–20.5) 

NOTE: Tumor stage based on International Union against Cancer (UICC) TNM Classification of Malignant Tumors, 4th, fully revised edition, 1987.

Abbreviations: HR (multi): HR (RCC-related death) corrected for age, sex, tumor size, stage, and grade; HR (sex/age): HR (RCC-related death) corrected for sex and age; OS, overall survival; P, P value for comparison of methylated (M) and unmethylated (U) cases.

This study was approved by the Medical Ethical Committee of the Maastricht University Medical Center.

Cell culture and drug treatment

Four ccRCC cell lines (SKRC1, SKRC10, SKRC52, SKRC59), kindly provided by Dr. E. Oosterwijk, Nijmegen Center for Molecular Life Sciences (NCMLS), The Netherlands, were cultured in RPMI-1640 media (Invitrogen) supplemented with 10% heat-inactivated FBS. The human kidney 2 cell line (HK-2) was purchased from ATCC and cultured in keratinocyte serum-free medium (K-SFM, GIBCO) supplemented with bovine pituitary extract (BPE) and human recombinant EGF. SKRC1, SKRC10, SKRC52, and SKRC59 cells were treated with 5-aza-2′-deoxycytidine (DAC) or trichostatin A (TSA) as previously described (22, 23). Approximately 10% confluent RCC cells were cultured in RPMI-1640 media containing 10% FBS with 5 μmol/L DAC (Sigma; stock solution: 1 mmol/L in PBS) for 96 hours, replacing media and DAC every 24 hours. Cell treatment with 300 nmol/L TSA (Sigma; stock solution: 5 mmol/L dissolved in ethanol) was performed for 18 hours, starting at approximately 30% cell confluence. Mock-treated cells were grown in parallel with the DAC treatment by adding equal volumes of PBS without drugs.

RNA expression microarray analyses

Microarray expression analyses were performed on mock-, DAC-, or TSA-treated cells (SKRC1, SKRC10, SKRC52, SKRC59), as previously described (22, 24). Total RNA was isolated using the RNeasy Kit (Qiagen) and quantified using the NanoDrop ND-100 followed by quality assessment with the 2100 Bioanalyzer [RIN = 10 (n = 10), 9.9 (n = 1), or 8.9 (n = 1), Agilent Technologies]. Sample amplification and labeling procedures were carried out using the Low RNA Input Fluorescent Linear Amplification Kit (Agilent Technologies). The labeled cRNA was purified using the RNeasy Mini Kit (Qiagen) and quantified. Cy3- or Cy5-labeled samples were mixed with control targets (Agilent Technologies), assembled on the Agilent Whole Human Genome 4 × 44 k microarray, hybridized, and processed according to the Agilent microarray protocol. Scanning was performed with the Agilent G2565BA microarray scanner. Array data were analyzed using R (version 2.10.0) and BioConductor, using the limma package (version 3.2.1; ref. 25). Median Cy5 and Cy3 signals read from the raw data and M-values, that is, log2(Cy5/Cy3), were loess-normalized. Basal expression was approximated using the single-channel Cy3 values, containing the mock-treated samples. Microarray data are retrievable from the Gene Expression Omnibus (GEO) database at the National Center for Biotechnology Information (NCBI; accession number GSE33916).

MBD affinity–based next-generation sequencing

DNA of SKRC1, SKRC10, SKRC52, SKRC59, and HK-2 cells were isolated using the Puregene DNA Purification Kit (BIOzym) according to the manufacturer's instructions. Cell line gDNA was fragmented using a COVARIS S2 system with AFA fiber micro-tubes to obtain fragments with an average length of 200 bp. Methylated DNA was captured by pull-down using the MethylCollector Kit (cancer cell lines, Active Motif) or MethylCap with High-Salt elution (HK-2, Diagenode), both on the basis of autologous methyl-binding domains, according to the manufacturer's protocols. Subsequently, fragments were sequenced using the Illumina Genome Analyzer II. The concentration of the fragmented and captured DNA was determined on a Fluostar Optima plate reader with the Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen P7589) 480/520 nm. For all samples together, the paired-end 40-bp sequence reads were mapped using BOWTIE (26) on the human reference genome (NCBI build 37.3). Coverage values were summarized using the Map of the Human Methylome developed by BioBix, which consists of putatively independently methylated regions (“methylation cores”) throughout the genome (http://www.biobix.be/map-of-the-human-methylome/mhm-version-2/). For each sample, and each methylation core, the maximum read coverage was used for further analysis. Subsequently, we focused on the broad promoter region for each gene (−2,000 to +500 bp). A Poisson background model was used to identify significantly methylated regions (null hypothesis: all signal is randomly distributed as background, α = 0.01), with λ estimated as the sum of the coverage values over all promoter methylation cores divided by the number of promoter methylation cores. This approach takes into account coverage differences between samples, although generally low coverage will result in low sensitivity. For absence of methylation in HK-2, we allowed a maximum coverage value of 2, corresponding with P ≤ 0.067. We have deposited the MBD affinity NGS data in the GEO database at NCBI (GSE33954).

Locus-specific DNA methylation and gene expression analysis

Promoter DNA methylation was determined by chemical modification of gDNA with sodium bisulfite and subsequent nested multiplex methylation-specific PCR (MSP) as described in detail elsewhere (27). Five hundred nanograms of DNA was modified by sodium bisulfite using the EZ DNA Methylation Kit (Zymo Research). All PCRs were performed with control samples containing only unmethylated DNA [DNA from normal lymphocytes or DNA from human umbilical vein endothelial cells (HUVEC)], only methylated DNA (normal lymphocyte DNA treated in vitro with SssI methyltransferase (New England Biolabs)], and a blank control without DNA. Primer sequences are provided in Supplementary Table S1. RNA of mock-, DAC-, and TSA-treated cells was isolated using the RNeasy Mini Kit (Qiagen) and treated with RNase-free DNase (Qiagen). Synthesis of cDNA was performed using the Iscript cDNA Synthesis Kit (Bio-Rad). Quantitative reverse transcription PCR was performed as described previously (28) using SYBR GreenPCR master mix (Bio-Rad). Cyclophilin A was used as a reference gene for normalization. Primer sequences for qRT-PCR are provided in Supplementary Table S2.

Clinical data analyses

Cause-specific survival was defined as the time from cancer diagnosis until renal cancer–related death or end of follow-up. Differences in clinicopathologic and angiogenesis characteristics between ccRCCs with and without candidate promoter methylation were evaluated with the Student t and Pearson χ2 tests, where appropriate. Kaplan–Meier analyses and log-rank tests were used to estimate the overall influence of candidate gene promoter methylation on cause-specific survival. HRs and corresponding 95% confidence intervals (CI) were assessed using Cox proportional hazards models. Assuming a 2-sided α of 0.05, a power of 0.80, and a median survival of 75 months in the unmethylated cases, a power analysis showed that the minimal detectable HR (MDHR) was 2.91 (15% methylated), 2.38 (25% methylated), and 2.02 (40% methylated) in the hospital-based series and 2.18 (15% methylated), 1.87 (25% methylated), and 1.72 (40% methylated) in the population-based series. Known prognostic factors for renal cancer were considered possible confounders if they influenced the crude HR by more than 10%. Possible confounders that were included in the model for both series were sex, age at diagnosis, cancer stage, tumor size, and nuclear grade. The proportional hazard assumption was tested using the Schoenfeld residuals and the log(−log) hazards plots. HarrellC statistic and the Akaike Information Criterion (AIC) were used to assess the predictive capacity and fit of the models. As a model with a 100% sensitivity and specificity would yield a HarrellC statistic of 1.00, the model with the highest HarrellC statistic was regarded as the best model. In general, a lower AIC is considered a better model. All analyses were performed with the statistical package STATA 11.0.

The findings were additionally validated in an independent set of ccRCC samples from The Cancer Genome Atlas (TGCA). Matched methylation and clinical follow-up data were available for 264 (Infinium HumanMethylation 450 K) or 467 unique patients (Infinium HumanMethylation 27 K and 450 K), depending on the available TCGA methylation analysis platform. All probes within a CpG island and within −1,000 to +500 bp of an annotated GREM1, LAD1, NEFH, or NEURL promoter were inspected. Histograms were generated for each one of the 28 selected probes to assess the distribution of methylation β-values. Using the bimodal distribution of the histograms, binary β-value cutoffs were derived per probe, calling samples either methylated or unmethylated (see Supplementary Fig. S1B). Cox proportional hazards models were fit for each probe. HRs were adjusted for age at diagnosis and sex and in the multivariable models also for pathologic stage and nuclear grade (tumor size was not available in the TCGA data). Because of the low number of patients with grade 1 ccRCC (n = 6), grade 1 and 2 ccRCCs were clustered in analyses. All TCGA analyses were executed in R.

All reported P values are 2-sided and P ≤ 0.05 was considered statistically significant.

Comprehensive characterization of the ccRCC methylome

To comprehensively determine the ccRCC methylome, we used an integrative approach combining data obtained by MBD affinity NGS (29) and expression microarray data of DAC-treated, globally demethylated, ccRCC cell lines (22).

Changes in gene expression following DAC or TSA treatment of 4 ccRCC cell lines are shown in Fig. 1A. As previously described, the characteristic spikes identify a zone in which gene expression did not increase upon treatment with TSA (TSA-negative; <1.4- and >0.71-fold), showed increased expression >1.4-fold after DAC treatment (DAC-positive), with no expression in the mock-treated cells (22, 23). In total, 1583 unique genes were identified by this approach, ranging between 508 and 731 genes per cell line. The number of cell lines expressing these genes is shown in Fig. 1B.

Figure 1.

Characteristics of gene discovery approaches in ccRCC cell lines. Gene expression changes for SKRC1, SKRC10, SKRC52, and SKRC59 cells treated with TSA (x-axis) or DAC (y-axis) are plotted by log-fold change, and each dot represents an individual gene (probe) on the array. Overlapping genes with MBD binding are indicated in red (without overlap are in black; A). Distribution and shared candidate methylated genes in the ccRCC cell lines are shown after DAC/TSA treatment (B) or MBD affinity–based NGS (C). Overlap in gene expression changes or MBD binding among 2, 3, or 4 cell lines is indicated. Integration of the 2 approaches to induce the likelihood to identify biological relevant candidate genes is depicted per cell line (D). P values indicate a significant overlap between the 2 approaches in all cell lines. M, methylated; U, unmethylated.

Figure 1.

Characteristics of gene discovery approaches in ccRCC cell lines. Gene expression changes for SKRC1, SKRC10, SKRC52, and SKRC59 cells treated with TSA (x-axis) or DAC (y-axis) are plotted by log-fold change, and each dot represents an individual gene (probe) on the array. Overlapping genes with MBD binding are indicated in red (without overlap are in black; A). Distribution and shared candidate methylated genes in the ccRCC cell lines are shown after DAC/TSA treatment (B) or MBD affinity–based NGS (C). Overlap in gene expression changes or MBD binding among 2, 3, or 4 cell lines is indicated. Integration of the 2 approaches to induce the likelihood to identify biological relevant candidate genes is depicted per cell line (D). P values indicate a significant overlap between the 2 approaches in all cell lines. M, methylated; U, unmethylated.

Close modal

Differentially methylated DNA sequences were identified by massive parallel sequencing of MBD-binding enriched DNA. MBD affinity NGS of the normal renal epithelial cell line HK-2 was performed to identify tumor-specific methylation. In total, 7,829 unique genes were found to have promoter methylation, with 2802 to 5,695 unique genes per cell line. The distribution of differentially methylated genes among the cell lines is depicted in Fig. 1C. The advantages of this integrated approach are demonstrated in Fig. 1D, which also takes into account the lower genomic coverage of expression arrays. Although we expected that treatment with TSA would exclude genes repressed by changes other than DNA methylation (22, 23), we found that there were still many genes with increased expression following treatment with DAC that did not contain promoter methylation. Interestingly, only 5.8% to 8.6% of the genes with methylated promoter regions were clearly re-expressed after DAC treatment.

To further identify candidate genes with functional methylation, we applied the following criteria: MBD binding in the broad promoter region (−2,000 to +500 bp) in conjunction with re-expression after demethylation by DAC in at least 3 of 4 ccRCC cell lines, but no methylation in the normal kidney cell line HK-2. To increase the possibility of the identification of (renal cell) cancer-specific genes, we focused on candidates involved in DNA repair, angiogenesis and hypoxia, carbohydrate metabolism, the TGFβ or Notch pathway based on Gene Ontology (GO) terms, and genes known to be mutated in RCC (11, 16) or frequently mutated in other tumors (ref. 30; Fig. 2). This resulted in 43 candidate genes which are, to the best of our knowledge, not previously reported to show promoter CpG island methylation in RCC, except for CST6 (31, 32), GREM1 (21, 31), NEFH (33), LAD1 (32), QPCT (17), and BMP2 (ref. 34; for a complete list, see Supplementary Table S3).

Figure 2.

Pipeline to select candidate methylation markers in ccRCC cell lines. Schematic view of developed pipeline applied to the identification of biologically relevant candidate genes. Genes identified in ccRCC cell lines as specifically upregulated after DAC treatment and enriched by MBD protein(s), but not in HK-2 cells, were considered as candidate promoter methylation markers. Those genes that met the downstream indicated requirements were selected for further analyses.

Figure 2.

Pipeline to select candidate methylation markers in ccRCC cell lines. Schematic view of developed pipeline applied to the identification of biologically relevant candidate genes. Genes identified in ccRCC cell lines as specifically upregulated after DAC treatment and enriched by MBD protein(s), but not in HK-2 cells, were considered as candidate promoter methylation markers. Those genes that met the downstream indicated requirements were selected for further analyses.

Close modal

Validation of promoter methylation and gene expression in ccRCC cell lines

To validate the used approach, we examined the methylation and mRNA expression status in ccRCC cell lines by MSP and qRT-PCR, respectively. We designed MSP primers for all 43 candidate genes at promoter regions which were identified by the MBD affinity NGS approach in RCC cell lines, but not in normal kidney cells. Validation of the expression microarray data resulted in a sensitivity and specificity of 80% (Supplementary Table S4). Validation of the MBD-NGS data resulted in a sensitivity of 90% and a specificity of 63% (Supplementary Table S4). Similarly, upregulation of gene expression after demethylation was validated by qRT-PCR in 74% using a subset (n = 22) of genes (Supplementary Fig. S2).

Promoter methylation in primary ccRCC

To determine whether promoter methylation of the candidate genes was also present in primary ccRCC, we first evaluated the 43 candidate genes on a pilot of 20 primary ccRCC samples and 20 normal kidney samples by MSP (Supplementary Table S3). Genes with a promoter methylation frequency of less than 20% in ccRCC (n = 15) or more than 10% in normal kidney (n = 14) were not further analyzed. This resulted in 14 genes (CSPG4, GREM1, KRT7, LAD1, NEFH, NEURL, QPCT, RGMA, SORL1, CD109, CHRDL2, FST, PLXDC1, RASGEF1A) that were further validated in a hospital-based series (n = 150) of primary ccRCCs. Of these, 5 genes (CD109, CHRDL2, FST, PLXDC1, and RASGEF1A) were additionally excluded as none or only few additional methylated cases were detected. Taken together, this resulted in 9 candidate genes (CSPG4, GREM1, KRT7, LAD1, NEFH, NEURL, QPCT, RGMA, SORL1) with promoter methylation frequencies ranging from 16% to 58% in primary ccRCC.

Promoter methylation of GREM1, LAD1, NEFH, and NEURL has prognostic value in primary ccRCC

The association between promoter methylation of the 9 candidate genes and tumor/patient characteristics and clinical outcome in the hospital-based series (n = 150) is shown in Table 1, Supplementary Table S5, Supplementary Fig. S3A–S3E, and Supplementary Fig. S4A, S4C, S4E, S4G.

The mean ccRCC-specific survival in the hospital-based series was 5.3 years. For 15% of the patients, a distant metastasis was reported at time of diagnosis (TNM stage IV). Distribution of age and sex did not differ by candidate gene methylation, except patients with NEFH methylation were significantly older than those without (P = 0.01; Table 1). In total, 51 patients died because of ccRCC during follow-up; 31 patients with localized disease (stage I–III) and 20 patients with metastasized disease (stage IV).

Promoter methylation of GREM1, LAD1, NEFH, and NEURL showed statistically significant associations with advanced disease and/or poor prognosis (Table 1). In addition to GREM1 (as we published previously; ref. 21), promoter methylation of LAD1 and NEURL showed an association with poor survival (log-rank P = 0.013 and <0.001, respectively; Supplementary Fig. S4C and S4G), tumor size, grade, and stage, tumor cell proliferation (TCP), endothelial cell proliferation (ECP), and microvessel density (MVD; Table 1). For NEFH, a trend toward poorer survival was seen (log-rank P = 0.1732; Supplementary Fig. S4E).

As shown in Fig. 3, the association between poor survival and promoter methylation of GREM1, LAD1, NEFH, and NEURL was also seen in the age- and sex-adjusted Cox proportional hazard analyses [HRGREM1, 2.30 (95% CI, 1.26–4.19); HRLAD1, 2.30 (95% CI, 1.24–4.26); HRNEFH, 1.50 (95% CI, 0.80–2.82); HRNEURL, 2.50 (95% CI, 1.38–4.53)], although again, this was not statistically significant for NEFH.

Figure 3.

HRs for RCC-related death for single markers and prognostic marker panels in the hospital-, population-, and TCGA-based ccRCC series. Forest plot depicting the increasing risk of RCC-related death by HRs (corrected for sex and age) for GREM1, LAD1, NEFH, and NEURL and the 2 and 3 best predictive marker panels according to Harrell C statistic as well as the 4-marker panel. Two- and 3-marker combinations hospital-based series: GREM1 and LAD1, and GREM1, LAD1, and NEURL, respectively; population-based series: GREM1 and LAD1, and LAD1, NEFH, and NEURL, respectively; and TCGA series: GREM1 and LAD1, and GREM1, LAD1, and NEURL, respectively.

Figure 3.

HRs for RCC-related death for single markers and prognostic marker panels in the hospital-, population-, and TCGA-based ccRCC series. Forest plot depicting the increasing risk of RCC-related death by HRs (corrected for sex and age) for GREM1, LAD1, NEFH, and NEURL and the 2 and 3 best predictive marker panels according to Harrell C statistic as well as the 4-marker panel. Two- and 3-marker combinations hospital-based series: GREM1 and LAD1, and GREM1, LAD1, and NEURL, respectively; population-based series: GREM1 and LAD1, and LAD1, NEFH, and NEURL, respectively; and TCGA series: GREM1 and LAD1, and GREM1, LAD1, and NEURL, respectively.

Close modal

The 2 marker combination of GREM1 and LAD1 [age-/sex-adjusted HR both genes methylated, 3.54 (95% CI, 1.58–7.94), HarrellC: 0.6415] and the 3 marker combination of GREM1, LAD1, and NEURL [age-/sex-adjusted HR 3 genes methylated, 3.61 (95% CI, 1.36–9.58), HarrellC: 0.6436] showed the best predictive capacity according to the Harrell C statistic. The model combining all 4 markers showed a slightly elevated age-/sex-adjusted HR of 3.64 (95% CI, 1.02–13.00) when compared with the best 3-marker panel, but the HarrellC statistic does not indicate a higher predictive capacity of this model (Harrell C: 0.6212), although the fit of the 4-marker model is better according to the AIC as compared with the 3-marker model (AIC, 321 and 346, respectively). For all 2- and 3-marker panel combinations, see Supplementary Table S6.

Figure 4 shows the overall cause-specific survival of patients with ccRCC using the model with the best predictive capacity according to the HarrellC statistic for the 2- (A) and 3- (C) marker panels as compared with the 4-marker (E) panel. In multivariate analyses adjusted for grade, stage, and tumor size, in addition to age and sex, however, the association between promoter methylation of GREM1, LAD1, NEFH, and/or NEURL and cause-specific survival was not statistically significant (Table 1).

Figure 4.

Cause-specific overall survival curves for marker panels in hospital- and population-based ccRCC series. Kaplan–Meier curves of the best 2-, 3-, and 4-marker panels in the hospital-based series: GREM1 and LAD1 (A); GREM1, LAD1, and NEURL (C); and GREM1, LAD1, NEFH, and NEURL (E), respectively, and population-based series: GREM1 and LAD1 (B); LAD1, NEFH, and NEURL (D); and GREM1, LAD1, NEFH, and NEURL (F), respectively. M, methylated; Number at risk, number of patients surviving at each time point; U, unmethylated.

Figure 4.

Cause-specific overall survival curves for marker panels in hospital- and population-based ccRCC series. Kaplan–Meier curves of the best 2-, 3-, and 4-marker panels in the hospital-based series: GREM1 and LAD1 (A); GREM1, LAD1, and NEURL (C); and GREM1, LAD1, NEFH, and NEURL (E), respectively, and population-based series: GREM1 and LAD1 (B); LAD1, NEFH, and NEURL (D); and GREM1, LAD1, NEFH, and NEURL (F), respectively. M, methylated; Number at risk, number of patients surviving at each time point; U, unmethylated.

Close modal

Promoter methylation of GREM1, LAD1, NEFH, and NEURL is an independent predictor of ccRCC survival

The results observed in the first series were validated in a second independent population-based series of ccRCCs (n = 185). Mean cause-specific survival in the population-based series was higher than in the hospital-based series (5.3 and 6.9 years) and 18% of the patients had a distant metastasis (stage IV). Age and sex did not differ according to methylation status of the candidate genes, except patients with LAD1 methylation were significantly older than those without (P = 0.02; Table 2). During follow-up, 96 patients in the population-based series died because of ccRCC; 68 patients with localized disease (stage I–III) and 28 patients with metastasized disease (stage IV).

Table 2, Figure 3, and Supplementary Fig. S4B, S4D, S4F, S4H show that, as reported previously for GREM1, promoter methylation of the single genes LAD1, NEURL, and NEFH was associated with increased malignancy and poor prognosis [age-/sex-adjusted HRGREM1, 2.30 (95% CI, 1.51–3.51); HRLAD1, 2.30 (95% CI, 1.47–3.60); HRNEFH, 2.80 (95% CI, 1.71–4.58); HRNEURL, 1.90 (95% CI, 1.24–2.91)] independent of other prognostic factors [multivariate HRGREM1, 2.30 (95% CI, 1.45–3.61); HRLAD1, 1.80 (95% CI, 1.10–3.06); HRNEFH, 2.62 (95% CI, 1.50–4.59); HRNEURL, 1.50 (95% CI, 0.96–2.42)], although for NEURL statistical significance was only borderline. The best 2-marker panel according to HarrellC statistic was the panel of GREM1 and LAD1 [age-/sex-adjusted HR both genes methylated 3.91 (95% CI, 2.16–7.08), HarrellC: 0.6837]. The best 3-marker panel however appeared to be the combination of LAD1, NEFH, and NEURL [age-/sex-adjusted HR three genes methylated 6.69 (95% CI, 2.82–15.86), HarrellC: 0.7108], although the second best predictive model, GREM1, LAD1, and NEFH [age-/sex-adjusted HR three genes methylated 5.76 (95% CI, 2.61–12.74)] only had a slightly lower HarrellC of 0.7073 but also a lower AIC of 551 as compared with 585 for the NEURL, LAD1, and NEFH combination. For all 2- and 3-marker panel combinations, see Supplementary Table S7.

Again, the 4-marker panel appeared to be superior in predicting prognosis of patients with ccRCC [HarrellC: 0.7129, AIC: 540, log-rank P < 0.0001, age-/sex-adjusted HR 4 genes methylated 7.54 (95% CI, 2.68–21.19)]. Figure 4 shows the overall cause-specific survival of patients with ccRCC using the model with the best predictive capacity according to the HarrellC statistic for the 2- (B) and 3- (D) marker panels as compared with the 4-marker (F) panel. The 5-year survival rate drastically decreased according to the presence of promoter methylation: from 100% in cases without methylation to 77%, 50%–52%, and 35% in cases with 1-, 2- or 3-, or 4-gene promoters methylated, respectively.

In contrast with the hospital-based population, the association with poor prognosis in the population-based series remained statistically significant when adjusting for other prognostic factors (tumor size, stage, and grade); multivariate HR best 2-marker panel 3.36 (95% CI, 1.71–6.60); multivariate HR best 3-marker panel 5.65 (95% CI, 2.03–15.72); and multivariate HR 4-marker panel 6.56 (95% CI, 2.09–20.52). Thus, patients with promoter methylation of the 4-marker panel represent a clinically distinct group with an up to 6.6-fold increased risk of cancer-related death, independent of known prognostic factors, and a 40% chance of dying from RCC within the first year, as compared with patients without methylation of these genes.

In line with the results from the hospital- and population-based series, methylation status of the 4 genes was associated with patient outcome in the TCGA series. The best probe for each one of the 4 genes was identified on the basis of HarrellC statistic [age-/sex-adjusted HRGREM1(cg08495115) 2.62 (95% CI, 1.73–3.95), HarrellC: 0.660; HRLAD1(cg23291161) 2.40 (95% CI, 1.60–3.60), HarrellC: 0.668; HRNEURL(cg03337393) 1.99 (95% CI, 1.31–3.02), HarrellC: 0.647; and HRNEFH(cg18129621)1.68 (95% CI, 1.10–2.56), HarrellC: 0.630; Fig. 3 and Supplementary Fig. S1A]. Using the same probes, the best 2-marker panel was the combination of GREM1 and LAD1 [age-/sex-adjusted HR 3.07 (95% CI, 1.93–4.90), HarrellC: 0.675, AIC 945] and the best 3-marker panel was the combination of GREM1, LAD1, and NEURL [age-/sex-adjusted HR 3.31 (95% CI, 1.97–5.57), HarrellC: 0.669, AIC 946]. The 4-marker panel appeared to be the best model; although it had a comparable HarrellC and the same fit according to AIC as the best 3-marker panel, it reached a higher HR [age-/sex-adjusted HR 3.60 (95% CI, 2.02–6.40), HarrellC: 0.667, AIC 946]. Adjusted for stage and nuclear grade in addition to age and sex, the 4-marker panel only reached borderline statistical significance but had a better HarrellC and better fit [multivariate HR 1.78 (95% CI, 0.96–3.32), P = 0.068, HarrellC: 0.755, AIC 895].

In the present study, we integrated global transcript expression microarray data of pharmacologically demethylated ccRCC cell lines and MBD affinity NGS data of the same samples to explore the ccRCC methylome and identify prognostic methylation markers. Although the stringent selection criteria used to identify candidate genes implies exclusion of other potentially interesting methylation markers, our strategy resulted in the identification of 9 genes showing tumor-specific promoter methylation in a substantial proportion (16%–58%) of tumors from the first series. Furthermore, using this approach, we identified genes previously reported to be methylated in RCC or even to be associated with patient prognosis, demonstrating the relevance of the strict selection criteria we applied.

We showed that promoter methylation of GREM1, LAD1, NEFH, and NEURL was individually associated with increased malignancy and/or poor survival in 2 independent series (n = 150 and n = 185) of primary ccRCC samples. The combination of these 4 genes in a promoter methylation marker panel strongly predicted ccRCC survival independent from age and sex both in the hospital- and population-based series. However, statistical significance was lost in multivariate analyses within the hospital-based population, whereas multivariate significance was observed for the marker panel in the population-based series. Importantly, the incidence of methylation of all 4 markers was more abundant in the population-based series (18%) as compared with the hospital-bases series (6%), which may easily have influenced the power to find significant results for the 4-marker panel in the hospital-based series. Furthermore, the hospital-based series has been retrospectively collected from 2 pathology archives and could represent a more selected group of patients than the population-based series which has been prospectively collected from a representative population of 120,852 men and women nationwide. The prognostic value of the 4-marker panel was further validated in a third, independent series of patients with ccRCC derived from the TCGA database.

In the past decade, there has been a surge of studies reporting on potential prognostic methylation markers in RCC; however, the reproducibility of the findings is often low, leading to failure during further validation attempts. Furthermore, it is becoming more evident that the choice of primer/probe location is important, as the level of methylation can vary, with also methylation region–specific effects on patient survival (21, 35). In this study, we showed that promoter methylation of GREM1, NEURL, LAD1, and NEFH predicts prognosis of patients with ccRCC in different patient series analyzed by 2 different techniques. Further inspection of the location of the probes that were significantly associated with patient survival in the TCGA data showed that they were located in close proximity to the MSP region, further illustrating the accuracy and robustness of these 4 methylation markers.

The biologic function of the identified genes increases our knowledge on ccRCC development and progression. For GREM1 and NEURL, a known role in embryogenesis and morphogenesis of the kidney is described (36–38). GREM1 inhibits the ligands BMP 2, 4, and 7, hindering downstream TGFβ signaling. However, GREM1 also has a proangiogenic function by directly binding to VEGFR2 in a BMP-independent manner (39). In many cancer types, increased angiogenesis is associated with poorer prognosis. In contrast, we have previously shown that silencing of GREM1 by promoter CpG island methylation is associated with lower microvessel density (MVD), but with higher tumor grade and shorter survival (21). It has been demonstrated, however, that renal tumors with low MVD have significantly lower expression of factors that are important for stabilization and maturation of the vasculature, such as placental growth factor (PLGF1), as compared with tumors with high MVD (40). Therefore, inactivation of GREM1 by promoter methylation could hinder maturation of the tumor vasculature, whereby tumor cells can more easily escape into circulation and metastasize, accounting for the poorer prognosis. In line with this hypothesis, Chen and colleagues also found that increased expression of GREM1 was associated with increased MVD and increased progression-free survival in pancreatic neuroendocrine tumors (41).

Neuralized (NEURL) is a highly conserved E3 ubiquitin ligase, which acts upon Notch ligands to regulate Notch pathway signaling (42). Promoter methylation and subsequent downregulation of the Notch pathway member NEURL has been reported in colorectal cancer in which overexpression caused reduced colony formation in vitro (22).

The relatively unknown LAD1 protein encoded by the Ladinin 1 gene is an anchoring filament that is a component of basement membranes. It may contribute to the stability of the association of the epithelial layers with the underlying mesenchyme (43). Recently, promoter methylation of LAD1 has been described in RCC and was associated with decreased survival of patients treated with antiangiogenic therapies (32).

The NEFH encoding neurofilament heavy chain is known to be one of the major components of the neuronal cytoskeleton neurofilaments (44), and promoter methylation of this gene has been described in several cancer types (45, 46). Recently, in line with our results, promoter CpG island methylation of NEFH has been reported in RCC and was also associated with poor progression-free survival (33).

LAD1 and NEFH could thus possibly be involved in cell structure and stability, but their role in (RCC) biology needs to be further elucidated.

This study has some important shortcomings. Specifically, we were not able to study intratumor heterogeneity (ITH). It is known that great genetic ITH exists in RCC which can lead to drug resistance and treatment failure and to difficulties in validating cancer biomarkers due to sampling bias (47, 48). However, Stewart and colleagues showed that despite great genetic ITH, there is relative homogeneity of the methylome, suggesting that the methylome is rather stable during tumor progression and that exploiting DNA methylation alterations as potential cancer biomarkers could therefore be very attractive (49). Furthermore, assuming that methylome heterogeneity is actually present, we believe that identification of DNA methylation alterations that are significantly associated with patient survival in samples that are intrinsically heterogeneous, only further endorses the biologic importance of these changes.

Second, because of lack of information on tumor necrosis and patient performance status in our patient series, we were not able to evaluate our model in comparison to other prognostic models such as the UISS and the SSIGN score. Finally, because of the (prospective-)retrospective study design, the highest level of evidence cannot be reached. Further validation in large, prospectively collected patient series in which also the additional prognostic value of the panel to standard clinicopathologic parameters and the recently identified prognostic mutations in BAP1 is assessed using multivariable prediction models is needed. Subsequently, a prospective clinical trial to evaluate the potential to use the 4-marker panel for therapeutic stratification is necessary.

Here, we demonstrated a 4-gene promoter methylation marker panel consisting of GREM1, NEURL, LAD1, and NEFH that predicts disease outcome in 3 independent patient series. Despite the fact that prognostic significance of the 4-marker panel was lost in multivariate analysis in one patient series, the 4-marker panel seems to be a promising and robust prognostic tool and might be used to improve current prognostic models. However, further research is needed to evaluate the potential use of the 4-marker panel in clinical practice.

W. Van Criekinge is a consultant/advisory board member for MDxHealth. No potential conflicts of interest were disclosed by the other authors.

Conception and design: I.J.H. van Vlodrop, L. Van Neste, M.M.L.L. Baldewijns, L.J. Schouten, W. Van Criekinge, M. van Engeland

Development of methodology: I.J.H. van Vlodrop, T. De Meyer, K.M. Smits, L. Van Neste, K.E. Schuebel, N. Ahuja, M. van Engeland

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): I.J.H. van Vlodrop, S.C. Joosten, L.J. Schouten, P.A. van den Brandt, J. Jeschke, W. Van Criekinge

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.C. Joosten, T. De Meyer, K.M. Smits, L. Van Neste, L.J. Schouten, J.M. Yi, N. Ahuja, F.T. Bosman, W. Van Criekinge, M. van Engeland

Writing, review, and/or revision of the manuscript: I.J.H. van Vlodrop, S.C. Joosten, T. De Meyer, K.M. Smits, L. Van Neste, V. Melotte, M.M.L.L. Baldewijns, L.J. Schouten, P.A. van den Brandt, J. Jeschke, N. Ahuja, J.G. Herman, M.J. Aarts, F.T. Bosman, M. van Engeland

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): I.J.H. van Vlodrop, L.J. Schouten

Study supervision: L. Van Neste, V. Melotte, M. van Engeland

We thank Dr. Egbert Oosterwijk from Nijmegen Center for Molecular Life Sciences (NCMLS), Radboud University Nijmegen Medical Center, Nijmegen, The Netherlands for providing the renal cell carcinoma cell lines. We acknowledge Prof. Dr. Evelyne Lerut and Prof. Dr. Hendrik van Poppel for providing RCC patient material from the University Hospital of Leuven and Geert Trooskens from Ghent University for bioinformatics assistance. We thank Sean Smith, Johanna Nedele, Musinu Zakari, Sarah De Keulenaer, and Jean-Pierre Renard for their technical support.

This work was supported by Ghent University (Multidisciplinary Research Partnership) “Bioinformatics: from nucleotides to networks” and the Susan G. Komen and Mary Kay Ash Foundation.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Ferlay
J
,
Soerjomataram
I
,
Ervik
M
,
Dikshit
R
,
Eser
S
,
Mathers
C
, et al
GLOBOCAN 2012 v1.0, Cancer incidence and mortality worldwide: IARC CancerBase No. 11 [Internet]
.
Lyon, France
:
International Agency for Research on Cancer
; 
2013
[cited 2014 Nov 18]. Available from
: http://globocan.iarc.fr.
2.
Rydzanicz
M
,
Wrzesinski
T
,
Bluyssen
HA
,
Wesoly
J
. 
Genomics and epigenomics of clear cell renal cell carcinoma: recent developments and potential applications
.
Cancer Lett
2013
;
341
:
111
26
.
3.
American Cancer Society
. 
Cancer facts & figures 2015
.
Atlanta, GA
:
American Cancer Society
; 
2015
.
4.
Kattan
MW
,
Reuter
V
,
Motzer
RJ
,
Katz
J
,
Russo
P
. 
A postoperative prognostic nomogram for renal cell carcinoma
.
J Urol
2001
;
166
:
63
7
.
5.
Frank
I
,
Blute
ML
,
Cheville
JC
,
Lohse
CM
,
Weaver
AL
,
Zincke
H
. 
An outcome prediction model for patients with clear cell renal cell carcinoma treated with radical nephrectomy based on tumor stage, size, grade and necrosis: the SSIGN score
.
J Urol
2002
;
168
:
2395
400
.
6.
Zisman
A
,
Pantuck
AJ
,
Dorey
F
,
Said
JW
,
Shvarts
O
,
Quintana
D
, et al
Improved prognostication of renal cell carcinoma using an integrated staging system
.
J Clin Oncol
2001
;
19
:
1649
57
.
7.
Zisman
A
,
Pantuck
AJ
,
Wieder
J
,
Chao
DH
,
Dorey
F
,
Said
JW
, et al
Risk group assessment and clinical outcome algorithm to predict the natural history of patients with surgically resected renal cell carcinoma
.
J Clin Oncol
2002
;
20
:
4559
66
.
8.
Kaelin
WG
 Jr
. 
The von Hippel-Lindau tumor suppressor gene and kidney cancer
.
Clin Cancer Res
2004
;
10
:
6290S
5S
.
9.
Schraml
P
,
Struckmann
K
,
Hatz
F
,
Sonnet
S
,
Kully
C
,
Gasser
T
, et al
VHL mutations and their correlation with tumour cell proliferation, microvessel density, and patient prognosis in clear cell renal cell carcinoma
.
J Pathol
2002
;
196
:
186
93
.
10.
Smits
KM
,
Schouten
LJ
,
van Dijk
BA
,
Hulsbergen-van de Kaa
CA
,
Wouters
KA
,
Oosterwijk
E
, et al
Genetic and epigenetic alterations in the von hippel-lindau gene: the influence on renal cancer prognosis
.
Clin Cancer Res
2008
;
14
:
782
7
.
11.
Varela
I
,
Tarpey
P
,
Raine
K
,
Huang
D
,
Ong
CK
,
Stephens
P
, et al
Exome sequencing identifies frequent mutation of the SWI/SNF complex gene PBRM1 in renal carcinoma
.
Nature
2011
;
469
:
539
42
.
12.
Guo
G
,
Gui
Y
,
Gao
S
,
Tang
A
,
Hu
X
,
Huang
Y
, et al
Frequent mutations of genes encoding ubiquitin-mediated proteolysis pathway components in clear cell renal cell carcinoma
.
Nat Genet
2012
;
44
:
17
9
.
13.
Pena-Llopis
S
,
Vega-Rubin-de-Celis
S
,
Liao
A
,
Leng
N
,
Pavia-Jimenez
A
,
Wang
S
, et al
BAP1 loss defines a new class of renal cell carcinoma
.
Nat Genet
2012
;
44
:
751
9
.
14.
Kapur
P
,
Pena-Llopis
S
,
Christie
A
,
Zhrebker
L
,
Pavia-Jimenez
A
,
Rathmell
WK
, et al
Effects on survival of BAP1 and PBRM1 mutations in sporadic clear-cell renal-cell carcinoma: a retrospective analysis with independent validation
.
Lancet Oncol
2013
;
14
:
159
67
.
15.
Hakimi
AA
,
Ostrovnaya
I
,
Reva
B
,
Schultz
N
,
Chen
YB
,
Gonen
M
, et al
Adverse outcomes in clear cell renal cell carcinoma with mutations of 3p21 epigenetic regulators BAP1 and SETD2: a report by MSKCC and the KIRC TCGA research network
.
Clin Cancer Res
2013
;
19
:
3259
67
.
16.
Dalgliesh
GL
,
Furge
K
,
Greenman
C
,
Chen
L
,
Bignell
G
,
Butler
A
, et al
Systematic sequencing of renal carcinoma reveals inactivation of histone modifying genes
.
Nature
2010
;
463
:
360
3
.
17.
Morris
MR
,
Ricketts
CJ
,
Gentle
D
,
McRonald
F
,
Carli
N
,
Khalili
H
, et al
Genome-wide methylation analysis identifies epigenetically inactivated candidate tumour suppressor genes in renal cell carcinoma
.
Oncogene
2011
;
30
:
1390
401
.
18.
Serre
D
,
Lee
BH
,
Ting
AH
. 
MBD-isolated Genome Sequencing provides a high-throughput and comprehensive survey of DNA methylation in the human genome
.
Nucleic Acids Res
2010
;
38
:
391
9
.
19.
van den Brandt
PA
,
Goldbohm
RA
,
van't Veer
P
,
Volovics
A
,
Hermus
RJ
,
Sturmans
F
. 
A large-scale prospective cohort study on diet and cancer in The Netherlands
.
J Clin Epidemiol
1990
;
43
:
285
95
.
20.
van Houwelingen
KP
,
van Dijk
BA
,
Hulsbergen-van de Kaa
CA
,
Schouten
LJ
,
Gorissen
HJ
,
Schalken
JA
, et al
Prevalence of von Hippel-Lindau gene mutations in sporadic renal cell carcinoma: results from The Netherlands cohort study
.
BMC Cancer
2005
;
5
:
57
.
21.
van Vlodrop
IJ
,
Baldewijns
MM
,
Smits
KM
,
Schouten
LJ
,
van Neste
L
,
van Criekinge
W
, et al
Prognostic significance of Gremlin1 (GREM1) promoter CpG island hypermethylation in clear cell renal cell carcinoma
.
Am J Pathol
2010
;
176
:
575
84
.
22.
Schuebel
KE
,
Chen
W
,
Cope
L
,
Glockner
SC
,
Suzuki
H
,
Yi
JM
, et al
Comparing the DNA hypermethylome with gene mutations in human colorectal cancer
.
PLoS Genet
2007
;
3
:
1709
23
.
23.
McGarvey
KM
,
Van Neste
L
,
Cope
L
,
Ohm
JE
,
Herman
JG
,
Van Criekinge
W
, et al
Defining a chromatin pattern that characterizes DNA-hypermethylated genes in colon cancer cells
.
Cancer Res
2008
;
68
:
5753
9
.
24.
Chan
TA
,
Glockner
S
,
Yi
JM
,
Chen
W
,
Van Neste
L
,
Cope
L
, et al
Convergence of mutation and epigenetic alterations identifies common genes in cancer that predict for poor prognosis
.
PLoS Med
2008
;
5
:
e114
.
25.
Smyth
GK
. 
limma: linear models for microarray data
.
In:
Gentleman
R
,
Carey
VJ
,
Huber
W
,
Irizarry
RA
,
Dudoit
S
,
editors
.
Bioinformatics and computational biology solutions using R and Bioconductor
.
New York, NY
:
Springer New York
; 
2005
.
p.
397
420
.
26.
Langmead
B
,
Trapnell
C
,
Pop
M
,
Salzberg
S
. 
Ultrafast and memory-efficient alignment of short DNA sequences to the human genome
.
Genome Biol
2009
;
10
:
R25
.
27.
Derks
S
,
Lentjes
MH
,
Hellebrekers
DM
,
de Bruine
AP
,
Herman
JG
,
van Engeland
M
. 
Methylation-specific PCR unraveled
.
Cell Oncol
2004
;
26
:
291
9
.
28.
Thijssen
VL
,
Brandwijk
RJ
,
Dings
RP
,
Griffioen
AW
. 
Angiogenesis gene expression profiling in xenograft models to study cellular interactions
.
Exp Cell Res
2004
;
299
:
286
93
.
29.
Brinkman
AB
,
Simmer
F
,
Ma
K
,
Kaan
A
,
Zhu
J
,
Stunnenberg
HG
. 
Whole-genome DNA methylation profiling using MethylCap-seq
.
Methods
2010
;
52
:
232
6
.
30.
Lin
J
,
Gan
CM
,
Zhang
X
,
Jones
S
,
Sjoblom
T
,
Wood
LD
, et al
A multidimensional analysis of genes mutated in breast and colorectal cancers
.
Genome Res
2007
;
17
:
1304
18
.
31.
Morris
MR
,
Ricketts
C
,
Gentle
D
,
Abdulrahman
M
,
Clarke
N
,
Brown
M
, et al
Identification of candidate tumour suppressor genes frequently methylated in renal cell carcinoma
.
Oncogene
2010
;
29
:
2104
17
.
32.
Peters
I
,
Dubrowinskaja
N
,
Abbas
M
,
Seidel
C
,
Kogosov
M
,
Scherer
R
, et al
DNA methylation biomarkers predict progression-free and overall survival of metastatic renal cell cancer (mRCC) treated with antiangiogenic therapies
.
PLoS One
2014
;
9
:
e91440
.
33.
Dubrowinskaja
N
,
Gebauer
K
,
Peters
I
,
Hennenlotter
J
,
Abbas
M
,
Scherer
R
, et al
Neurofilament Heavy polypeptide CpG island methylation associates with prognosis of renal cell carcinoma and prediction of antivascular endothelial growth factor therapy response
.
Cancer Med
2014
;
3
:
300
9
.
34.
Mitsui
Y
,
Hirata
H
,
Arichi
N
,
Hiraki
M
,
Yasumoto
H
,
Chang
I
, et al
Inactivation of bone morphogenetic protein 2 may predict clinical outcome and poor overall survival for renal cell carcinoma through epigenetic pathways
.
Oncotarget
2015
;
6
:
9577
91
.
35.
van Vlodrop
IJ
,
Niessen
HE
,
Derks
S
,
Baldewijns
MM
,
van Criekinge
W
,
Herman
JG
, et al
Analysis of promoter CpG island hypermethylation in cancer: location, location, location!
Clin Cancer Res
2011
;
17
:
4225
31
.
36.
Chen
L
,
Al-Awqati
Q
. 
Segmental expression of Notch and Hairy genes in nephrogenesis
.
Am J Physiol Renal Physiol
2005
;
288
:
F939
F52
.
37.
Michos
O
,
Goncalves
A
,
Lopez-Rios
J
,
Tiecke
E
,
Naillat
F
,
Beier
K
, et al
Reduction of BMP4 activity by gremlin 1 enables ureteric bud outgrowth and GDNF/WNT11 feedback signalling during kidney branching morphogenesis
.
Development
2007
;
134
:
2397
405
.
38.
Michos
O
,
Panman
L
,
Vintersten
K
,
Beier
K
,
Zeller
R
,
Zuniga
A
. 
Gremlin-mediated BMP antagonism induces the epithelial-mesenchymal feedback signaling controlling metanephric kidney and limb organogenesis
.
Development
2004
;
131
:
3401
10
.
39.
Mitola
S
,
Ravelli
C
,
Moroni
E
,
Salvi
V
,
Leali
D
,
Ballmer-Hofer
K
, et al
Gremlin is a novel agonist of the major proangiogenic receptor VEGFR2
.
Blood
2010
;
116
:
3677
80
.
40.
Baldewijns
MM
,
Thijssen
VL
,
Van den Eynden
GG
,
Van Laere
SJ
,
Bluekens
AM
,
Roskams
T
, et al
High-grade clear cell renal cell carcinoma has a higher angiogenic activity than low-grade renal cell carcinoma based on histomorphological quantification and qRT-PCR mRNA expression profile
.
Br J Cancer
2007
;
96
:
1888
95
.
41.
Chen
MH
,
Yeh
YC
,
Shyr
YM
,
Jan
YH
,
Chao
Y
,
Li
CP
, et al
Expression of gremlin 1 correlates with increased angiogenesis and progression-free survival in patients with pancreatic neuroendocrine tumors
.
J Gastroenterol
2013
;
48
:
101
8
.
42.
Le Borgne
R
,
Bardin
A
,
Schweisguth
F
. 
The roles of receptor and ligand endocytosis in regulating Notch signaling
.
Development
2005
;
132
:
1751
62
.
43.
Moll
R
,
Moll
I
. 
Epidermal adhesion molecules and basement membrane components as target structures of autoimmunity
.
Virchows Archiv
1998
;
432
:
487
504
.
44.
Lee
MK
,
Cleveland
DW
. 
Neuronal intermediate filaments
.
Annu Rev Neurosci
1996
;
19
:
187
217
.
45.
Kim
MS
,
Chang
X
,
LeBron
C
,
Nagpal
JK
,
Lee
J
,
Huang
Y
, et al
Neurofilament heavy polypeptide regulates the Akt-beta-catenin pathway in human esophageal squamous cell carcinoma
.
PLoS One
2010
;
5
:
e9003
.
46.
Calmon
MF
,
Jeschke
J
,
Zhang
W
,
Dhir
M
,
Siebenkas
C
,
Herrera
A
, et al
Epigenetic silencing of neurofilament genes promotes an aggressive phenotype in breast cancer
.
Epigenetics
2015
;
10
:
622
32
.
47.
Gerlinger
M
,
Rowan
AJ
,
Horswell
S
,
Larkin
J
,
Endesfelder
D
,
Gronroos
E
, et al
Intratumor heterogeneity and branched evolution revealed by multiregion sequencing
.
N Engl J Med
2012
;
366
:
883
92
.
48.
Gerlinger
M
,
Horswell
S
,
Larkin
J
,
Rowan
AJ
,
Salm
MP
,
Varela
I
, et al
Genomic architecture and evolution of clear cell renal cell carcinomas defined by multiregion sequencing
.
Nat Genet
2014
;
46
:
225
33
.
49.
Stewart
GD
,
Powles
T
,
Van Neste
C
,
Meynert
A
,
O'Mahony
F
,
Laird
A
, et al
Dynamic epigenetic changes to VHL occur with sunitinib in metastatic clear cell renal cancer
.
Oncotarget
2016
;
18
:
25241
50
.

Supplementary data