- Open Access
Methylome profiling reveals functions and genes which are differentially methylated in serrated compared to conventional colorectal carcinoma
Clinical Epigenetics volume 7, Article number: 101 (2015)
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.
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.
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.
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) . Criteria for its histologic diagnosis have been proposed  and recently validated in a series of 81 cases . 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) . Accordingly, our group observed that SAC, compared to CC, displays a higher frequency of adverse histological features at the invasive front  that differs in the expression pattern of adhesion molecules  and oncogene mutation . Additionally, two previous studies on mRNA profiling have revealed that SAC showed a higher representation of morphogenesis-, hypoxia-, cytoskeleton- and 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 morphogenesis-related 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) . Sporadic CRC characterized by MSI-H and BRAF mutation displays typical histological features compared to SAC [10–13] and is considered a different endpoint of the serrated pathway .
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 , there is no study analysing methylome profiling in SAC.
In this study, we have analysed the methylome profile of SAC with the aim of:
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.
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 well-established 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  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.
The differentially methylated genes that we found encode transcription factors (FOXD2), kinases (RIOK3), G protein-coupled receptor and GTPases (OR4N5, OR51G1, LRRK2) which are involved in hormone regulation (DIO3), cytoskeleton and vesicle transport (RIOK3, WASF3, ATP6V1C1), apoptosis (XKR4), morphogenesis (FOXD2), regulation of nervous cells (FAM19AS, QPRT) and telomere maintenance (POT1) amongst others (Table 2).
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 methylation-specific 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 also successful showing that MSI-H expressed less FOXD2 than MSI-L/MSS colon carcinomas (p = 0.001) (Additional file 6).
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 non-tumoural 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)  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) (Table 5). In addition, SAC showed higher D3 and FOXD2 intensity staining than CC (p = 0.051, p = 0.054), being statistically significant when strong intensity was compared to combined weak and moderate staining (D3, 45.2 vs. 22.4 %; p = 0.018, FOXD2, 54.8 vs. 30.6 %; p = 0.020).
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)
Altered methylation patterns not only serve as important diagnostic and prognostic markers but are also aetiological factors in colonic carcinogenesis . 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–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  and found that higher CIMP levels (four or more markers positive) were significantly more frequent in advanced stages of the serrated pathway. However, no previous studies have analysed the SAC methylome by means of high-throughput techniques.
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  and 11  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–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  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 non-tumoural 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 . 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 . 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 . 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 is a result of the canonical Wnt/β-catenin pathway activation . 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 , nor it also observed in a serrated mice model , 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 . 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 . 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 , and in regulating sensitivity to cAMP in T lymphocytes . 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  and, in the other, the protein was found to be more highly expressed in prostate cancer and lymph node metastases compared to normal prostate . 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 . 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 .
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.
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 ) and so for CCs . Frozen samples of 39 SACs were retrieved from the Santa Lucia University Hospital, Cartagena Spain (n = 21) and the Oulu University, Finland (n = 18) for CpG methylation profiling and, in addition, adjacent normal mucosa was analysed from 16 SACs (12 Spanish and four Finnish cases). A control group of 34 CCs frozen samples matched with SACs for gender, age, location, Dukes’ stage, WHO grade and mucinous pattern was selected from the same tumour banks (25 Spanish and nine Finnish cases) accompanied by mucosal tissue in 14 of these (10 Spanish and 4 Finnish). MSI-H was present in 20 % and 2.9 % of SAC and CC cases, respectively. The mutation rate in KRAS and BRAF oncogenes, which was previously assessed , was 48.5 and 3 % for CC, 33.3 and 33.3 % for SAC and 11.1 and 66.7 % for hMSI-H, respectively. Validation by qPCR was performed on the methylome series and on an additional set of 24 SACs and 12 CCs frozen tumour cases. Paraffin blocks of 42 SAC and 49 CC cases, included in previous studies [6, 8], were used for immunohistochemical validation of the microarray results. Clinico-pathological features of the cases are shown in Table 1. Additionally, a previously described series of 13 CRCs [6, 8] showing MSI-H molecular and histological features (mucinous, signet-ring cell, and medullary carcinoma, tumour infiltrating and peritumoural lymphocytes, “Crohn-like” inflammatory response, poor differentiation, tumour heterogeneity and “pushing” tumour border ) termed hMSI-H were also studied. These hMSI-H were found mostly in females (69.2 %), mean age 70.5 (±10.0), located in proximal colon, showing Dukes’ B (38.5 %) and C (61.5 %) stage, morphological WHO low grade and 15.4 % displaying a mucinous pattern in more than 50 % of the tumour. MSI status was previously determined, according to the manufacturer’s instructions, in all CRC cases using the MSI Analysis System, version 1.2 provided by Promega (Madison, USA) which includes fluorescence-labeled primers for co-amplification of seven markers including five mononucleotide repeat markers (BAT-25, BAT-26, NR-21, NR-24 and MONO-27) and two pentanucleotide ones (Penta C and Penta D) .
A volume of approximately 10 mm3 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
HumanMethylation450K BeadChip (Illumina, Inc., San Diego, CA), using Infinium HD Methylation assay for genome-wide DNA methylation screening, was employed. In brief, genomic DNA (1000 ng) from each sample was bisulfite converted with the EZ DNA Methylation Kit (Zymo Research, Orange, CA) according to the manufacturer’s recommendations. Bisulfite-treated DNA was isothermally amplified at 37 °C (20–24 h), and the DNA product was fragmented by an endpoint enzymatic process, then precipitated, resuspended, applied to an Infinium Human Methylation450K BeadChip (Illumina, San Diego, CA, USA) and hybridized at 48 °C (16–24 h). The fluorescently stained chip was imaged by the Illumina i-SCAN and Illumina’s Genome Studio program (Methylation Module) was used to analyse BeadArray data to assign site-specific DNA methylation β values to each CpG site.
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 . Then, wateRmelon  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  intra-sample normalization procedure, included in the wateRmelon R-package, was applied to correct the bias of type II probe values.
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)  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 , using Cohen’s kappa value  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  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.
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).
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 DyNAmo 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.
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 . 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 non-tumoural) 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) . 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.
colorectal carcinoma showing typical molecular and histological features of MSI-H
high-grade microsatellite instability
quantitative polymerase chain reaction
World Health Organisation
Hamilton SR, Bosman FT, Boffetta P, Ilyas M, Morreau H, Nakamura SI, et al. Carcinoma of the colon and rectum. In: Bosman FT, Carneiro F, Hruban RH, Theise ND, editors. WHO classification of tumors of the digestive system. Lyon: IARC; 2010. p. 134–46.
Mäkinen MJ. Colorectal serrated adenocarcinoma. Histopathology. 2007;50:131–50.
García-Solano J, Pérez-Guillermo M, Conesa-Zamora P, Acosta-Ortega J, Trujillo-Santos J, Cerezuela-Fuentes P, et al. Clinicopathologic study of 85 colorectal serrated adenocarcinomas: further insights into the full recognition of a new subset of colorectal carcinoma. Hum Pathol. 2010;41:1359–68.
García-Solano J, Conesa-Zamora P, Trujillo-Santos J, Mäkinen MJ, Pérez-Guillermo M. Tumor budding and other prognostic pathological features at invasive margins in serrated colorectal adenocarcinoma: a comparative study with conventional carcinoma. Histopathology. 2011;59:1046–56.
García-Solano J, Conesa-Zamora P, Trujillo-Santos J, Torres-Moreno D, Mäkinen MJ, Pérez-Guillermo M. Immunohistochemical expression profile of β-catenin, E-cadherin, P-cadherin, laminin-5γ2 chain, and SMAD4 in colorectal serrated adenocarcinoma. Hum Pathol. 2012;43:1094–102.
García-Solano J, Conesa-Zamora P, Carbonell P, Trujillo-Santos J, Torres-Moreno D, Pagán-Gómez I, et al. Colorectal serrated adenocarcinoma shows a different profile of oncogene mutations, MSI status and DNA repair protein expression compared to conventional and sporadic MSI-H carcinomas. Int J Cancer. 2012;131:1790–9.
Laiho P, Kokko A, Vanharanta S, Salovaara R, Sammalkorpi H, Järvinen H. Serrated carcinomas form a subclass of colorectal cancer with distinct molecular basis. Oncogene. 2007;26:312–20.
Conesa-Zamora P, García-Solano J, García-García F, Turpin Mdel C, Trujillo-Santos J, Torres-Moreno D, et al. Expression profiling shows differential molecular pathways and provides potential new diagnostic biomarkers for colorectal serrated adenocarcinoma. Int J Cancer. 2013;132:297–307.
Huang CS, Farraye FA, Yang S, O'Brien MJ. The clinical significance of serrated polyps. Am J Gastroenterol. 2011;106:229–40.
Bellizzi AM, Frankel WL. Colorectal cancer due to deficiency in DNA mismatch repair function: a review. Adv Anat Pathol. 2009;16:405–17.
Tuppurainen K, Mäkinen JM, Junttila O, Liakka A, Kyllönen AP, Tuominen H, et al. Morphology and microsatellite instability in sporadic serrated and non-serrated colorectal cancer. J Pathol. 2005;207:285–94.
Stefanius K, Ylitalo L, Tuomisto A, Kuivila R, Kantola T, Sirniö P, et al. Frequent mutations of KRAS in addition to BRAF in colorectal serrated adenocarcinoma. Histopathology. 2011;58:679–92.
Redston M. Epithelial neoplasms of the large intestine. In: Odze RD, Goldblum JR, Crawford JM, editors. Surgical Pathology of the GI tract, liver, biliary tract, and pancreas. Philadelphia: Saunders; 2004. p. 441–72.
Bettington M, Walker N, Clouston A, Brown I, Leggett B, Whitehall V. The serrated pathway to colorectal carcinoma: current concepts and challenges. Histopathology. 2013;62:367–86.
Network CGA. Comprehensive molecular characterization of human colon and rectal cancer. Nature. 2012;487:330–7.
Esteller M. Cancer epigenomics: DNA methylomes and histone-modification maps. Nat Rev Genet. 2007;8:286–98.
O'Brien MJ, Yang S, Mack C, Xu H, Huang CS, Mulcahy E, et al. Comparison of microsatellite instability, CpG island methylation phenotype, BRAF and KRAS status in serrated polyps and traditional adenomas indicates separate pathways to distinct colorectal carcinoma end points. Am J Surg Pathol. 2006;30:1491–501.
Thaker AI, Rao MS, Bishnupuri KS, Kerr TA, Foster L, Marinshaw JM, et al. IDO1 metabolites activate β-catenin signaling to promote cancer cell proliferation and colon tumorigenesis in mice. Gastroenterology. 2013;145:416–25.
Ishizuya-Oka A, Shi YB. Molecular mechanisms for thyroid hormone-induced remodeling in the amphibian digestive tract: a model for studying organ regeneration. Dev Growth Differ. 2005;47:601–7.
Dentice M, Luongo C, Ambrosio R, Sibilio A, Casillo A, Iaccarino A, et al. β-Catenin regulates deiodinase levels and thyroid hormone signaling in colon cancer cells. Gastroenterology. 2012;143:1037–47.
Davies EJ, Marsh Durban V, Meniel V, Williams GT, Clarke AR. PTEN loss and KRAS activation leads to the formation of serrated adenomas and metastatic carcinoma in the mouse intestine. J Pathol. 2014;233:27–38.
Dentice M, Luongo C, Huang S, Ambrosio R, Elefante A, Mirebeau-Prunier D, et al. Sonic hedgehog-induced type 3 deiodinase blocks thyroid hormone action enhancing proliferation of normal and malignant keratinocytes. Proc Natl Acad Sci U S A. 2007;104:14466–71.
Baqui M, Botero D, Gereben B, Curcio C, Harney JW, Salvatore D, et al. Human type 3 iodothyronine selenodeiodinase is located in the plasma membrane and undergoes rapid internalization to endosomes. J Biol Chem. 2003;278:1206–11.
Wu SC, Grindley J, Winnier GE, Hargett L, Hogan BL. Mouse mesenchyme forkhead 2 (Mf2): expression, DNA binding and induction by sonic hedgehog during somitogenesis. Mech Dev. 1998;70:3–13.
Johansson CC, Dahle MK, Blomqvist SR, Grønning LM, Aandahl EM, Enerbäck S, et al. A winged helix forkhead (FOXD2) tunes sensitivity to cAMP in T lymphocytes through regulation of cAMP-dependent protein kinase RIalpha. J Biol Chem. 2003;278:17573–9.
Sulman EP, White PS, Brodeur GM. Genomic annotation of the meningioma tumor suppressor locus on chromosome 1p34. Oncogene. 2004;23:1014–20.
van der Heul-Nieuwenhuijsen L, Dits NF, Jenster G. Gene expression of forkhead transcription factors in the normal and diseased human prostate. BJU Int. 2009;103:1574–80.
Jass JR. Serrated route to colorectal cancer: back street or super highway? J Pathol. 2001;193:283–5.
Davis S, Du P, Bilke S, Triche T, Jr., Bootwalla M. methylumi: Handle Illumina methylation data. 2014. R package version 2.10.0.
Schalkwyk LC, Pidsley R, Wong CC, Touleimat wfcbN, Defrance M, Teschendorff A, Maksimovic J. wateRmelon: Illumina 450 methylation array normalization and metrics. 2013. R package version 1.4.0.
Teschendorff AE, Marabita F, Lechner M, Bartlett, Tegner J, Gomez-Cabrero D, et al. A beta-mixture quantile normalisation method for correcting probe design bias in Illumina Infinium 450k DNA methylation data. Bioinformatics. 2013;29:189–96.
Al-Shahrour F, Carbonell J, Minguez P, Goetz S, Conesa A, Tárraga J, et al. Babelomics: advanced functional profiling of transcriptomics, proteomics and genomics experiments. Nucleic Acids Res. 2008;36:W341–6.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.
Al-Shahrour F, Arbiza L, Dopazo H, Huerta-Cepas J, Mínguez P, Montaner D, et al. From genes to functional classes in the study of biological systems. BMC Bioinformatics. 2007;8:114.
R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2014. http://www.R-project.org/. Last access date: November, 2014.
Cohen JA. A coefficient of agreement for nominal scales. Educ Psychol Meas. 1960;20:37–46.
This work was supported by a grant from Instituto de Salud Carlos III, Ministerio de Sanidad, Spain (ref: PI12-1232). We are grateful to the Madrid division of Spanish Genotyping National Center (Cegen) in the National Oncology Research Center (CNIO), Madrid, Spain; to Fundación para la Formación e Investigación Sanitarias from Healthcare Council of Murcia Region, Spain; to Ginés Luengo from Centro Regional de Hemodonación, Murcia, Spain for statistical assistance and to Dr. Mike V. Tobin for reviewing the English version of this manuscript.
The authors declare that they have no competing interests.
JGS, DTM and JW participated in data analysis. MCT and PSL in data analysis and generation of figures. EE in data interpretation and generation of figures, AT and MM in data collection and interpretation and MPG and AC in writing the manuscript and PCZ in the study design, data collection and interpretation and writing of the manuscript. All authors had final approval of the submitted and published versions.
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 cluster, bar length the number of genes mapping this function and the percentage their contribution to the total of genes included in the function.
Significant Gene Ontology cellular components differentially methylated between SAC and CC. In the GOplot (A), each node shows a significant function and its size the grade of significance. Red contours around the nodes indicate that this function is more represented in SAC whereas blue signify CC. Nodes are grouped in clusters showing similar functions. The number for each cluster shows the amount of unique genes for this cluster. The different functions are grouped according to the concordance Kappa value based on the number of shared genes between functions (only lines representing a Kappa > 0.2 are depicted and line thickness indicates higher Kappa). In the box plot (B), each colour indicates the functions belonging to the same cluster. Bar length represents the number of genes mapping this function, and the percentage indicates the contribution of these genes to the total of genes included in the GO function.
Primer sequences used and PCR conditions for the validation by MSP (including a representative result), CpG pyrosequencing (including the result table), quantitative PCR and details on immunohistochemistry and statistical analysis.
Comparison of non-tumoural compared to tumoural cases in terms of methylation percentage of the nine CpGs at the 3′ UTR FOXD2 region (A) and mRNA DIO3 and FOXD2 expression (B). Note the different standard deviation from these two groups in A. Asterisk indicates statistical significance.
Multivariate analysis of the clinicopathological and molecular factors associated with FOXD2 methylation (A); DIO3 and FOXD2 mRNA expression in hMSI-H tumoural and normal specimens (B) and external validation using TCGA database showing that MSI-H expressed less FOXD2 than MSI-L/MSS colon carcinomas (C).
Positive correlation between mRNA expression of DIO3 and FOXD2 in both non-tumoural and tumoural tissues.
Distribution of the CpG sites studied in the FOXD2 gene and correlation table of the methylation percentages at the 3´ UTR CpG sites with FOXD2 mRNA expression.