Chordoma Characterization of Significant Changes of the DNA Methylation Pattern

Chordomas are rare mesenchymal tumors occurring exclusively in the midline from clivus to sacrum. Early tumor detection is extremely important as these tumors are resistant to chemotherapy and irradiation. Despite continuous research efforts surgical excision remains the main treatment option. Because of the often challenging anatomic location early detection is important to enable complete tumor resection and to reduce the high incidence of local recurrences. The aim of this study was to explore whether DNA methylation, a well known epigenetic marker, may play a role in chordoma development and if hypermethylation of specific CpG islands may serve as potential biomarkers correlated with SNP analyses in chordoma. The study was performed on tumor samples from ten chordoma patients. We found significant genomic instability by Affymetrix 6.0. It was interesting to see that all chordomas showed a loss of 3q26.32 (PIK 3CA) and 3q27.3 (BCL6) thus underlining the potential importance of the PI3K pathway in chordoma development. By using the AITCpG360 methylation assay we elucidated 20 genes which were hyper/hypomethylated compared to normal blood. The most promising candidates were nine hyper/hypomethylated genes C3, XIST, TACSTD2, FMR1, HIC1, RARB, DLEC1, KL, and RASSF1. In summary, we have shown that chordomas are characterized by a significant genomic instability and furthermore we demonstrated a characteristic DNA methylation pattern. These findings add new insights into chordoma development, diagnosis and potential new treatment options.


Introduction
Chordomas are malignant tumors with a phenotype that recapitulates the notochord. These tumors arise within the bones of the axial skeleton and show a destructive growth [1,2]. Chordomas are typically largely resistant to conventional chemoand radiotherapy and therefore surgery remains the main treatment option. However, the critical anatomic location and the commonly large tumor size rarely allow a wide curative excision. Therefore recurrent disease is a common event and even metastases have been reported in up to 40% of cases [3]. The molecular and genetic events involved in the development and progression of chordomas are not well understood and biomarkers do not exist. Although chordomas harbor common chromosomal gains and losses [4] they lack balanced or unbalanced chromosomal exchanges. Those lead to the creation of fusion genes and also screening for mutations in brachyury (a nuclear transcription factor highly expressed in chordomas) and other common cancer associated genes like KRAS and BRAF which failed to show a consistent genetic profile. DNA methylation is a tightly regulated process during normal development and it becomes deregulated during neoplastic transformation and disease development [5].
DNA methylation is relatively stable in body fluids like serum or plasma and can therefore be easily detected by sensitive PCRbased assay [6]. Hypomethylation and/or hypermethylation of specific gene loci, including tumor suppressor genes are strongly associated with disease development [7]. DNA methylation of cytosine at CpG islands can function as transcription repressor, which subsequently leads to the silencing of the associated genes.
To the best of our knowledge epigenetic data on chordomas are not available. Therefore we decided to explore if DNA methylation, a well known epigenetic marker, may play a role in chordoma development and if hypermethylation of specific CpG islands may serve as potential biomarkers correlated with single nucleotide polymorphisms (SNP) analyses in chordoma.

