Multilayer OMIC Data in Medullary Thyroid Carcinoma Identifies the STAT3 Pathway as a Potential Therapeutic Target in RET^M918T Tumors

Medullary thyroid carcinoma (MTC) is a malignant tumor of the thyroid gland. It is composed of C cells and accounts for up to 2% of all thyroid cancers (1). Around 25% of MTCs arise due to germline mutations in the "rearranged during transfection" (RET) gene as part of multiple endocrine neoplasia type 2 syndrome, whereas the remaining forms are sporadic. The genomic landscape of recurrent driver alterations found in sporadic disease seems to be mainly restricted to activating mutations affecting the RET (2) or RAS oncogenes (3, 4). There is a distinct genotype–phenotype correlation for each of these alterations (3, 5, 6). The RET^M918T mutation seems to trigger the most aggressive disease, showing poor prognosis with a higher prevalence of distant metastases during follow-up (5, 7). MTC cases with advanced disease are treated with tyrosine kinase inhibitors (TKI), such as vandetanib or cabozantinib, which have shown to have therapeutic benefits. However, many patients suffer secondary toxicities that lead to dose reductions or treatment interruption (8, 9). A comprehensive understanding of the molecular mechanisms driving the different genetic classes of this disease is required to identify new therapeutic opportunities for these patients.

High-throughput technologies have already been used to define genomic signatures linked to particular driver mutations in cell lines and other thyroid cancer subtypes (10–15), uncovering molecular events associated with progression and recurrence. However, data on MTC are scarce due to the disease's low prevalence and, in consequence, the difficulty in collecting an informative sample set. On top of that, C cells represent less than 1% of the thyroid gland, and therefore it is very challenging to portray their normal genomic landscape. To date, only a handful of studies using mRNA expression arrays have been published, reporting that MTC expression profiles are mutation-specific (16–18). Interestingly, the overexpression of genes related to the tumor growth factor-β pathway (17) and tumor invasion and metastases (16, 18) were notable among MTCs caused by the RET^M918T mutation. On the other hand, the miRNA profiling studies published so far have focused on patient outcomes rather than on genetic subtypes of the disease but have nevertheless uncovered some molecular events related to metastasis and worse outcomes (19–21).

Some of the differentially expressed genes identified in the previous profiling studies are undoubtedly regulated by aberrant methylation. However, explorations of this epigenetic mechanism in MTC, the deregulation of which is a well-known hallmark of cancer, have been limited, either focused on specific candidate genes, such as RAS association domain family protein 1 (RASSF1; ref. 22) and Sprouty 1 (SPRY1; ref. 23) or exploring the whole methylome in very few samples (15). Thus, it remains poorly characterized.

Therefore, we studied the role of aberrant DNA methylation in the pathology of genetic subtypes of MTC. We quantitatively profiled the largest cohort of MTCs published to date, composed of 48 samples, on the basis of the DNA methylation levels of more than 27,000 CpGs across the genome. We observed significant differences between the methylation patterns among the samples carrying distinct mutations. Taking advantage of multilayer genomic data consisting of the transcriptome, miRNome, and methylome, we uncovered genes negatively regulated by methylation. Moreover, pathway enrichment analysis underscored JAK/Stat pathway involvement in the RET^M918T-driven tumors. Using a series of 101 formalin-fixed, paraffin-embedded (FFPE) MTC samples, we observed that phosphorylation of the JAK/Stat pathway effector STAT3 (Tyr705) was more evident among RET^M918T-positive tumors. In vitro pharmacologic inhibition of STAT3 phosphorylation in 2 human MTC cell lines (MZ-CRC-1 and TT) revealed a higher dependency of the RET^M918T-mutated MZ-CRC-1 cell line on STAT3 signaling. Furthermore, we observed an apparently synergistic effect between the action of STAT3 inhibitor and Vandetanib, which could be exploited to treat patients with MTC.

