Methylome profiling reveals functions and genes which are differentially methylated in serrated compared to conventional colorectal carcinoma

Background Serrated adenocarcinoma (SAC) is a recently recognized colorectal cancer (CRC) subtype accounting for 7.5–8.7 % of CRCs. It has been shown that SAC has a worse prognosis and different histological and molecular features compared to conventional carcinoma (CC) but, to date, there is no study analysing its methylome profile. Results The methylation status of 450,000 CpG sites using the Infinium Human Methylation 450 BeadChip array was investigated in 103 colorectal specimens, including 39 SACs and 34 matched CCs, from Spanish and Finnish patients. Microarray data showed a higher representation of morphogenesis-, neurogenesis-, cytoskeleton- and vesicle transport-related functions and also significant differential methylation of 15 genes, including the iodothyronine deiodinase DIO3 and the forkhead family transcription factor FOXD2 genes which were validated at the CpG, mRNA and protein level using pyrosequencing, methylation-specific PCR, quantitative polymerase chain reaction (qPCR) and immunohistochemistry. A quantification study of the methylation status of CpG sequences in FOXD2 demonstrated a novel region controlling gene expression. Moreover, differences in these markers were also evident when comparing SAC with CRC showing molecular and histological features of high-level microsatellite instability. Conclusions This methylome study demonstrates distinct epigenetic regulation patterns in SAC which are consistent to previous expression profile studies and that DIO3 and FOXD2 might be molecular targets for a specific histology-oriented treatment of CRC. Electronic supplementary material The online version of this article (doi:10.1186/s13148-015-0128-7) contains supplementary material, which is available to authorized users.