Patient samples
The Caucasian study-group included ten chordoma specimens obtained from four male and six female patients. The age of patients at time of diagnosis was between 25 to 75 years (average age 59.7). Tumors were located in the skull, the sacrum/coccyx and the mobile spine. Tumor-volume ranged from 1.5 to 668.2 cm 3 (average 146). All ten chordomas were morphological and histological classified as classic chordomas. The follow-up period ranged from 1 to 113 months (average 41.9). All patients included in the present study were treated by surgery. Seven patients had an intralesional resection, two patients a wide, and one patient a marginal resection. Three out of ten patients received an irradiation-therapy. During the follow-up half of the patients developed a chordoma recurrence. Two patients showed lung metastases. At the end of the follow-up period four patients were DOD (death of disease), one patient suffered a DOC (death of other cause), three patients were AWD (alive with disease), and two patients had NED (no evidence of disease). The research is an original one, presently not under consideration for publication elsewhere, free of conflict of interest and conducted by the highest principles of human subjects. The study protocol and the consent of the informed patients were approved by the ethics committee of the Medical University Graz (vote #18-192ex06/07; valid until 17.04.2013). No research outside Austria was conducted. All patients were informed in detail and have given their written approval.
Affymetrix SNP 6.0 array processing and analysis Genomic DNA was isolated from chordoma tumor tissue and primary peripheral blood cells using the QIAmp DNA Kit (Qiagen, Hilden, Germany). Affymetrix GeneChip Human Mapping SNP 6.0 arrays were performed as described in the Genome-Wide Human SNP Nsp/Sty 6.0 User Guide (Affymetrix Inc., Santa Clara, CA). SNP 6.0 data were imported and normalized using the Genotyping Console 4.0 program default settings. All samples passing QC criteria were subsequently genotyped using the Birdseed (v2) algorithm. We used 60 raw HapMap data generated with the Affymetrix Genome-Wide Human SNP Array 6.0 as reference. Data were obtained from Affymetrix web site and used for normalization. For visualization of Copy Number state and LOH Chromosome Analysis Suite 1.1 software was used.

DNA methylation analyses
The digestion of 600 ng genomic DNA with methylationsensitive restriction enzymes (MSRE) was performed overnight at 37uC by employing a mixture of 6 units of each AciI (New England Biolabs, Frankfurt, Germany), Hin6I (Fermentas, St. Leon-Rot, Germany) and HpaII (Fermentas). Completion of digestion was confirmed by using a control PCR covering known differentially methylated and cancer gene regions (DMRs; H19, IGF2, ABL1, PITX2, XIST and FMR1) as published [8]. Then restriction enzymes were heat inactivated at 65uC for 20 min and digested DNA was amplified in 16 multiplex reactions covering a total of 360 59UTR targets using biotinylated reverse primers. Amplicons of the 16 multiplex PCRs were pooled and upon agarose-gel-control mixed with hybridization buffer and hybridized onto the AIT-CpG360 microarray, presenting triplicate spots of amplicon-specific DNA probes. Upon hybridization and stringency washings, the hybridized amplicons were detected via streptavidin-Cy3 fluorescence. Microarrays were scanned and intensity data extracted from images using Genepix6.0 software  Table 3. doi:10.1371/journal.pone.0056609.g001 (AXON). Then data were subjected to statistical analysis using BRB-AT (see section ''data analysis''). Detailed information on AIT-CpG360 design and analyses is available as supplemental info (Suppl. S1); DNA sequences of primers and probes are published [9].
High throughput quantitative PCR analysis for confirming DNA methylation changes qPCR was performed on MSRE-digested DNA for confirmation of AIT-CpG360 microarray analyses in a nanoliter microfluidics device (running 48 qPCR assays of 48 DNA samples in parallel) using the BioMark system (Fluidigm Corporation, San Francisco, CA). qPCR confirmation was conducted upon preamplification of methylation sensitive restriction enzyme digested DNA using a pool of 48 primer pairs. Pre-amplification products were subjected to single gene-specific qPCRs in a BioMark Instrument using the 48.48 nanoliter qPCR devices (Fluidigm Corporation, CA) as outlined in ''Methods S1''. The qPCR ct values were extracted with Real-Time PCR Analysis Software of the BioMark instrument (Fluidigm Corporation). Transformed ''45-Ct'' values were used for data analyses.

Data analysis
Statistical analysis of microarray and qPCR experiments was performed using the BRB-ArrayTools software 3.8.1 developed by Dr. Richard Simon and the BRB-ArrayTools Development Team (http://linus.nci.nih.gov/brb). Values of AIT-360-CpG-arrays were log 2 -transformed and a global normalization was used to median center the log intensity values within one experiment. To identify genes, differentially methylated between patient-sample classes, a random-variance t-test for paired samples was applied to both data sets [10] [11]. Genes were considered statistically significant, if the parametric p-value was less than 0.01. Significance of differentially methylated genes was ranked using the p-value of the univariate test. In addition the false discovery rate (FDR) was calculated using the method of Benjamini and Hochberg as provided within BRB-ArrayTools software. For defining classifiers with potentially diagnostic value, ''class prediction'' analyses were conducted in BRB and classifiers defined by leaving one out cross validation (see also the BRB website: http://linus.nci.nih.gov/brb/TechReport.htm).

Results
Chromosomal copy number variation (CNV) analysis using Affymetrix 6.0CNV/SNP arrays Ten chordoma samples were tested for copy number (CN) and LOH using Affymetrix 6.0 CNV/SNP Arrays. Most common (.50% of the samples) chromosomal CN gains were observed for 1q21.1-q44, 7q36.3, 14q32.33, and 22q11.22 and losses for 3p26.3-q29, 9, 13q12.11-q22.1, and 22q12-q13.2. The most common chromosome loss involved chromosome 3 where 6 of 10 patients showed a loss of 3p25.2 (RAF1) and all 10 patient showed a loss of 3q26.32 (PIK 3CA) and 3q27.3 (BCL6). Table 1 summarizes the most common CN in ten chordoma patients. The cross-linking of interesting genes is shown in Figure 1 using IPA (Ingenuity Pathway Analysis) software. Furthermore, copy numbers were matched with methylation data and presented in Figure 2 to see whether a chromosome is particularly affected by CN-variation or hyper/hypo methylation pattern.

Identification of DNA methylation changes in chordoma
We analysed 36 DNA samples and 3 negative controls using the AITCpG360 methylation assay. The aim was to identify biomarkers for serum-based patient testing. Therefore we also included healthy blood samples from volunteers in our analyses. For the identification of genes differentially methylated in chordoma versus normal blood we used ''class comparison'' using a cut off value on the single gene level of p,0.01 elucidated 20 genes. Four of them showed p-values below 0.001 (HIC1, CTCFL, ACTB, RASSF1). Based on the geometric mean of the chip intensities from the class of blood samples and chordoma samples the fold change between classes ranged from 0.024-3.82. Values below zero indicate hypermethylation in chordoma versus peripheral blood (inverted values range from 41.66 to 0.026 fold increase in intensities in chordoma ( Table 2). It is of utmost interest for serum-cfDNA methylation based diagnostic testing of clinically suspected patients suffering from chordoma to elucidate a classifier for proper distinction between the methylation pattern of chordoma and blood-DNA to avoid false positives due to the background blood-DNA which is very likely to be the most abundant DNA population present in cell free serum. For identification and building a classifier for ''prediction'' of novel samples we performed ''class prediction''. A feature selection was set to include only genes significantly different between the classes at p,0.01 significance level, and the ''Leave-one-out crossvalidation'' method was used to compute misclassification rate. 94% of samples were correctly classified (sensitivity = 100.0% and specificity = 88.9%; AUC = 0.94) by the gene methylation classifier derived from the ''diagonal linear discriminant analysis'' and also from the ''1-nearest neighbor'' classifier. Other prediction methods (compound covariate predictor, 3-nearest neighbor, nearest centroid, support vector machines and Bayesian compound covariate predictor) made 89% correct classification possible. The classifier genes including summary statistics are listed in Table 3.

qPCR confirmation of DNA methylation changes in chordoma
Analytical qualification of MSRE-coupled qPCR. To reconfirm the microarray-hybridization based analyses we subjected both the undigested and MSRE-digested DNA samples to qPCR analyses using nanoliter scaled microfluidic qPCR arrays in a Fluidigm 48.48 array for quantification of DNA methylation. PCR reactions were redesigned for covering at least 3 MSRE cut sites. On average 6 MSRE sites were present in amplicons and qPCR reactions were qualified according to MIQE guidelines (data not shown). Optimised qPCR conditions enabled parallel analyses of the 20 methylation marker candidate genes. By using the entire capacity of the 48.48 PCR array 28 genes were analysed in addition to the 20 classifier genes.
Of the 48 genes tested, 39 were significant (p,0.05) with an overall mean dCt between digested and undigested sample DNAs of 2.8218.6 (corresponding to 7-416000 fold change) indicating proper digestion for qPCR based elucidation of methylation differences. The amplicons for H19, CDKN2A, IGF2, C3, SRGN, PIWIL4, GBP2, IRF4 showed 0.23-0.36 (in the enlisted order of genes; p = 0.057-0.260) fold differences. DNAJA4 was only minimally changed (0.75 fold), which is in line with the RRBS (reduced representation bisulphite sequencing) and 450 k Infinium Comparison of ''blood and chordoma'' using MSRE coupled qPCR methylation analyses. Group wise comparison of blood (n = 7; four female, three male) and chordoma (n = 10; five male and five female) amplicon Ct values derived from qPCR upon MSRE digestion revealed 10 genes with higher than 2-fold changes corresponding to hypermethylation in tumors (fold change ranged from 322670); in addition to another 10 genes showed more than 2 fold hypermethylation in peripheral blood (a factor of 0.520.13 corresponds to 2-7fold; Table S1).
qPCR-confirmation of the ''classifier derived from chip based screening.''. For confirmation of the 20 classifier genes derived from chip based screening qPCR-Ct values were used for class prediction. Using different classification algorithms, 88-94% of samples were correctly classified; one chordoma and one peripheral blood sample were frequently misclassified by the different prediction tools (Table S2).
For exemplification the performance of the Support Vector Machine Classifier enables correct classification of 94% samples at a sensitivity of 0.889 and a specificity of 1 (one chordoma sample was not correctly classified). The receiver operating characteristics (ROC) derived from the Bayesian Compound Covariate Predictor provides an area under the curve AUC of 0.952. Although the parametric p-values of several single gene qPCR ct values were below p,0.05, the classification success is very impressive.
Generation of a novel classifier from the entire set of 48 qPCR amplicons applying the feature selection criteria ''Genes with univariate misclassification rate below 0.2'' for class prediction elucidates a classifier of 23 genes enabling perfect classification of the entire set of study samples (AUC = 1) by the Compound Covariate Predictor, the 1-Nearest Neighbor and the Bayesian Compound Covariate Predictor. Correct classification of 94% was obtained by using the Diagonal Discriminant, the Nearest Centroid, and the Support Vector Machines analyses. The 3-Nearest Neighbor classification success was 88% (Table S3). For reducing the classifier to a lower number of genes feature selection by ''univariate p-value ,0.05 and 2 fold -change between classes'' was applied and class prediction was performed again on the entire set of all the 48 amplicons used for qPCR. Thereby a classifier for distinction between peripheral blood and chordoma was generated. This classifer consisted of qPCR-ct methylation measures of RASSF1, KL, C3, HIC1, RARB, TACSTD2, XIST, and FMR1 (Table 4). That classifier enabled perfect classification of the set of study samples (AUC = 1) by the 1-Nearest Neighbor method. Correct classification of 94% was obtained by using the Compound Covariate Predictor and the Support Vector Machines. The classification success was 88% achieved by the Diagonal Discriminant Analyses, the Nearest Centroid, and analyses and the 3-Nearest Neighbors classifier. The Bayesian Compound Covariate Predictor allowed also perfect classification. Two samles, however, could not be classified (indicated as ''NA'' in Table S4).

Discussion
The presented study was conducted to investigate SNPs as well as hyper/hypomethylated genes in ten chordoma patients. The first part of the study focused on SNP analyses. We confirmed copy number alterations known in chordomas and describe novel gains and losses. In accordance to Le et al. who reported array comparative genomic hybridization (CGH) in 21 sporadic chordomas and found large copy number losses on chromosomes 1p, 3,4,9,10,13,14, and 18, we were also able to demonstrate Table 3. Composition of the classifier derived from class prediction (Sorted by t -value): HIC1 presented by two different probes on the CpG360 array is present twice in two lines. The letters (A-S) can be found in Figure 2, where SNP data are combined with methylation data. The column ''%CV-support'' in the table indicates the percentage of the cross-validation training sets in which each gene was selected during ''leave-one out cross validation''. 100% means that the gene is so strong that it was selected in all of the cross-validated training sets. ''Geom. mean of intensities'' is derived from chip intensity-values. doi:10.1371/journal.pone.0056609.t003 common losses including chromosome 1,4,9,10,13,14,18,20, and 22 as well as common gains in 7, 12, and 19. Interestingly, in 100% of our cases chromosome 3 (3p26.3-q29) was lost. These findings are also comparable with data previously reported by Hallor et al. and Le et al. [4,12]. Copy number gains/losses present in more than $50% of patients of well-known cancer associated genes for example CDKN2A, CDKN2B, RAF1 and PTEN were found and may play an important role in chordoma development. Overall, chromosome 3 shows the majority of genetic alterations. The common involvement of this chromosome has been already described in early studies by demonstrating a 3p loss [13,14] In addition losses on chromosome 3 relating PIK3CA and BCL6 were obtained. Currently there are several inhibitors of the PI3K (phosphatidylinositol 3-kinase) pathway under investigation in solid tumours [15]. Although cross talk of the PI3K pathway with other pathways in particular the RAS/RAF/MEK pathways have been reported, inhibition of the PI3K pathway could be an attractive therapeutic target and is definitely worth further investigations. BCL6 is a transcriptional repressor binding DNA through zinc fingers and regulates transcription through interacting with other factors like Jun proteins and histone deacetylase family proteins [16,17]. Usually BCL6 is associated with normal and abnormal B-cell development. However, Chamdin et al. showed that BCL6 arrests the differentiation of neural crest cells in neuroblastoma (NB) and may therefore play a similar role in chordoma development [18]. By merging the data, it's apparent that also RB1 (retinoblastoma) signalling plays a central role in chordoma oncogenesis [4,12]. We were able to show that chordomas are characterized by significant genomic instability. Although a common pattern of genetic changes could be demonstrated, a consistent genetic change in all samples was not identified.
The second part of the study provides the first evidences that DNA methylation of tumor suppressor genes exit in chordomas and may serve as a marker for early tumor detection.
Early tumor detection is extremely important for chordoma patients, because these tumors are resistant to chemotherapy and irradiation. Surgical excision remains the main treatment option and based on the challenging anatomic location early detection is important to allow complete resection and to reduce the high incidence of the local recurrence. Therefore, the aim was to identify hypermethylated genes that could serve as biomarkers for early tumor detection to optimize patients' treatment. We used blood from healthy volunteers as comparison, due to the fact that notochord as comparatively tissue was not available.
DNA methylation has already provided useful biomarkers for diagnosing cancer, monitoring treatment and predicting the prognosis. Aberrant DNA hypermethylation of CpG islands in the promoter region of genes is well established as a common mechanism for the silencing of tumor suppressor genes in cancer and serve as an alternative mechanism of functional inactivation.
RASSF1, KL, and HIC1 are known to be tumor suppressor genes. The inactivation of tumor suppressor genes is usually accompanied by a copy of the gene mutations and loss of the corresponding allele [19]. RASSF1 encodes a protein similar to the RAS effector proteins. In normal cells RASSF1 (Ras association domain family1 protein) a tumor suppressor gene is involved in controlling cell cycle and in repairing DNA [20]. RASSF1 has been shown to be transcriptionally silenced by promoter methylation and are frequently methylated in various tumor types. Especially in breast and colorectal cancer [21,22], inactivation of this gene was found to be correlated with CpGisland promoter region hypermethylation.
Another tumor suppressor gene KL (klotho) is a single pass type I transmembrane protein that is localize at the plasma membrane as well as in the cytoplasm. It was initially identified as antisenescence gene [23]. Recently, reduced KL gene expression was shown to contribute to tumorigenesis. KL has been found to function as tumor suppressor in various cancers like breast, pancreas, lung, and cervix [24]. This transmembrane protein can be shed, act as circulating hormone and is a modulator of the IGF-1 (insulin-like growth factor IGF-1) and the FGF (fibroblast growth factor) pathways. Those have recently been demonstrated to be activated in chordomas [25,26]. KL potently inhibits liganddependent activation of the insulin and IGF-1 pathways [27,28] and binds to FGFR (fibroblast growth factor receptors) [28,29].
Another tumor suppressor gene, the HIC1 (Hypermethylated in Cancer 1) gene, is a transcriptional target of p53 and is frequently deleted or hypermethylated in various solid tumors, including colon, lung, breast, brain, and kidney [30]. We have applied this methylation assay to several other cancerous diseases [31][32][33] and could delineate several candidate-biomarker panels for the different settings. Nikolaidis et al. showed the impact of DNA methylation-based assays in the diagnosis of cytologically occult lung neoplasms [34]. HIC1 methylation could be used as a target for pharmacologic DNA-methyltransferase and could therefore suit as a potential new target to treat chordoma patients. In summary, we have shown that chordomas are characterized by significant changes of the DNA methylation pattern. A multigene DNA methylation based classifier suitable to distinguish healthy blood and chordoma DNA presented here will add a new dimension for chordoma diagnosis and treatment. We believe that our findings should be explored to circulating tumor cells or circulating cell free DNA found in peripheral blood, serum or plasma of patients, to improve chordoma diagnoses and disease monitoring. Although validation of results has to be conducted on additional patient sample cohorts and serum cfDNA, we think that the DNA methylation classifiers elucidated here, could be useful novel biomarkers advancing diagnostic workup for patients.