Prognostic CpG methylation biomarkers identified by methylation array in esophageal squamous cell carcinoma patients.

BACKGROUND: Esophageal squamous cell carcinoma (ESCC) is an aggressive cancer with poor prognosis. We aimed to identify a panel of CpG methylation biomarkers for prognosis prediction of ESCC patients.
METHODS: Illumina's GoldenGate methylation array, supervised principal components, Kaplan-Meier survival analyses and Cox regression model were conducted on dissected tumor tissues from a training cohort of 40 ESCC patients to identify potential CpG methylation biomarkers. Pyrosequencing quantitative methylation assay were performed to validate prognostic CpG methylation biomarkers in 61 ESCC patients. The correlation between DNA methylation and RNA expression of a validated marker, SOX17, was examined in a validation cohort of 61 ESCC patients.
RESULTS: We identified a panel of nine CpG methylation probes located at promoter or exon1 region of eight genes including DDIT3, FES, FLT3, NTRK3, SEPT5, SEPT9, SOX1, and SOX17, for prognosis prediction in ESCC patients. Risk score calculated using the eight-gene panel statistically predicted poor outcome for patients with high risk score. These eight-gene also showed a significantly higher methylation level in tumor tissues than their corresponding normal samples in all patients analyzed. In addition, we also detected an inverse correlation between CpG hypermethylation and the mRNA expression level of SOX17 gene in ESCC patients, indicating that DNA hypermethylation was responsible for decreased expression of SOX17.
CONCLUSIONS: This study established a proof-of-concept CpG methylation biomarker panel for ESCC prognosis that can be further validated by multiple cohort studies. Functional characterization of the eight prognostic methylation genes in our biomarker panel could help to dissect the mechanism of ESCC tumorigenesis.

Introduction

Esophageal cancer is one of the most aggressive cancers with an increasing incidence worldwide. Clinical outcome depends on the disease stage 1. Most patients with esophageal cancer are diagnosed at advanced stage frequently 2. There are two predominant subtypes of esophageal cancer, squamous cell carcinoma and adenocarcinoma, which differ both pathologically and epidemiologically. Esophageal squamous cell carcinoma (ESCC) preferentially affects males and exhibits a marked geographical variation in incidence. For example, the regions of highest occurrence are localized in the so-called Asian Esophageal Cancer Belt 3,4, which extends from northern China to northwestern Iran, including Japan and India. In spite of the advances in diagnostic and therapeutic approaches, the prognosis of ESCC remains poor. Therefore, development of sensitive and specific molecular biomarkers for predicting survival outcome is immensely important.

Previous studies have reported that DNA methylation, which usually occurs in CpG dinucleotides, is a major epigenetic modification in mammalian genome. Hypermethylation of CpG islands is an important mechanism to inactivate tumor suppressor genes 5,6. In addition, DNA methylation is frequently not restricted to a single CpG island but affects multiple independent loci, which is reflective of a widespread deregulation of DNA methylation pattern in different types of tumors 7,8. However, only a few hypermethylated genes have been identified in ESCC and characterized for their role in ESCC tumorigenesis 9-11 and there is a limited number of studies focusing on prognostic methylation biomarker by genome-wide methylation assay in ESCC 12,13.

In this study, we identified nine CpG methylation biomarkers located at eight different genes for prognosis prediction of ESCC by genome-wide methylation analysis of a screen cohort of ESCC patients and validated the methylation array results by pyrosequencing quantitative methylation assay in a new cohort of ESCC patients. In addition, we confirmed that this panel of eight-gene methylation biomarkers showed significantly increased methylation levels in tumor tissues compared to their corresponding normal samples. A risk score defined using the eight-gene panel could predict poor prognosis of a subset of ESCC patients.

Materials and Methods

Clinical samples of ESCC patients