Forty-eight fresh-frozen MTC tumor samples were collected at the Spanish National Cancer Research Center (CNIO) in collaboration with the CNIO Tumor Bank. Written informed consent was obtained from all study participants, and the study was approved by the relevant institutional review board (Comité de bioética y bienestar animal of the Instituto de Salud Carlos III). Sections of each sample were evaluated by a pathologist. Only samples in which at least 80% of cells were cancerous were used in this study, and cases with high amyloid content were excluded. A total of 101 paraffin-embedded (FFPE) MTC samples from unrelated individuals were used for validation studies. Twenty-five were used for the validation of DNA functional methylation markers, and the entire collection was used for an immunohistochemical (IHC) study. Hematoxylin and eosin–stained sections of each FFPE tumor sample were examined by a pathologist to select MTC areas representative of each tumor and these were used to construct tissue microarrays (TMA). Three TMAs containing 2 cores from each of 94 tumor samples were constructed. The remaining 7 tumors were evaluated as complete sections.

MZ-CRC-1 and TT cell lines were provided by Dr. James Fagin (Human Oncology & Pathogenesis Program, Memorial Sloan Kettering Cancer Centre, New York, NY) and authenticated by RET sequencing. The MZ-CRC-1 cell line is derived from a metastatic MTC and harbors the RET^M918T mutation, whereas TT cell line harbors a mutation in the 634^th codon of RET. MZ-CRC-1 and TT were cultured in DMEM GlutaMAX (Invitrogen) and F12 Nut Mix medium with GlutaMAX (Invitrogen), respectively, supplemented with 10% (v/v) FBS (PAA Laboratories), 1% (v/v) penicillin/streptomycin, and 0.6% (v/v) Fungizone (Gibco).

Genomic DNA from all samples (frozen and FFPE tumors and cell lines) was extracted using the DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer's protocol. Total RNA was extracted from 5 × 10⁶ of MZ-CRC-1 and TT cells and 1 mL of TRIzol (Life Technologies) using the standard conditions. For immunoblotting, cells lysates were produced in RIPA lysis buffer (Sigma) containing protease inhibitor mixture (Sigma) and Phosphatase Inhibitor Cocktail 2&3 (Sigma).

All samples were screened for hotspot mutations in exons 8, 10, 11, 13, 14, 15, and 16 of RET by Sanger sequencing. If negative, the hotspot codons 12, 13, and 61 of exons 2 and 3 of all RAS genes were screened using the same technique. Primer sequences and examples of the identified mutations can be found in Supplementary Table S1. Tumors were classified as “wild-type” (WT) if no mutation was found.

Genomic DNA was bisulfite-converted using the EZ DNA Methylation Kit (Zymo Research) following the manufacturer's recommended procedures. Genome-wide promoter DNA methylation was determined using the Illumina Infinium HumanMethylation 27K Platform (Illumina) as described previously (24). This assay generates DNA methylation data for 27,578 CpG dinucleotides covering 14,473 genes. For each CpG site, methylation levels were quantified using β-values, which represent the proportion of methylation, calculated as M/(M + U), where M is the methylated probe intensity and U is the unmethylated probe intensity. β-values range from 0 to 1, with 0 being completely unmethylated and 1 being completely methylated. We excluded probes designed for sequences on either the X or the Y chromosome (1,085 or 7 probes, respectively), as well as 46 probes with missing value in at least one sample. The data used in these analyses have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE72729 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72729). The applied analysis workflow is schematically depicted in Supplementary Fig. S1.