Background
Serrated adenocarcinoma (SAC) has been accepted in the latest WHO classification of tumours of the digestive system as a new subtype of colorectal cancer (CRC) [1]. Criteria for its histologic diagnosis have been proposed [2] and recently validated in a series of 81 cases [3]. Its frequency ranges from just 7.5 to 8.7 % of all CRCs [2,3], but according to 2012 Globocan statistics and given the high incidence of CRC from all cancers (9.7 %), SAC would have an incidence similar to that of multiple myeloma (0.8 %) (http://globocan.iarc.fr). It has been shown that SAC has distinct histological and molecular features and a worse prognosis than conventional carcinoma (CC) [3]. Accordingly, our group observed that SAC, compared to CC, displays a higher frequency of adverse histological features at the invasive front [4] that differs in the expression pattern of adhesion molecules [5] and oncogene mutation [6]. Additionally, two previous studies on mRNA profiling have revealed that SAC showed a higher representation of morphogenesis-, hypoxia-, cytoskeletonand vesicle-transport-related functions and also an over-expression of HIF-1α, fascin1 (actin-bundling protein associated with invasion) and the antiapoptotic gene hippocalcin and a downregulation of the morphogenesisrelated proteins EPHB2 and PTCH [7,8]. Moreover, SAC differs from those CRCs displaying histological features of high-level microsatellite instability (hMSI-H) in terms of oncogene mutation prevalence, MSI status and MLH1 expression [6,8].
A two-arm model has been proposed to explain the progression of the serrated pathway, both of which may progress to SAC. Cancers developing from this pathway may show high-or low-level microsatellite instability (MSI-H and MSI-L, respectively) or may be microsatellite stable (MSS) [9]. Sporadic CRC characterized by MSI-H and BRAF mutation displays typical histological features compared to SAC [10][11][12][13] and is considered a different endpoint of the serrated pathway [14].
In the serrated pathway, inhibition of apoptosis and subsequent inactivation of DNA repair genes by promoter methylation appear to play an important role [1,2]. Although aberrant cytosine-phospho-guanine (CpG) methylation has been proposed as the typical leading mechanism for serrated carcinogenesis [14], there is no study analysing methylome profiling in SAC.
In this study, we have analysed the methylome profile of SAC with the aim of: Evaluating the differentially methylated biological functions of SAC and comparing them with previous studies including mRNA profiling [7,8].
Validating the genes differentially methylated in SAC in comparison with CC at the CpG, mRNA and protein level. Studying the presence of the identified biomarkers in hMSI-H and comparing them with those found in SAC.

Results
As a prerequisite for a matched group, CC did not show significant differences in terms of demographic and clinico-pathologic features with SAC in the training and validation sets included in the study (Table 1).

Differentially methylated functions
Noisy methylation measurements and the high number of genes tested in the array, which imposes strong p value corrections, usually hamper the detection of significant methylation differences at the gene level. As this may be the case in our study, we opted for an alternative, yet wellestablished approach to identify regulated functions in high "omics" experiments, namely the Gene Set Enrichment Analysis approach (GSEA) which relies on the ordering of genes according to a molecular phenotype, as differential expression or methylation, rather than on a gene selection based on a pre-defined p value cutoff. We used the GSEA-related method Fatiscan [8] to analyse our data. This approach revealed a considerable number of GO terms differentially methylated in SAC vs. CC: 86 GO biological processes (BP), 31 GO CC and 69 GO molecular functions (MF). Differentially methylated activities were related to ion binding, intracellular transport, actin binding, GTPases and kinase signaling, neural markers, DNA repair and VEGF signaling amongst others. Figure 1 shows the FatiScan annotated function corresponding to the GO biological process and molecular functions, and Additional file 1 represents box plots indicating the number and the percentage of genes belonging to each of these functions. Additional file 2 shows the GO plots and box plots corresponding to the GO cellular component category. In contrast, only 7 GO biological processes, 8 cellular components and 11 molecular functions were differentially methylated when comparing Spanish SAC and Finnish SAC tumour cases (data not shown).

Differentially methylated genes
The analysis of the methylome microarray data identified 15 differentially methylated genes, 14 of which were more methylated in SAC than in CC ( Table 2). No significant methylated genes were observed when comparing normal SAC and CC mucosa or Spanish and Finnish serrated cases.

Validation of methylated sites by MSP and pyrosequencing
Based on the extent of the differential methylation grade, the importance of the biological functions, the design of suitable primers and the availability of antibodies, we decided to validate DIO3 and FOXD2 by methylationspecific PCR (MSP) and pyrosequencing, respectively, and by quantitative polymerase chain reaction (qPCR) and immunohistochemistry (IHC) for both. As regards the first two, the distribution of CpG sites and the flanking sequence of the DIO3 and FOXD2 CpG islands make the validation by MSP more feasible in the case of DIO3 (CpG sites more condensed) and by pyrosequencing in the case of FOXD2 (CpGs more scattered).
Consistent with the microarray results, MSP revealed that DIO3 CpGs were more methylated in SAC than in CC cases (36/59, 61.1 vs. 14/44, 31.8 %; p = 0.003), Additional file 3 showing representative results. When the correlation of DIO3 methylation was evaluated with oncogene mutation and MSI status, no significant associations were found (data not shown) except for a tendency of BRAF mutation with DIO3 methylation (BRAF mutated, 10/13 (79.9 %) vs. BRAF WT, 24/51 (47.1 %); p = 0.054). Concerning FOXD2, CpG sequence analysis revealed three different clusters of CpG sites, two in the 5′ UTR and one in the 3′ UTR. From these three regions, only the 3′ UTR, which comprised nine CpG sites, was found to be substantially methylated (Fig. 2a). The mean and SD of the methylation percentage for each CpG site and study group is shown in Additional file 1. Interestingly, normal tumour-adjacent mucosa showed a lower percentage of methylation (mean = 31.19) than the tumoural group (mean = 52.81), this difference being statistically significant (F (1.106) = 28.636, p < 0.001, partial = 0.213). The non-tumoural cases had a more similar methylation pattern than the tumoural cases (Levene's F > 16, p < 0.001) (Additional file 4A). Figure 3 and Additional file 3 show that SAC displays a higher percentage of methylation compared to that observed for CC in each of the nine CpG sites, although not reaching statistical significance. In order to test which variable was associated with higher FOXD2 methylation, a bivariate analysis was performed. Higher FOXD2 methylation was associated with tumoural specimens, KRAS, BRAF and exon 20 PIK3CA mutation and MSI-H status (Table 3). Therefore, specimen status, KRAS, BRAF and MSI were used as predictors for the multivariate analysis whereas exon 20 PIK3CA mutation was excluded as the sample size was less than 15. Additional file 5 indicates that only tumoural status and MSI remained significant independent predictors. For the normal and MSS cases, the model predicts a 32.8 % occurrence of methylation in R3 with an increase of 17.1 points for tumoural cases and 13.7 for MSI cases. Moreover, external validation with the TCGA database was

Validation by qPCR
With the aim of finding out whether the DIO3 and FOXD2 hypermethylation affected gene activation, an analysis of the mRNA expression by quantitative PCR was performed. As shown in Table 4, all the compared medians were significantly different. The tumoural and nontumoural groups differed in the expression of DIO3 and FOXD2, the non-tumoural group showing a lower median for both variables (0.154 ± 0.18 vs. 0.054 ± 0.134; p = 0.043 and 0.281 ± 0.473 vs. 0.095 ± 0.301; p = 0.013, respectively, Additional file 4B). Amongst tumoural cases, SAC displayed lower expression medians than CC for both genes (DIO3, 0.007 ± 0.095 vs. 0.131 ± 0.155; p = 0.02 and FOXD2, 0.01 ± 0.253 vs. 0.174 ± 0.331; p = 0.041) (Fig. 4). Additionally, a positive relationship between DIO3 and FOXD2 expression was observed. For the whole group (tumoural and non-tumoural, serrated and conventional), Pearson's correlation coefficient was r = 0.792 (p < 0.001, n = 52); this value indicating that the expression of both genes shares almost 63 % of their variance (R 2 = 0.627). If we consider only the tumoural cases, it was r = 0.751 (p < 0.001, n = 43), both genes sharing 56 % of their variance (R 2 = 0.564) (Additional file 6). In fact, CRC data deposited in the TCGA data sets confirmed this direct correlation as significant (p = 0.037) [15] as retrieved via www.cbioportal.org.
In order to test if 3′ UTR FOXD2 CpG sites might be regulating FOXD2 expression, we performed a  Spearman's correlation analysis which showed that all these sites revealed negative r values indicating an inverse relationship between methylation and mRNA expression and statistically significant correlations except for CpG2 (p < 0.01) and CpG9 sites (p = 0.084). This finding was also observed when representing the mean of the methylation of the nine CpG positions with the qPCR results (Fig. 2b, Additional file 7).

Validation by immunohistochemistry
In order to investigate whether differential methylated status of DIO3 and FOXD2 in SAC compared to CC could have an effect on protein expression within the tissue cells, immunohistochemistry for both D3 and FOXD2 was performed. Both markers showed a granular cytoplasmic staining which was more evident at the luminal border (Fig. 5). Whereas no significant differences were observed in the staining distribution of D3 when SAC vs. CC were compared, FOXD2 staining was more diffuse in SAC (p = 0.013) (

Biomarker comparison between SAC and hMSI-H
As shown in Fig. 3, hMSI-H tumour cases were more methylated in FOXD2 than SAC and CC which was not unexpected since the MSI-H feature was associated with, and was an independent factor for, FOXD2 methylation (Table 3, Additional file 5). Consequently, qPCR analysis revealed that FOXD2 expression was less in hMSI-H than that in SAC (p = 0.030) or CC (p < 0.001) and that of DIO3 less in hMSI-H than that in CC (p = 0.002) but without reaching statistical significance in the comparison with SAC (p = 0.243) (Table 4). Interestingly, immunohistochemical analysis revealed that hMSI-H showed less D3 and FOXD2 expression than SAC and CC, both in staining distribution and intensity, especially for FOXD2 (Table 5; Fig. 4)

Discussion
Altered methylation patterns not only serve as important diagnostic and prognostic markers but are also aetiological factors in colonic carcinogenesis [16]. Whereas the pathogenic sequence leading to conventional carcinoma has been extensively studied, the molecular events responsible for the serrated pathway are poorly understood. SAC is considered as one endpoint of the serrated pathway and, since its recognition by the WHO as a different CRC entity, several studies have aimed to characterize it [2][3][4][5][6][7][8]. O' Brien et al. investigated the methylation status of five genes (hMLH1, MGMT, MINT1, MINT2, p16) comprising the so-called CpG-island methylation phenotype (CIMP) in serrated carcinomas and precursor polyps [17] and found that higher CIMP levels (four or more markers positive) were significantly more frequent in advanced Our study suggests that SAC shows differential methylated functions and genes compared to CC and, moreover, the 186 GO differentially methylated functions found are consistent with previously reported differentially expressed functions (GTPases, cytoskeleton, morphogenesis, apoptosis, vesicle transport and neuronal regulation functions) obtained from two microarray studies comparing the mRNA profile of 5 [7] and 11 [8] SACs with that of CCs. In addition, these functions are also supported by different studies characterizing the histological, immunohistochemical and molecular features of SAC [3][4][5][6]12]. Noteworthy, Spanish and Finnish SAC cases only showed 26 differentially methylated GO functions and no differentially methylated gene, thus supporting the criteria for its diagnosis [11] and for a common pathogenic mechanism despite the differences in environmental and genetic background.
As regards the differentially methylated genes, it is important to emphasize that none were observed in the comparison of non-tumoural SAC vs. matched nontumoural CC, therefore supporting the absence of a distinct and relevant baseline methylation alteration. The comparison between SAC and CC showed 15 genes differentially methylated, 14 of which were more methylated in SAC, thus highlighting the involvement of aberrant methylation in SAC compared to CC. In fact, SAC, in comparison to CC, was more frequently CIMP high than CC (11/37, 30 % vs. 3/32, 9 %; p = 0.0035) (unpublished observations). The differentially methylated genes found are consistent with the functions previously assigned to SAC in microarrays studies [7,8] like GTPase activity (OR51G1 and LRRK2), cytoskeleton (RIOK3 and WASF3), vesicle transport (ATP6V1C1), mTOR pathway (PRR5), apoptosis (XKR4), membrane-associated functions (NIPAL4 and DIO3), morphogenesis (DIO3 and FOXD2), immune response (FAM19A5 and WASF3) and neural markers (QPRT).
The only gene more methylated in CC than in SAC was QPRT which codes for an enzyme that was found to detoxify quinolinate, a potent endogenous neuron toxin that is elevated in the brain of patients with neurodegenerative disorders such as Alzheimer's disease [18]. In this respect, previous studies from our group reported the differential expression of neural regulation processes and the upregulation of specific genes such as hippocalcin, a sensor protein that participates in preventing the calcium-driven apoptosis in neurons, of SAC compared to CC [8]. Therefore, this finding suggests that SAC tumoral cells use neuroprotective mechanisms to prevent apoptosis, a common hallmark of this tumour which is responsible for its saw-shaped growth pattern. The validation of DIO3 and FOXD2 at the DNA, mRNA and protein level was carried out based on its level of significance, CpG site distribution and reagent availability. The MSP study confirmed that DIO3 was more methylated in SAC than in CC, and qPCR confirmed that this might result in a decrease in the mRNA level of expression. DIO3 protein, also known as D3, is the type 3 deiodinase which inactivates the tissue effects of thyroid hormones, these playing a pivotal role in the regulation of cell proliferation and morphogenesis of the intestine as observed in animal models [19]. In fact, D3 is elevated in various colon cancer cell lines as well as in human colorectal cancer tissues compared to their normal counterparts and its expression  [20]. Therefore, the less-methylated status and the increased expression of DIO3 in CC cases suggest an involvement of D3 in the downstream transcriptional effects of β-catenin which characterize the CC, but not the SAC carcinogenic process. As opposed to CC, nuclear βcatenin is not expressed in the tumour invasive front [6], nor it also observed in a serrated mice model [21], thus suggesting a lack of implication of this protein in the epithelial-mesenchymal transition in SAC. Other reports have demonstrated that D3, under the control of the Sonic Hedgehog (Shh) pathway, promotes cell proliferation in basal cell carcinomas [22]. Interestingly, whereas the increase in DIO3 methylation in SAC was confirmed with decreasing mRNA expression, the IHC study revealed a higher intensity cytoplasmic staining of D3 in this type of tumour compared to CC. D3 is an integral membrane protein most of which is extracellular thus allowing ready access to circulating thyroid hormones. It has been reported that the way it is inactivated is by internalization from plasma membrane to early endosomes [23]. The granular cytoplasmic staining observed in our immunohistochemical experiments could suggest a negative feedback mechanism of methylation-driven silencing of DIO3 gene when D3 accumulates in the cytoplasm or alternative mRNA transcripts which may not be regulated by this CpG methylation. Possible impairment of D3 trafficking to the membrane or accelerated internalization in SAC could be responsible for this observation. In fact, vesicle transport and cytoskeleton are amongst the differentially expressed functions obtained from the comparison between SAC and CC.
Much less is known about FOXD2, a member of the family of the well-conserved winged helix forkhead transcriptional factors. Previous reports in animal models have demonstrated a role for it in morphogenesis, possibly through Shh signaling [24], and in regulating sensitivity to cAMP in T lymphocytes [25]. Noteworthy, the only two studies of FOXD2 in humans are related to cancer; in one, a locus including the FOXD2 was found to be deleted in meningioma [26] and, in the other, the protein was found to be more highly expressed in prostate cancer and lymph node metastases compared to normal prostate [27]. In our study, we validated the microarray result by qPCR showing that FOXD2 was less expressed in SAC than in CC. By quantifying the relative methylation of FOXD2 CpG island, there was a consistent trend although not reaching statistical significance. Nevertheless, we observed that a region of the 3′ UTR, but not 5′ UTR, appears to be involved in the control of FOXD2 expression and, in fact, we have shown that methylation in this region is strongly correlated with decreased mRNA expression. Moreover, this FOXD2 methylation is related to microsatellite instability, thus the reason why the methylation results did not reach significance as there were MSS CRCs amongst SAC cases. Differences on FOXD2 RNA expression are more profound and consistent than those observed in its methylation, thus suggesting that additional factors to methylation such gene deletion could be influencing FOXD2 expression. Interestingly, FOXD2 maps to chromosome 1p32-34, a locus frequently deleted in the serrated pathway [28]. Whether FOXD2 could be acting as a tumour suppressor gene in MSI-H CRC deserves future studies. FOXD2 showed an increased cytoplasmic expression in SAC compared to CC and, despite being a transcriptional factor, it was not detected in the nucleus of any cases. Hence, as with D3, a negative feedback of cytoplasmic FOXD2 on gene silencing may be implicated. It is important to stress that protein expression assessed by immunohistochemistry does not give information on protein functionality or isoform type. The increased staining of D3 and FOXD2 could be a result of the translation of alternative gene transcripts with a different regulation and stability. In this context, the positive correlation observed between DIO3 and FOXD2 expression and prior evidence of the participation of these two proteins in the Shh pathway might justify a common behaviour in protein expression and functionality. In fact, according to TCGA database, DIO3 expression is associated with only 12 genes, FOXD2 being one of them [15].
Our results also show that SAC and hMSI-H are different CRC entities, the latter showing a higher FOXD2 methylation percentage and lower FOXD2 mRNA expression. The immunohistochemical study showed that hMSI-H, as opposed to SAC, displayed a significant lack of staining for both D3 and FOXD2, thus suggesting a direct effect of methylation upon protein expression and highlighting further molecular differences between these two histological CRC subtypes considered as two endpoints of the serrated pathway.

Conclusions
-SAC and CC are different CRC in terms of methylated functions, thus reinforcing the concept that SAC is characterized by distinct defects in the regulation of morphogenesis, apoptosis, neural markers and Wnt/β-catenin pathway activation. -DIO3 and FOXD2 are two novel genes that are differentially methylated and expressed in CC, SAC and another endpoint of the serrated polyps pathway as it is the CRC showing molecular and histology features of MSI-H, and thus, they may serve as novel interesting diagnostic markers, pathogenic routes tracers or molecular targets. -This study also reports for the first time the association of FOXD2 methylation with CRC development and shows which CpGs seem to be involved in gene expression regulation.

Patients and tumour samples
The clinico-pathological features of the patients have been previously reported [3,4]. Approval for the study was granted by the Santa Lucia University Hospital Ethical Board, and informed consent was obtained from the patients. SACs were diagnosed on the basis of prior established criteria (epithelial serrations, clear or eosinophilic cytoplasm, abundant cytoplasm, vesicular nuclei, absence of, or less than 10 % necrosis of the total surface area, mucin production and cell balls and papillary rods in mucinous areas of a tumour [2]) and so for CCs [13]. and MONO-27) and two pentanucleotide ones (Penta C and Penta D) [6].

DNA extraction
A volume of approximately 10 mm 3 was extracted from each frozen tissue using the disposable sterile biopsy punch. DNA was extracted following the manufacturer's instructions (QIAGEN, Hilden, Germany). Briefly, tissue was disrupted and homogenized in ATL buffer using a Tissueruptor (QIAGEN) incubated with proteinase K, and the homogenate was subjected to automatic DNA extraction using the Qiacube equipment and the QiaAmp DNA Mini Kit (cat n°:51306), both provided by QIAGEN.
Bisulfite treatment and DNA methylation assay

Preprocessing of Methylation data
Processing of raw data was done using R packages. Probes with a low detection p value (p < 0.01) in more than 95 % of the samples and those measuring SNPs or mapping in X or Y chromosomes were removed and normalization followed a three-step procedure. Firstly, a colour bias adjustment was applied using the methylumi R-package [29]. Then, wateRmelon [30] R-package was used to perform between-sample normalization by equalization of type I and type II backgrounds followed by separated quantile normalization of methylated and unmethylated intensities. Finally, A BMIQ [31] intra-sample normalization procedure, included in the wateRmelon R-package, was applied to correct the bias of type II probe values.

Functional profiling
Functional profiling of the differentially methylated genes was performed using the FatiScan method included in the Babelomics [32,33] web suite. The Gene Ontology (GO) [34] database was used in this functional profiling analysis. Results of functional profiling were represented by two graphs for each considered GO category: biological process, molecular function and cellular component. Significant GO terms of each category are clustered applying hclust function included on R [35], using Cohen's kappa value [36] as a measure of the similarity between two terms. This value measures the agreement between two GO terms as regards shared and exclusive items between them.

Differential methylation analysis
The analysis of differentially methylated genes was performed using limma [32] R-package. Data were fitted to a linear model and differential methylated genes were identified by using the empirical Bayes method included in the package. If the comparison was done between paired samples, a moderated paired t test was applied. A FDR-corrected p value of 0.05 was used as the threshold to select differentially methylated genes.

MSP
Converted DNA from the study cases was subjected to MSP CpG islands analysis. Information on primer sequences and annealing temperatures is provided as Additional file 3. After amplification, PCR products were subjected to electrophoresis using the QIAxcel equipment and the QIAxcel DNA High Resolution Cartridge (ref: 929002, QIAGEN, Hilden, Germany).
Pyrosequencing DNA methylation signatures of three different CpG island regions were analysed and quantified using pyrosequencing. Details of primer sequences and PCR are provided in Additional file 3. PCR products were verified using the QIAxcel DNA high-resolution electrophoresis system. Pyrosequencing of methylated sites was performed using the PyroMark Q24 (QIAGEN) according to the manufacturer's protocol. The methylation level was assessed using the PyroMark Q24 2.0.6 Software (QIAGEN) by which the methylation percentage (mC/mC + C) for each CpG was calculated. The results are presented as the percentage (mean ± SD) of the different CpG sites studied for each of the three regions analysed whose sequences and relative positions are shown as Additional file 7.
Quantitative PCR for assessing mRNA expression RNAs from 18 SACs and 25 CCs, including those from the training set, were extracted with the miRNeasy kit (ref:217004, QIAGEN) and used for validation by qPCR. The retrotranscriptase reaction was performed from a total of 1 μg of DNAseI-treated RNA using the Dy-NAmo cDNA synthesis Kit (ref:F470L) provided by Thermo Scientific (Rockford, IL), and information on the qPCR experiment is provided as Additional file 3. The relative quantitation was done by the 2-ΔCt method using β-actin as the housekeeping gene.

Immunohistochemistry
The validation subset consisted of 63 SC and 64 CC cases matched for gender, age and location, and a representative area of each tumour was selected by one of us (JGS) for a tissue microarray (TMA) construction as previously described [5]. Details on the immunohistochemistry procedure are provided as Additional file 3.
TMA sections of 2.5 μm were stained with anti-D3 and anti-FOXD2 antibodies after confirming a homogeneous tumour cell staining in whole tissue sections.
These markers were evaluated by considering a staining intensity (1 = none or weak staining, 2 = moderate, 3 = strong) and a staining area score (A < one third, B = between one and two thirds, C > two thirds) in a given area. For statistical analysis, both intensity and distribution were considered.

Statistical analysis of validation data
For the analysis of quantification of methylated DNA sequences, the data correspond to a split-plot design with one between-subject factor defining six independent groups of cases (SAC, CC, hMSI-H; tumoural and nontumoural) and one within-subject factor (CpG sites) defining nine repeated measures for every case. Accordingly, we performed two ANOVA SPF-p-q. The first compared the means of the tumoural vs. non-tumoural groups in each of the nine different CpG sites, and the second the means of the six different groups in these sites. For checking the relationship between methylation percentage and binary variables, the t test for independent samples and the Mann-Whitney's U test were used. Statistical significance in the immunohistochemistry study was assessed using Pearson χ 2 or Fisher´s exact test when indicated. Descriptive statistics were computed for real-time PCR. External validation of FOXD2 expression and DIO3-FOXD2 correlation was performed using the TCGA database for colon carcinomas (n = 324) [15]. Statistical analysis was performed using the SPSS (Version 22, Chicago, IL) package.

Availability of supporting data
The data set supporting the results of this article are available in the GEO repository, GSE68060 in http:// www.ncbi.nlm.nih.gov/geo/info/linking.html.

Additional files
Additional file 1: Box plots representing the GO biological process (A) and molecular functions (B) differentially methylated between SAC and CC. Each colour indicates functions belonging to the same