A total of 61 surgically resected ESCC patients were recruited from the Cancer center, National Cheng Kung University Hospital, Tainan, Taiwan after obtaining appropriate institutional review board permission and informed consent from the patients. Among them, 40 patients were analyzed for genome-wide methylation and candidate genes selection as a training cohort. A validation cohort of 61 patients, including the original 40 patients and the additional 21 patients, was analyzed using pyrosequencing quantitative methylation assay for further prognostic analyses. These 61 patients were also included in the analysis of SOX17 methylation and mRNA expression. Overall survival was calculated from the day of surgery to the date of death or the last follow-up. The end of the follow-up period was defined as October 2012. The mean follow-up period for all patients was 19.6 months (range 1.5-68.0 months). Tumor typing and disease staging were performed according to the World Health Organization classification and the TNM classification system, respectively. Information on the sex, age, and smoking history of the patients were obtained from hospital records. Surgically resected tumor tissue and corresponding normal tissue samples were collected for DNA methylation and mRNA expression analyses.

Genomic DNA extraction and sodium bisulfite conversion

Genomic DNA from primary esophageal tissue samples were extracted using proteinase K digestion and phenol-chloroform extraction. For methylation analyses, 1 µg genomic DNA was used for bisulfite conversion using the EpiTect Bisulfite kit (Qiagen, Duesseldorf, Germany) according to the protocols provided by the manufacturer. All unmethylated cytosines were converted to uracil, then to thymidine during the subsequent PCR step.

The genome-wide methylation analysis platform

Illumina's GoldenGate Genotyping Assay (1505 CpG dinucleotides for 807 genes; Illumina, San Diego, CA, USA) was adapted for DNA methylation detection. Using a four-probe design, the assay is able to differentiate between methylated and unmethylated sequences based on Bead Array technology. DNA targets were generated using a pair of allele-specific and locus-specific oligonucleotides followed by allele-specific extension and ligation. Universal primers were then used for amplification and the resulting products were hybridized to a bead array at sites bearing complementary sequences. The hybridized targets contained a fluorescent label denoting methylated or unmethylated state for a given locus. This value, also known as the β-value, ranges from 0 (unmethylated) to 1 (fully methylated) 14.

Pyrosequencing assay

To quantify cytosine methylation in individual CpG sites of candidate methylation probes identified by methylation array, bisulfite-converted DNA was analyzed using a pyrosequencing system (PyroMark Q24, Qiagen). Pyrosequencing treats each individual CpG site as a C/T polymorphism and generates quantitative data for relative proportion of methylated versus unmethylated allele. Specific pyrosequencing primer and PCR primer were designed for ''target'' CpG sites in the probes to be analyzed. Bisulfite-modified DNA was dissolved in 20 μl H2O, and 1 µl of DNA template was used for PCR amplification. Hot-start PCR was performed using the PyroMark PCR Kit (Qiagen), and pyrosequencing was carried out in accordance with the manufacturer's protocol (Qiagen). The target CpG sites were evaluated by converting the resulting pyrograms to numerical values for peak heights. Primer sequences are listed in Table . The percentage of methylation was calculated as the mean of all CpG analyzed 15. Methylation levels below and above the mean level of all tumor tissues analyzed were defined as low methylation and high methylation, respectively.

Calculation of risk score by multivariate Cox regression

A stepwise multivariate Cox regression analysis, in which the regression coefficients were obtained for eight genes showing significant correlation with patient survival in the training cohort. The risk score of each patient was derived by a summation of methylation level of all eight genes each multiplied by the corresponding coefficient. The equations used are as follows: Risk score = (DDIT3 methylation value × coefficients) + (FES methylation value × coefficients) + (FLT3 methylation value × coefficients) + (NTRK3 methylation value × coefficients) + (SEPT5 methylation value × coefficients) + (SEPT9 methylation value × coefficients) + (SOX1 methylation value × coefficients) + (SOX17 methylation value × coefficients). The SOX17 methylation value was the average methylation level of two CpG probes in the gene. The risk score was then used to classify patients into high risk group or low risk group using the mean of the overall risk scores as the cutoff value.

RNA extraction and quantitative reverse transcription-PCR (qRT-PCR)