Unsupervised hierarchical clustering was carried out using Cluster 3.0 software with “average linkage” (Pearson correlation, uncentered metrics). The clusters were subsequently visualized using Treeview (http://rana.stanford.edu/software). Principal component analysis (PCA) was performed using R CRAN version 2.15.3 (R, 2013). Differences in DNA methylation status among tumor groups (based on the driver mutations present in the samples) were tested using POMELO II, applying linear models (limma; ref. 25). Only groups containing at least 5 tumors were considered. Each genetic class was compared with WT tumors and with each of the other classes. To account for multiple hypothesis testing, P values were adjusted using Benjamini and Hochberg false discovery rate (FDR) correction. We defined a probe to be hypomethylated or hypermethylated when it displayed a mean β-value difference (Δβ-value) <−0.2 or >0.2, respectively, calculated as mean difference among distinct tumor groups, and had an FDR < 0.05. With the extracted lists of differentially methylated genes, we performed Venn diagram analyses using BioVenn (http://www.cmbi.ru.nl/cdd/biovenn/) and Pathway enrichment analysis using the functional annotation clustering tool implemented in DAVID Bioinformatics Resources 6.7 (http://david.abcc.ncifcrf.gov/home.jsp).

Because we had both mRNA (18) and miRNA expression (GEO series accession number GSE72728; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72728] data for 33 and 31 of the tumors used in this study, respectively, we aimed to identify genes whose expression is correlated with the methylation status of their corresponding promoter regions. Expression of all the available protein-coding genes identified as differentially methylated by tumor genotype (FDR < 0.05, |Δβ-value| > 0.2) was examined as previously described (12, 26). The integration of miRNA data required a genome-wide approach, as there was only one microRNA gene promoter showing differential methylation. More specifically, there were 254 probes included on the methylation platform that mapped to promoters of 110 microRNA genes previously annotated by the manufacturer (24). To this number, we added 423 probes belonging to putative promoter regions of 151 miRNAs identified using the PROmiRNA method (27). These additional probes and their corresponding miRNA genes are listed in Supplementary Table S2. Correlation was measured by the Spearman coefficient using R CRAN version 2.15.3 (R, 2013) and genes with a negative correlation coefficient and P < 0.05 were considered significantly correlated.

To confirm the effect of DNA methylation on gene expression, MZ-CRC-1 cells at 60% of confluence were treated with 50 μg/mL 5-aza-2′-deoxycytidine (Sigma). This treatment was repeated every 48 hours over 6 days and then the cells were collected for RNA and DNA extraction.

Bisulfite sequencing of MTC cell lines using the primers listed in Supplementary Table S3 was carried out according to protocols described elsewhere (12) to assess methylation level of the selected promoters (belonging to 3 protein-coding genes and 3 miRNA genes).

For the assessment of expression of the selected protein-coding genes, 1 μg of total RNA was reverse-transcribed using Superscript II (Invitrogen) and an oligo dT₁₄ primer following manufacturer's instructions. The amounts of DKK4, MMP20, and PLCB2 mRNA were quantified by real-time PCR with QuantStudio 6 Flex Real-Time PCR system (Applied Biosystems), using primers designed to be specific for the 3 genes (Supplementary Table S3) and probes from the Universal ProbeLibrary Set, Human (Roche). Normalization was carried out with the internal standard β-actin (ACTB).

For miRNA gene expression quantification, 10 ng of total RNA was used for first-strand cDNA synthesis using the miRCURY LNA Universal RT miR PCR System (Exiqon), and LNA miR-PCR primer/SYBR Green Mix (Exiqon) was used for subsequent quantification of miR-10a, -30a, and -200c according to the manufacturer's recommendations. miR-16 was selected as the reference gene for normalization.

Negative controls were included in all PCR series and assays were carried out in triplicate. The ΔΔC_t method was used for the calculation of mRNA/miRNA content (28).

Three of the most differentially methylated probes, all showing negative correlation with expression of the corresponding protein-coding gene, and with biologic functions potentially relevant to MTC were selected for validation. The technical validation of microarray results in a subset of the original discovery series (comprising 6 RET^M918T, 6 RET^C634, 3 RAS, 2 WT tumors, and 6 tumors bearing other RET mutations) was performed using bisulfite sequencing, as described elsewhere (12). The candidate functional methylation markers were then validated by pyrosequencing in 25 independent FFPE MTC samples (8 RET^M918T, 4 RET^C634, 3 RAS, and 10 WT tumors). The region of interest, including the CpG assessed by the Infinium array, was first amplified from the bisulfite-treated DNA by a set of primers designed using the PyroMark assay design software (version 2.0.01.15). Following PCR amplification, pyrosequencing was performed using PyroMark Q24 reagents, vacuum prep workstation, equipment, and software (Qiagen). The primers used in these steps are listed in Supplementary Table S3.

Three paraffin-embedded TMAs and 7 complete sections were used for the evaluation of a JAK/Stat pathway effector pSTAT3 (Tyr705). Monoclonal antibody was obtained from Cell Signaling Technology (#9145). Two independent experienced pathologists (X. Matias-Guiu and M. Santacana) evaluated the staining by visual examination under a microscope. Since each TMA included 2 different tumor cores from each case, IHC scoring was done after examining both samples. When the 2 cores gave inconsistent results, the complete section of the corresponding tumor was evaluated. Samples were scored from 0 (no immunostaining) to 300 (intensive immunostaining) depending on both the intensity and extension of staining.

For the study of proliferation inhibition by TKI and the inhibition of the phosphorylation of STAT3, vandetanib (Seleckchem) and LLL12 (EMD Millipore; ref. 29) were made up as a stock solution of 5 mmol/L in DMSO. MZ-CRC-1 and TT cells were plated in 200 μL of medium at concentrations of 2 × 10⁴ cells per well. After 36 hours, cells were treated with increasing concentrations of the indicated compounds (0, 0.005, 0.05, 0.1, 0. 5, 1, 2.5, 5 μmol/L). A concentration of 0.1% DMSO was used in all conditions. This treatment was repeated every 4 days. Proliferation was measured at 10 days using a cell proliferation kit (WST-1 reagent, Roche). The concentration that led to 50% growth inhibition (IC₅₀) was calculated using GraphPad Prim. All experiments were performed in triplicate. For the combinational treatment, we also used the above-mentioned experimental scheme. A fixed concentration of 0.04 μmol/L of LLL12 was combined with increasing concentrations of vandetanib (0, 0.005, 0.05, 0.1, 0.5, 1, 2.5, 5 μmol/L). The suboptimal LLL12 concentration (0.04 μmol/L) used causes only 10% growth inhibition of MZ-CRC-1 cells. In addition, for Western blot analysis, MZ-CRC-1 cells were treated with increasing concentrations of vandetanib (0.01, 0.1, 0.5 μmol/L), 0.04 μmol/L LLL12, and a combination of LLL12 and vandetanib (0.04 and 0.1 μmol/L, respectively) for 24 hours. Cells were then washed with ice-cold PBS, collected by scraping, and cell lysates were prepared using the above-mentioned protocol. The following primary antibodies were used for immunoblotting: p-ERK, pAKT, pSTAT3, and GAPDH (Cell Signaling; #9101, #4060, #9145, and #2118, respectively). Absolute band intensities of the indicated proteins were captured and quantified with the Image Lab software v4.1 (BioRad). GAPDH was used as loading control normalizer. Data were expressed in relative units using the control condition as reference.

Among the tumors whose genome-wide DNA methylation levels were measured, 36 (75%) harbored a RET mutation (13 carried the RET^M918T mutation, 14 RET^C634, and 9 less frequent RET mutations) and 6 (12.5%) tumors harbored a RAS mutation. The remaining 6 tumors (12.5%) did not carry any alteration in the studied genes and were considered WT. The results of the mutational screening together with the main clinicohistopathologic characteristics of the samples are summarized in Supplementary Table S4.

We identified 11,915 probes that were constitutively unmethylated (β < 0.2) in all samples. According to DAVID functional enrichment analysis (30), the majority (97.5%) of these were located within CpG islands adjacent to housekeeping genes (Benjamini–Hochberg adjusted P = 7.1 × 10⁻¹⁰). On the other hand, 179 probes were constitutively methylated (β-values > 0.8 for all samples). Both constitutively unmethylated and methylated probes were excluded from further analyses. A PCA of the remaining 14,346 probes (belonging to 9,216 consensus coding sequences) did not identify any batch effects affecting the data (Supplementary Fig. S2).

Unsupervised hierarchical analysis using the 851 probes with the variance >0.2 identified 2 principal clusters that differed in their underlying genetics (Fig. 1). Only 2 samples did not fall into either cluster. Interestingly, cluster B showed higher levels of methylation when compared with cluster A (P = 3.0 × 10⁻⁹) and was enriched with WT cases and those harboring the RET^C634 mutation, whereas cluster A was composed mostly of RET^M918T-positive tumors (P < 0.02). We did not observe any clear clustering of the RAS-mutated samples.

Supervised analysis according to driver mutation allowed us to identify a list of differentially methylated probes associated with specific mutations, which was especially long for RET^M918T-related tumors. These findings confirmed those from the unsupervised clustering, as in this group, there were more hypomethylated probes (Fig. 2A and B and Supplementary Table S5). The protein-coding gene DKK4 (Dickkopf WNT signaling pathway inhibitor 4), previously described to be upregulated in RET^M918T-related MTCs (18), was one of the genes affected by hypomethylation among RET^M918T-positive carcinomas. Also consistent with the unsupervised clustering result, the lists of differentially methylated probes characteristic of the other genetic conditions (RET^C634- and RAS-related tumors) when comparing to WT, were considerably shorter (Fig. 2B). Nevertheless, hypomethylation of GAL was detected among RET^C634-positive tumors when comparing with WT, and this event could cause the increased expression of this molecule reported elsewhere (18).

Hypermethylation more frequently affected probes located within CpG islands (P = 0.0014) as well as those located near stem cell PolyComb Group target genes (P < 0.0001). Hypomethylated probes, on the other hand, were enriched with CpGs that are heavily methylated in embryonic stem cells (P < 0.0001; Fig. 2C).

The above findings indicated that some epigenetic changes can have a functional consequence, which encouraged us to systematically explore genes (protein and miRNA coding) regulated by methylation at their promoters (see Supplementary Fig. S1 for further details). mRNA expression data were available for more than 91% of protein-coding genes that were differentially methylated on the basis of tumor genotype. A negative correlation between expression and methylation was observed for 5% of these (Supplementary Fig. S3). PLCB2, DKK4, and MMP20 were selected for further study, as they have been previously described to have a role in the promoting tumor growth (18, 31–33). Moreover, a genome-wide assessment of methylation and miRNA expression of matching samples identified 5 miRNA genes (6.9%) showing negative correlation (Supplementary Fig. S3). MiR-10a was also found as differentially methylated in the previous supervised analyses and was selected for further studies. Two additional miRNA genes (miR-30a and -200c) showing negative correlation between expression and methylation, and a tendency toward differential methylation (FDR < 0.15) in supervised clustering, were also studied further due to their biologic function (21, 34, 35). As assessed by bisulfite sequencing, the levels of methylation of the 6 genes (3 protein-coding and 3 microRNAs) differed between the 2 available MTC cell lines (Fig. 3A). Apart from 2 genes that were not expressed by either of the cell lines (MMP20 and miR-30a), the expression of the other 4 genes confirmed the results of in silico predictions of inverse correlation between DNA methylation and gene expression (Fig. 3B). Because MZ-CRC-1 showed systematically higher levels of DNA methylation of these 6 loci, we treated it with the demethylating agent 5-aza-2′-deoxycytidine. A higher expression of all genes in the MZ-CRC-1 cell line was detected in the treated condition as compared with untreated one (Fig. 3C), further confirming the functional effect of these epigenetic changes.

CpGs from the promoters of 3 differentially methylated genes (DKK4, MMP20, PLCB2) were selected for validation. These genes also showed a negative correlation between DNA methylation levels and gene expression (Supplementary Fig. S3). As depicted in Fig. 4, the results of bisulfite sequencing in a subset of the original series revealed a high correlation with the array-based measures (R² ranging from 0.66 to 0.80). Moreover, we were able to replicate the findings by bisulfite pyrosequencing of an independent series of 25 FFPE MTCs (Fig. 4), although the difference between the DNA methylation levels of the compared groups was smaller, probably because the DNA source was paraffin tissue.

Functional annotation clustering by DAVID (30) of the genes with aberrant methylation (Supplementary Fig. S4) and by Gene Set Enrichment Analysis (GSEA) of differentially expressed genes (18) among RET^M918T-positive carcinomas uncovered a significant enrichment for signaling pathways such as cytokine–cytokine receptor interaction (P = 1.5 × 10⁻⁷) and JAK/Stat (P = 2.3 × 10⁻⁴). STAT3 is a JAK/Stat pathway effector and a putative downstream target of RET (36). Thus, we assessed by IHC its active phosphorylated form (at Tyr705) in 101 MTC samples. RET^M918T-related tumors were more frequently and intensively positively stained than the other genetic classes (P < 0.02, Fig. 5A) confirming activation of the JAK/Stat cascade and the in silico prediction. We aimed to exploit these differences to study in vitro the functional relevance of STAT3 activation in MTC. As revealed by proliferation assay, MZ-CRC-1 cells were more sensitive to inhibition of STAT3 phosphorylation than TT cells [IC₅₀ (MZ-CRC-1) = 0.06 μmol/L, IC₅₀ (TT) = 0.23 μmol/L, Fig. 5B]. In fact, the RET^M918T-mutated cell line showed a higher sensitivity to this treatment than to the standard of care, TKI vandetanib [IC₅₀ (MZ-CRC-1) = 0.62 μmol/L, Fig. 5B].

Because STAT3 is a downstream point of convergence for many growth factor pathways, synergy may be expected when combined with other inhibitors (37). A suboptimal dose of pSTAT3 inhibitor, when combined with vandetanib, apparently showed a synergistic effect on cell growth inhibition of RET^M918T-positive cells (Fig. 5C). Moreover, treatment with vandetanib produced a dose-dependent inhibition of pSTAT3 in MZ-CRC-1 and suboptimal combinational treatment of the 2 drugs (0.04 μmol/L LLL12 with 0.1 μmol/L vandetanib) potentiated the inhibitory effect on both canonical (pERK, pAKT) and noncanonical (pSTAT3) effectors of RET downstream signaling (Fig. 5D).

The current genomic era offers great opportunities to answer many clinically relevant questions concerning various types of cancer through comprehensive and integrative analyses of OMIC data. Some cancer types are being neglected by these efforts, generally due to their low prevalence, and thus small overall public health impact. MTC is one of the least prevalent subtypes of thyroid cancer but is responsible for a large proportion of thyroid cancer–related deaths (38), mainly because treatment with cytotoxic drugs and/or standard radiotherapy has been proven to be ineffective (39).

Using the largest cohort of MTC samples reported to date, we have applied a multilayer OMIC data integration focused on the effects of DNA methylation in the etiology of this tumor. The current work is the first to report not only that underlying genetics relates closely to genome-wide DNA methylation fingerprints but also that the aberrant methylation events affect specific genomic loci (e.g., PolyComb targets) pointing to the existence of an epigenetic progenitor cell signature in MTC (40). Furthermore, we identify both protein-coding and microRNA genes, the expression of which is negatively correlated with the methylation status of their promoters, and we confirm some of these findings in MTC cell lines. We propose promoter methylation of MMP20, DKK4, and PLCB2 as functional epigenetic markers of MTC genetic classes. The integrative approach allowed us to uncover the involvement of the JAK/Stat pathway in RET^M918T-related MTCs, which apparently sustains the proliferative features of the MTC cell model for this genetic class.

DNA methylation is traditionally believed to be a regulator of gene expression. In the case of protein-coding genes, the classical view of the negative effect of methylation on expression is being challenged by recent technical advances. Many cancer-orientated studies, including the current one, report a surprisingly low percentage of genes for which this inverse relationship applies (26). The situation with microRNAs is even more complex, given the difficulties with identifying their promoter regions (41). Recently, a methodology became available that recognizes putative promoter regions of miRNAs using both sequence- and histone-based techniques (27). By integrating this strategy in our analysis, we were able to identify CpGs belonging to additional putative miRNA promoters included on the 27K array. The overall portion of miRNA genes for which a negative correlation between methylation and expression was observed was similar to that for protein-coding genes. Nevertheless, for the majority of the genes selected for validation of the inverse correlation between expression and methylation in MTC cell lines, we were able to functionally confirm the in silico predictions. One of these genes, DKK4, has previously been linked to MTC pathogenesis (18), and here we propose aberrant methylation as the regulatory mechanism responsible for its overexpression in RET^M918T-related tumors. On the other hand, PLCB2 plays a role in multiple transmembrane signal transduction pathways involving inositol lipid metabolism. In cancer, its overexpression has been associated with mitosis promotion, migration, and poor outcome (31, 42). Our finding of hypomethylation in the RET^M918T-positive tumor samples, which was negatively associated with increased PLCB2 expression, is consistent with this gene having a cancer-promoting role. Even though the prognostic utility of MMP20, a member of the family of metalloproteinases, has already been demonstrated in cancer (32), we did not detect its expression in the studied cell lines. Nevertheless, its validation in an independent series of MTC carcinomas confirms its potential utility as a methylation marker of RET^M918T-related tumors.

Regarding miRNA genes, the expression of miR-200c has been recently described as causative of the metastatic potential of human MTCs (21). Interestingly, according to our results in MTC cell lines, the aberrant methylation of the putative promoter of miR-200c could underlie its differential expression among patients with distinct disease outcomes. MiR-10a mostly exerts an oncogenic effect (43, 44) and we found little DNA methylation of its promoter in both MTC cell lines. There is evidence from hepatocellular carcinoma that miR-10a is negatively regulated by DNA methylation (45), which was further validated by our findings. MiR-30a warrants further study due to its tumor-suppressive features (34, 35). PROM1, a cancer stem cell marker found overexpressed in RET^M918T-positive tumors (18), is one of its predicted targets, according to targetscan tool (http://www.targetscan.org/). Our findings indicate that a hypermethylation-driven downregulation of miR-30a could lead to an upregulation of its target PROM1.

Our integrated analysis pointed to the involvement of signaling networks such as cytokine–cytokine interaction and JAK/Stat in RET^M918T tumors. These pathways have previously been related to malignant tumor behavior, which is consistent with the observation that these patients with MTC have the worst clinical outcomes (5). The genetic class–specific implication of JAK/Stat pathways was confirmed by IHC in a large number of human tumors. This is in line with previous in vitro findings from 293T cells transfected with WT RET or RET^M918T, where only the latter induced constitutive phosphorylation of STAT3. On top of that, RET^M918T was far more efficient in recruiting STAT3 protein than both RET^C634 and wild-type RET in vitro (46). Constitutive activation of STAT3 has been reported at high frequency in large numbers of malignant cell lines, in vivo animal models, and human tumors (reviewed in ref. 37). In general, activation of STAT3 is associated with a worse prognosis due to its role in uncontrollable proliferation and modulation of microenvironment (reviewed in ref. 47). Thus, STAT3 is a very attractive molecule for targeted therapies (reviewed in ref. 48). Indeed, while no drug has shown promising results in clinical trials (49), it has been suggested that inhibition of STAT3 signaling in combination with other therapies (targeted therapies or radiotherapy) could lead to clinical benefit (50). We found that low doses of the inhibitor of pSTAT3 increased the sensitivity of MTC cells to vandetanib treatment. This could have an important impact on the management of patients with MTC, as vandetanib treatment leads to frequent toxicities, which require dose reductions or even discontinuation (8).

In summary, this comprehensive integrative genomic study performed using the largest cohort of MTC samples reported to date, provides insights into the involvement of DNA methylation in the etiology of this disease. We confirmed the regulatory role of DNA methylation for DKK4, PLCB2, MMP20, miR-10a, miR-30a, and miR-200c using MZ-CRC-1 and TT cell lines. Moreover, hypomethylation may induce activation of key pathways related to the malignant behavior of RET^M918T-related MTCs. Inhibition of one of these—the STAT3 signaling pathway—sensitizes RET^M918T-mutated MTC cells to standard vandetanib treatment.

No potential conflicts of interest were disclosed.

Conception and design: V. Mancikova, C. Montero-Conde, M. Robledo

Development of methodology: V. Mancikova, C. Montero-Conde, J. Perales-Paton, E. Castelblanco, X. Matias-Guiu, M. Robledo

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): V. Mancikova, C. Montero-Conde, K. Jodkowska, E. Castelblanco, S. Borrego, M. Encinas, X. Matias-Guiu, M. Robledo

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): V. Mancikova, C. Montero-Conde, J. Perales-Paton, A. Fernandez, M. Santacana, L. Inglada-Pérez, X. Matias-Guiu, M. Fraga, M. Robledo

Writing, review, and/or revision of the manuscript: V. Mancikova, C. Montero-Conde, K. Jodkowska, L. Inglada-Pwska, X. Matias-Guiu, M. Fraga, M. Robledo

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): V. Mancikova, X. Matias-Guiu, M. Robledo

Study supervision: M. Robledo

We thank Manuel Morente and María Jesús Artiga from the Spanish National Tumor Bank Network (CNIO) for their hard work collecting the tumor samples used in this study. We also thank Magdolna Djurec and Teresa Blasco for their help with the Western blot analysis. The genotyping service was carried out at CEGEN-PRB2-ISCIII; it is supported by grant PT13/0001, ISCIII-SGEFI/FEDER.

This work was supported by the Fondo de Investigaciones Sanitarias (FIS) project PI14/00240, cofinanced by FEDER, and the Comunidad de Madrid (Grant S2011/BMD-2328 TIRONET) to M. Robledo. L. Inglada-Pérez is supported by the Centro de Investigacion Biomédica en Red de Enfermedades Raras (CIBERER). V. Mancikova was supported by a predoctoral fellowship from the "la Caixa"/CNIO international PhD programme. C. Montero-Conde is supported by a postdoctoral fellowship from the Fundación AECC.

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