Total RNA was extracted from 61 tissue samples and all cell lines using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). A total of 4 µg RNA was reverse-transcribed into cDNA using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Carlsbad, CA, USA) according to the manufacturer's protocols. qRT-PCR was performed to detect SOX17 gene mRNA expression level using the Corbett Rotor-Gene 3000 (Qiagen) with β-actin as an internal control. The primers used for SOX17 are sense 5' ACGCTTTCATGGTGTGGGCTAAG 3' and antisense 5' GTCAGCGCCTTCCACGACTTG 3'; for β-actin are 5' GGCGGCACCACCATGTACCCT 3' and antisense 5' AGGGGCCGGACTCGTCATACT 3'. cDNA was amplified using the QuantiFast SYBR Green PCR Kit (Qiagen) and thermal cycling conditions comprised of 95°C for 5 minutes followed by 45 cycles at 95°C for 10 seconds and 60°C for 30 seconds. Cycle threshold (Ct), the fractional cycle number at which the amount of amplified target reached a fixed threshold, was determined. The concentration ratios of SOX17 were calculated using 2-⊿Ct (⊿Ct = Ct-SOX17 - Ct-β-actin). Data were presented as fold differences relative to SOX17 expression in tumor versus corresponding normal tissues based on calculations of 2-⊿⊿Ct (⊿⊿Ct = ⊿Ct-tumor - ⊿Ct-normal) 16. The low expression of SOX17 mRNA was defined as 2-⊿⊿Ct value lower than mean of all patients.

Data processing and Statistical analysis

Type III censoring was performed on subjects who were still alive at the end of the study. Overall survival curves were calculated according to the Kaplan-Meier method after the log-rank test. Hazard ratio for survival risk was determined using the Cox regression analysis. The statistical analyses of promoter methylation and RNA expression were performed using Statistical Package for the Social Sciences version 17.0 (SPSS Inc., Headquarters Chicago, IL, USA). Paired t-test was used to compare the relationship between DNA methylation and mRNA level in tumor and corresponding normal tissues from ESCC patients. P<0.05 was considered statistically significant.

Results

Selection of Prognostic CpG Methylation Biomarkers from Methylation Microarray

We collected surgically dissected tumor and corresponding normal tissues from 40 ESCC patients for genome-wide methylation array analysis. Illumina's GoldenGate Genotyping Assay was adapted for methylation detection in DNA isolated from these tissue specimens as a screen cohort. To obtain DNA methylation markers that can be used as predictor of ESCC patients' survival, several procedures were performed as described below. First, we followed the protocol of Byun and associates 17 to screen for qualified probes by removing probes containing single nucleotide polymorphism (SNP), repeat sequence, probes not in CpG island, or probes in X chromosome or non-differential probes having minimum of β>0.9 and maximum of β<0.1 in all samples. The β-value represents the methylation level of each CpG probes on array, ranging from 0 (unmethylated) to 1 (fully methylated). Second, we selected 141 informative probes with large variance of βTumor (top-30% ranked), or probes with larger mean of (Δβ) (Δβ=βTumor-βNormal of the same patient) among the patients. Third, for each probe, a supervised principal components (Superpc) analysis 18 was applied to compare the survival distributions between patients with methylation levels above and below the mean level of all tumor tissues analyzed; 50 probes were chosen based on these tests. Fourth, the hazard ratio was obtained from the univariate Cox regression analyses of 50 probes in 40 patients analyzed in methylation array. Furthermore, among them 28 probes were chosen according to the hazard ratio >1.5 (). Fifth, multivariate Cox regression and Kaplan-Meier survival analyses were performed for the 28 selected probes in 40 patients (Figure . Finally, we obtained nine DNA methylation probes which were significantly associated with poor prognosis of patients (Figure ). These CpG methylation probe biomarkers were located in eight genes including DDIT3, FES, FLT3, NTRK3, SEPT5, SEPT9, SOX1, SOX17(1), and SOX17(2), in which SOX17 gene with two CpG sites showed significance.

Pyrosequencing Validation of Prognostic CpG Methylation Biomarkers in ESCC Patients

Pyrosequencing is a highly reliable, quantitative, and high-throughput method for the analysis of DNA methylation for further clinical application 15. We thus validated the methylation array results by pyrosequencing quantitative methylation assay for the nine selected probes. A high concordance of the methylation level between microarray and pyrosequencing platforms was found in the screen cohort (Figure and ). Therefore, we performed pyrosequencing of the nine selected probes in a total of 61 ESCC patients including those analyzed in methylation array for a more intuitive comparison of the biomarkers. The Kaplan-Meier survival curves in Figure showed the prognostic significance of nine probes of eight genes in the validation cohort of 61 ESCC patients classified by the mean of methylation level of all tumor tissues (P<0.001~P=0.033, log-rank test). Note that the survival curves of SOX17 were shown by the average of methylation level from two probes.

CpG hypermethylation of all eight genes showed an association with poor survival outcome, suggesting that these eight genes may function as tumor suppressor-like genes in ESCC tumorigenesis. Therefore, we performed pyrosequencing of the selected genes in the paired tumor and normal tissues of 61 ESCC patients. A significantly higher methylation level in tumor tissues compared to their corresponding normal samples was found in all genes analyzed (P<0.0001~P=0.026,), suggesting that these prognosis markers involved in ESCC tumorigenesis.

The Risk Score Calculation and Prognostic Prediction Based on the Eight-Gene Panel

Furthermore, we performed the stepwise logistic regression to verify the multi-genes panel with high discrimination accuracy by obtaining the regression coefficients from the Cox model in the training cohort of 40 patients (Figure ). The regression coefficients for each gene were: DDIT3, -0.299; FES, 0.039; FLT3, 2.098; NTRK3, 0.849; SEPT5, 0.552; SEPT9, 2.379; SOX1, 0.557; SOX17, 0.520. The risk score of each patient was derived by a summation of methylation value multiplied by the corresponding coefficient of each gene. The equations used are as follows: Risk score = DDIT3 methylation value × (-0.299) + FES methylation value × 0.039 + FLT3 methylation value × 2.098 + NTRK3 methylation value × 0.849 + SEPT5 methylation value × 0.552 + SEPT9 methylation value × 2.379 + SOX1 methylation value × 0.557 + SOX17 methylation value × 0.520. Moreover, we applied the risk score calculation to prognosis prediction in the validation cohort of 61 patients. Mean of the overall risk scores was chosen as the cutoff value to classify patients into two groups whereby patients below and above the mean of the overall risk scores were defined as better survivors and poorer survivors, respectively. Such a calculation significantly predicted a subset of patient with a high risk score showing poorer survival (P=0.012, Figure ).

In addition, we performed the multivariate Cox regression model to demonstrate that the power of prognostic discrimination with the eight-gene panel remained high in subgroups stratified according to other clinicopathological parameters of the ESCC patients. The Table shows that the methylation level of eight-gene panel was an independent prognosis factor in ESCC patients even after adjusting for other variables including tumor stage, age and cancer recurrence status by the multivariate analysis (hazard ratio, HR, 1.996; 95% confidence interval, CI, 1.153-3.456; P=0.014).

Promoter Methylation is the Predominant Mechanism of SOX17 Gene Low mRNA Expression in ESCC Patients

Note that SOX17 gene was verified by two probes using both microarray and pyrosequencing assays, suggesting the biological importance of DNA methylation of SOX17 gene in ESCC tumorigenesis. Given that aberrant methylation can result in transcriptional silencing of the target gene 19, we further examined the mRNA expression level of SOX17 in 61 ESCC patients. The qRT-PCR results indicated that mRNA level of SOX17 was significantly lower in the high methylation group than the low methylation group of tumor tissues from ESCC patients (P=0.018, Figure ). Importantly, an inverse correlation between SOX17 DNA methylation and mRNA expression was observed (r= -0.418, P<0.01, Figure ), suggesting that DNA hypermethylation of SOX17 resulted in low mRNA expression in ESCC patients. In addition, patient with mRNA expression level was grouped by the mean of the mRNA expression ratio between tumor and corresponding normal tissues of all patients analyzed. Levels below and above the mean were defined as low and normal mRNA expression, respectively. The Kaplan-Meier analysis showed that patients with low mRNA expression group of SOX17 had significantly poorer overall survival compared with normal expression group (P=0.041, Figure ). The similar prognostic results of SOX17 mRNA expression were observed in another cohort of patients (P=0.042, ), confirming that SOX17 DNA methylation and mRNA low expression are both potential prognostic biomarkers for ESCC patients.

Discussion

Despite the recent advancements in diagnostic techniques and multimodality treatments, ESCC still remains one of the most aggressive and lethal diseases. In this study, we collected a cohort of paired tumor and normal tissues of ESCC patients to conduct genome-wide DNA methylation array analyses. The array data were validated for the methylation level of candidate probes by pyrosequencing to obtain more quantitative results in more patients. We identified eight prognostic methylation genes in ESCC. These genes included DDIT3, FES, FLT3, NTRK3, SEPT5, SEPT9, SOX1, SOX17(1) and SOX17(2), in which SOX17 gene with two CpG sites showed significance. The risk score calculated using the eight-gene panel could serve as an independent prognosis biomarker by Cox regression model and the multivariate analysis. Among the eight-gene methylation biomarkers, we evaluated the inverse correlation between DNA methylation and mRNA expression of SOX17 gene. We provided novel prognostic methylation biomarkers in ESCC that can be further validated by multiple cohort studies.

Lima et al. had recently shown the use of Illumina's GoldenGate methylation array to establish methylation profiles of ESCC 20. They identified a total of 37 CpG sites corresponding to 34 genes to be differentially methylated in tumors compared with surrounding tissues in ten ESCC patients. They further compared the methylation patterns of ESCC patients with esophageal mucosa from four healthy individuals and the TFF1 promoter methylation was speculated to represent an early event in the development of ESCC 20. In our study, prognostic methylation biomarkers were derived from Illumina's GoldenGate methylation microarray, stepwise bioinformatics and statistics analyses of 40 ESCC patients. We found 9 CpG loci across eight different genes that could predict ESCC survival risk. The risk score of eight-gene panel was revealed as an independent prognostic marker even after adjusting for other clinicopathological parameters. Indeed, some of our selected prognostic methylation biomarkers had been individually shown to be methylated in other cancers. For example, SEPT9 is a frequently methylated tumor suppressor gene in head and neck squamous cell carcinoma 21,22 In addition, SEPT9 is a methylation biomarker in plasma for colorectal cancer 23. The promoter of DDIT3 gene is frequently hypermethylated in chronic myeloid leukemia patients 24. Note that SOX1 is a DNA methylation marker for cervical cancer and ovarian cancer 25,26. We analyzed the alteration mechanisms of SOX17 at DNA and RNA levels in ESCC patients and found an inverse correlation between SOX17 DNA methylation and mRNA expression in clinical model. These results suggested that DNA hypermethylation is a predominant mechanism for SOX17 gene silencing, which is in agreement with previous reports of esophageal cancer patients determined by methylation-specific PCR method 27,28.

In conclusion, our study provided a proof-of-concept CpG methylation biomarker panel for ESCC prognosis prediction. The panel of eight-gene prognostic methylation biomarkers also showed a significantly higher methylation level in tumor tissues compared to their corresponding normal samples, supporting their role in ESCC tumorigenesis. Further characterization of novel genes validated in our methylation biomarker panel could help to dissect the mechanism of ESCC tumorigenesis and progression. In addition, the high concordance of methylation level between array and pyrosequencing platforms suggested that our selected panel of DNA methylation biomarkers validate for further clinical application by pyrosequencing, which is a highly reliable, quantitative, and high-throughput method for the analysis of DNA methylation at multiple CpG sites. We propose that the prognostic eight-gene panel may allow clinicians to risk‐stratify their ESCC patients for close monitoring of long‐term outcomes. The use of this set of methylation markers combined with endoscopic ultrasound features would identify patients in dire need of intensive surveillance and further adjuvant therapy to seek for possible improvement of their clinical outcome.

Supplementary Table and Supplementary Figures.

Click here for additional data file.

Table 1

The Primers Used for Pyrosequencing Analysisa.

Gene	Primer name	Sequence ( 5'→3' )	PCR size (bp)
DDIT3	Biotin-Forward	GGGACACCGCTGATCGTTTAGGTAAGATTTTTTGGTAAGTGATGG	197
	Reverse	ACCAAACAACAAATCCTATAACTCC
	Sequencing	AACCAAATCTCATCC
FES	Biotin-Forward	GGGACACCGCTGATCGTTTAGTTTAGGGGTTTGGGTTAATTGAAAT	223
	Reverse	ACCCCCACCTAAATAAATCTCTA
	Sequencing	CCCCCACCAACAACTACTCCTAAA
FLT3	Forward	GGAGAGGTTGGGTAGAGT	295
	Biotin-Reverse	GGGACACCGCTGATCGTTTATCAACCCCTTTCTCAAAACCTCAAAAATCC
	Sequencing	GGTTTAGTAGGAGAAGTTGTAT
NTRK3	Biotin-Forward	GGGACACCGCTGATCGTTTAAGGGGGGGAAGGAAATAA	288
	Reverse	CCACCCTCTATTACAACCCTACAAAATCT
	Sequencing	TCAAACCCTAAAAAATTTTACATCT
SEPT5	Biotin-Forward	GGGACACCGCTGATCGTTTAGGTTGGTTAGTGTTTAGTATTGTT	330
	Reverse	TACAAACCCTAACCCTACCTT
	Sequencing	CCTACAACCTACCAAAT
SEPT9	Forward	GGTTGGGGTTAGTTTAGGATAGAG	154
	Biotin-Reverse	GGGACACCGCTGATCGTTTAAAAAACCAACTCACAAAATCAAATT
	Sequencing	GGTTAGTTTAGGATAGAGGA
SOX1	Forward	GGGGTGGAGGTAAAGATAATAAAAAAAGGT	393
	Biotin-Reverse	GGGACACCGCTGATCGTTTATCCCCCCAACACCTAAAACCAACACATA
	Sequencing	GTGTAGTAGTGGTTAAGTTGTA
SOX17(1)	Forward	GTTTGGGAGGGTTGATTGTA	147
	Biotin-Reverse	GGGACACCGCTGATCGTTTAAAATAAACTAAATCACCCACCACTA
	Sequencing	GGTTGATTGTATTTTGGAAT
SOX17(2)	Forward	GGGGTTTTTTTTGTATAGATGTGG	379
	Biotin-Reverse	GGGACACCGCTGATCGTTTAAAACTAAATCACCCACCACTAA
	Sequencing	TTT TGTATAGATGTGGTTAATG

Pyrosequencing was performed for the nine methylation probes identified by methylation array. Pyrosequencing-PCR was amplified by forward and reverse primers with one labeled with biotin to facilitate further purification. Pyrosequencing-sequencing primer was used to generate the quantitative methylation level of fragment generated by Pyrosequencing-PCR.

Table 2

Multivariate Cox regression assay of eight-gene panel and clinicopathological parameters of the validation cohort of 61 ESCC patients.

Characteristics	N	HR (95%CI)a	P value
8-genes panel (classification by mean)
Low risk score (<mean)	36	1.00
High risk score (>mean)	25	1.996 (1.153-3.456)	0.014
Tumor stage
Early (I/II)	15	1.00
Late (III/IV)	46	1.913 (0.968-3.781)	0.062
Age
<65	43	1.00
>65	18	0.747 (0.405-1.378)	0.351
Recurrence
No	25	1.00
Yes	36	1.695 (0.924-3.109)	0.088

a HR: hazard ratio; CI: confidence interval.