DNA methylation markers for determining the risk of metastatic progression of thyroid cancer

Abstract

The present invention generally refers to a method of determining the methylation level of a set of biomarkers in a subject suspected of having thyroid cancer, the method comprising: a) from an isolated or extracted genomic DNA obtained from a biological sample from the subject suspected of having thyroid cancer, determining the methylation level of one or more biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462,137-85,462,137), cg18305394 (located in the 5' region of PCK1 gene, chr20:56,134,788-56,134,788), cg11925561 (located in an intergenic region of chromosome 5, chr5:3,304,588- 3,304,588), and cg14672100 (located in an intergenic region of chromosome 5, chr5:6,344,259-6,344,259) from the extracted genomic DNA. The present invention also refers to diagnostic, prognostic and patient selection uses of said method.

Description

[0001] DNA METHYLATION MARKERS FOR DETERMINING THE RISK OF METASTATIC PROGRESSION OF THYROID CANCER

[0002] Field of the invention

[0003] The present invention relates to conditions characterized by differentially methylated CpG dinucleotides and, in particular, to diagnostic and prognostic methods that exploit the presence of such CpGs that exhibit altered DNA methylation levels.

[0004] Background art

[0005] Thyroid cancer is the most common endocrine neoplasia (Siegel RL, et al. Cancer statistics, 2024. CA Cancer J Clin. 2024; 74(1):12-49) and its incidence has been steadily increasing over the last few decades (Lim H, et al. Trends in Thyroid Cancer Incidence and Mortality in the United States, 1974- 2013. JAMA. 2017;317(13):1338; Kim J, et al. Geographic influences in the global rise of thyroid cancer. Nat Rev Endocrinol.

[0006] 2020;16(1):17-29). Most patients with differentiated thyroid cancer (DTC), including both papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC), have an excellent prognosis and can be efficiently cured by standard treatments based on surgery, radioiodine (RAI) and thyroid hormone suppression therapy (Haugen BR, et al. 2015 American Thyroid Association Management Guidelines for Adult Patients with Thyroid Nodules and Differentiated Thyroid Cancer: The American Thyroid Association Guidelines Task Force on Thyroid Nodules and Differentiated Thyroid Cancer Thyroid.

[0007] 2016;26(1):1-133). However, a subgroup of them progresses and develops distant metastases (DM), which represent the main cause of thyroid cancer-related death (Schlumberger M, Leboulleux S. Treatment of distant metastases from follicular cell-derived thyroid cancer. F1000Prime Rep. 2015;7:22.; Vallejo Casas JA, et al. Initial clinical and treatment patterns of advanced differentiated thyroid cancer: ERUDIT study.Eur Thyroid J. 2022;11(5):e210111). Specifically, 1-4% of DTC patients present DM at initial diagnosis (synchronous DM), while 7-23% of patients develop metachronous DM, even many years after initial diagnosis and treatment (Wang LY, et al. Multi-Organ Distant Metastases Confer Worse Disease-Specific Survival in Differentiated Thyroid Cancer. Thyroid. 2014;24(11): 1594-1599). Early identification of patients with DTC at high risk of developing DM would allow the clinician to offer these patients a more aggressive initial treatment as well as an intensive follow-up to prevent rapid progression and increase survival.

[0008] Ultrasound (US) and fine-needle aspiration biopsy (FNAB) are the gold standard in the evaluation of thyroid nodules. Several pathological factors have been associated with DM in DTC such as tumor size, vascular invasion and lateral nodal metastasis; however, their impact in the development of DM over time remains controversial (Kim H, et al. Prediction of follicular thyroid carcinoma associated with distant metastasis in the preoperative and postoperative model. Head Neck. 2019;41(8):2507-2513; Wang LY, et al. Multi-Organ Distant Metastases Confer Worse Disease-Specific Survival in Differentiated Thyroid Cancer. Thyroid. 2014;24(11):1594-1599; Chen D, et al. Innovative analysis of distant metastasis in differentiated thyroid cancer. Oncol Lett. 2020 Mar; 19(3): 1985-1992; Vuong HG, et al. Role of molecular markers to predict distant metastasis in papillary thyroid carcinoma: Promising value of TERT promoter mutations and insignificant role of BRAF mutations — a meta-analysis. Tumor Biol.

[0009] 2017;39(10): 101042831771391; Yip L, et al. Risk assessment for distant metastasis in differentiated thyroid cancer using molecular profiling: A matched case-control study. Cancer. 2021;127(11):1779-1787). Currently, the most widely used marker to detect DM is the measurement of serum thyroglobulin (Tg), but it has some limitations such as the low capacity to discriminate between benign and malignant nodules, and loss of Tg expression in some aggressive DTC cases. Imaging methods other than US (RAI wholebody scintigraphy, CT, MRI, PET / CT, etc.) are also of great value in tumor andmetastasis detection, but are not predictive. Thus, the identification of predictors of DM is of utmost importance to improve the therapeutic intervention success rates.

[0010] In recent years, genomic studies have been carried out in DTC (Yoo S-K, et al. Integrative analysis of genomic and transcriptomic characteristics associated with progression of aggressive thyroid cancer. Nat Commun. 2019;10(1):2764), especially in PTC, which have led to important advances in understanding the molecular landscape of low- and intermediate-risk DTC. However, high-risk DTC, including those developing DM, are still poorly understood. DTC has one of the lowest mutation rates of all solid tumors, and most of the known driver mutations affect MAPK and PI3K pathways. BRAF(V600E) is the most common mutation (50-60% PTC) followed by mutations in RAS (20-40% FTC and 15-20% follicular variant of PTC -fvPTC) (Song YS, Park YJ. Genomic Characterization of Differentiated Thyroid Carcinoma. Endocrinol Metab.

[0011] 2019;34(1):1; Acuña-Ruiz A, et al. Genomic and epigenomic profile of thyroid cancer. Best Pract Res Clin Endocrinol Metab. 2023;37(1):101656). Other less common mutations (TERT promoter, PIK3CA, PTEN, IDH1, TP53 or AKT1) have been associated with more aggressive tumors and more advanced stages, but the prognostic value of these genetic alterations in DTC is yet to be fully elucidated (Marotta V, et al. Application of molecular biology of differentiated thyroid cancer for clinical prognostication. Endocr Relat Cancer. 2016;23(11): R499-R515). The currently commercial genetic tools based on next-generation sequencing (NGS) (e.g. Afirma XA, Thyroseq v3, ThyGeNext / ThyraMIR) that have been designed to improve the diagnostic capacity of thyroid nodule molecular testing have also introduced opportunities to provide prognostic information. However, their ability to identify patients that will progress has yet not been determined (Patel J, Klopper J, Cottrill EE. Molecular diagnostics in the evaluation of thyroid nodules: Current use and prospective opportunities. Front Endocrinol (Lausanne). 2023; 14:1101410. Published 2023 Feb 24. doi:10.3389 / fendo.2023.1101410). In addition, these methodologies have a high costwhich implies an important economic impact on national healthcare systems. Therefore, there is an urgent need to develop new cost-effective tools to improve patient prognosis that can be easily incorporated into the clinical setting.

[0012] DNA methylation is a well-known epigenetic mechanism involved in cancer initiation and progression (Ehrlich M. DNA hypomethylation in cancer cells. Epigenomics.

[0013] 2009;1(2):239-259; Esteller M. Epigenetics in cancer. N Engl J Med.

[0014] 2008;358(11):1148-1159). Changes in this epigenetic mark have rapidly gained clinical attention as biomarkers for diagnostic, prognostic and predictive purposes in many cancers (Davalos V, Esteller M. Cancer epigenetics in clinical practice. CA Cancer J Clin.

[0015] 2023;73(4):376-424. doi:10.3322 / caac.21765). DTC is also characterized by DNA methylation alterations (Zafon C, et al. DNA methylation in thyroid cancer. Endocr Relat Cancer. 2019;26(7): R415-R439). DNA methylation mostly occurs in the cytosines within the CpG dinucleotides and is involved in many biological processes such as regulation of gene expression and DNA structure (Mattei AL, et al. DNA methylation: a historical perspective. Trends in Genetics. 2022;38(7):676-707; Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet.

[0016] 2012;13(7):484-92). Most of the human genome is methylated, except some CpG-rich regions - known as CpG islands (CGIs) -which usually remain unmethylated and are associated with promoters. DNA methylation is often considered a repressive mark. However, its function varies with the genomic context: in regulatory elements (promoters and enhancers) and repetitive sequences, DNA methylation represses transcription, changing the degree of compaction of chromatin; in the body of genes, it can enhance transcriptional elongation and affect splicing (Mattei AL, et al. DNA methylation: a historical perspective. Trends in Genetics. 2022;38(7):676-707; Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet.

[0017] 2012;13(7):484-92). In comparison with normal cells, tumor cells harbor both gains and losses of DNA methylation (hypermethylation and hypomethylation, respectively)(Ehrlich M. DNA hypomethylation in cancer cells. Epigenomics. 2009;1(2):239-259; Esteller M. Epigenetics in cancer. N Engl J Med. 2008;358(11):1148-1159). Hypermethylation is usually local and affects regulatory elements, while hypomethylation has been detected both locally and globally (affecting extensive domains of the genome that include repetitive and unique sequences). All these epigenetic alterations have been widely associated with downstream transcriptional alterations that trigger a profound transformation of the cellular phenotype.

[0018] Several genome-wide DNA methylation studies in DTC, including one by our group (Mancikova V, et al. DNA methylation profiling of well-differentiated thyroid cancer uncovers markers of recurrence free survival. Int J Cancer. 2014;135(3)), have revealed that DTC methylomes are specifically associated with histology and BRAF and RAS mutations (Beltrami CM, et al. Integrated data analysis reveals potential drivers and pathways disrupted by DNA methylation in papillary thyroid carcinomas. Clin Epigenetics. 2017;9(1):45. 753; Ellis RJ, et al. Genome-wide methylation patterns in papillary thyroid cancer are distinct based on histological subtype and tumor genotype. J Clin Endocrinol Metab. 2014;99(2): E329-37; Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell.

[0019] 2014;159(3):676-690; Marczyk VR, et al. Classification of Thyroid Tumors Based on DNA Methylation Patterns. Thyroid. 2023;33(9): 1090-1099). However, to the best of our knowledge, no studies have characterized the methylome of primary DTC tumors that develop DM (metastatic primary DTC [mDTC]). In a previous study by our group exploring the role of DNA methylation, we found that mDTC, together with poorly differentiated thyroid cancer (PDTC) and anaplastic thyroid cancer (ATC) - the most aggressive thyroid tumors - were affected by global hypomethylation of Alu repeats (Klein Hesselink EN, et al. Increased Global DNA Hypomethylation in Distant Metastatic and Dedifferentiated Thyroid Cancer. J Clin Endocrinol Metab. 2018;103(2):397-406). These results not only suggest the involvement of this epigenetic alteration in thyroidcancer progression and dedifferentiation but also its potential as prognostic marker. In fact, the 2022 World Health Organization (WHO) Classification of Thyroid Neoplasms proposes global hypomethylation as a marker of adverse biology in DTC (Baloch ZW, et al. Overview of the 2022 WHO Classification of Thyroid Neoplasms. Endocr Pathol.

[0020] 2022;33:27-63).

[0021] The rapidly increasing incidence of DTC, mainly due to overdiagnosis of low-risk DTC, has led to a shift towards conservative treatment approaches (active surveillance, deferred intervention, minimalist surgery and / or no radioiodine — RAI). However, these approaches can lead to the undertreatment of high-risk patients, increasing the risk of poor outcomes and future reinterventions. In fact, many clinicians are concerned about the risk of DM, and to prevent progression they tend to overtreat patients (extensive surgery, high RAI doses) and perform intense follow-ups (ultrasound, scans and analytics), which may have adverse effects, including psychological stress. Therefore, there is an urgent need for quantifiable biomarkers that predict the evolution of each DTC patient.

[0022] Since the role of DNA methylation in progressive thyroid cancer remains unclear, we aimed to study DNA methylomes along DTC metastatic progression to identify new biomarkers that predict the development of DM.

[0023] It is expected to develop an easy-to-implement epigenetic tool that will enhance preoperative risk stratification. It will identify DTC patients at high risk of DM allowing for more personalized treatment and follow-up, while helping to reduce the overtreatment of patients at low risk of DM. In addition, its incorporation into clinical practice will reduce healthcare costs.Brief description of the invention

[0024] In summary, the present invention refers to an in vitro method for determining the DNA methylation level of a set of biomarkers in a subject having thyroid cancer, the method comprising:

[0025] a. from an isolated or extracted genomic DNA obtained from a biological sample from the subject, determining the methylation level of one or more biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), cg11925561 (located in an intergenic region of chromosomes, chr5:3, 304, 588-3, 304, 588), and cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259) from the isolated or extracted genomic DNA.

[0026] In some embodiments, the step a) comprises determining the methylation level of at least cg14672100 (located in an intergenic region of chromosome 5, chr5:6,344,259-6,344,259).

[0027] In some embodiments, the step a) comprises determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259.

[0028] In some embodiments, the step a) further comprises determining the methylation level of at least one or more biomarkers selected from the group consisting of cg25901805 (located in PLAGL1 gene, chr6:144, 286, 162-144, 286, 162), cg10780897 (located in RP1 gene, chr8:55, 528, 099-55, 528, 099), cg07530798 (located in chr1:119,535,589-119,535,589) and cg04690464 (located in LINC02691 gene, chr14: 104, 711,952-104,711,952).In some embodiments, step a) comprises determining the methylation level of at least all biomarkers selected from the group consisting of cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259), cg25901805 (located in PLAGL1 gene, chr6:144, 286, 162-144, 286, 162), cg10780897 (located in RP1 gene, chr8:55, 528, 099-55, 528, 099), cg07530798 (located in chr1:119,535,589-119,535,589) and cg04690464 (located in LINC02691 gene, chr14:104, 711,952-104,711,952).

[0029] In some embodiments, the methylation level is determined by a quantitative assay based on bisulfite pyrosequencing to measure the DNA methylation level.

[0030] Additionally, the present invention provides an in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:

[0031] a. from an isolated or extracted genomic DNA obtained from a biological sample from the subject, determining the DNA methylation level of a set of biomarkers as defined in step a) of the previous embodiments; and

[0032] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by a classification method that identifies the subject as a subject at risk of developing metastasis based on the biomarkers identified in step (a), preferably by a predictive model which correlates the methylation level of one or more of the biomarkers identified in step (a) with representative methylation level profiles from samples obtained from subjects having thyroid cancer previously identified as suffering from metastasis.

[0033] In some embodiments, step a) comprises determining the DNA methylation level of a set of biomarkers as defined in step a) of the previous embodiments.

[0034] In some embodiments, the methylation level is determined by a quantitative assay selected from the group consisting of bisulfite-based methods, methylation-sensitive restriction enzyme-based methods, microarrays-based methods, next-generationsequencing (NGS) methods, methylation-specific melting curve analysis, mass spectrometry, immunoprecipitation methods, pyrosequencing, bisulfite pyrosequencing, digital polymerase chain reaction (digital PCR), quantitative polymerase chain reaction (qPCR) or any combination thereof to measure the DNA methylation level.

[0035] The present invention also relates to an in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:

[0036] a. from an isolated or extracted genomic DNA obtained from a biological sample from the subject, determining the DNA methylation level of a set of biomarkers as defined in step a) of the previous embodiments; and b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by a classification method that identifies the subject as a subject at risk of developing metastasis based on the biomarkers identified in step (a), preferably by a predictive model which correlates the methylation level of one or more of the biomarkers identified in step (a) with representative methylation level profiles from samples obtained from subjects having thyroid cancer previously identified as suffering from metastasis.

[0037] In some embodiments, step a) comprises determining the DNA methylation level of a set of biomarkers as defined in step a) of the previous embodiments.

[0038] In some embodiments, the biological sample is a sample containing DNA derived from thyroid tumor.

[0039] The present invention is also directed to an in vitro method for classifying a subject as a subject having thyroid cancer at risk of metastasis, preferably of distant metastasis (DM), the method comprising:

[0040] a. determining whether the patient is at risk of metastasis, preferably of distant metastasis (DM), following the method of the previous embodiments; andb. classifying said patient as a subject having thyroid cancer at risk of metastasis, preferably of distant metastasis (DM), if the patient is identified as a subject at risk of metastasis, preferably of distant metastasis (DM), according to step a).

[0041] This invention further provides a composition suitable for the treatment of a subject having thyroid cancer at risk of metastasis, for use in the treatment of thyroid cancer in a subject classified as having thyroid cancer at risk of metastasis according to step b) of the previous in vitro method for classifying a subject.

[0042] In some embodiments, the composition comprises one or more therapies selected from the group consisting of radiation therapies, hormonal therapies, targeted therapies and chemotherapies, or any combination thereof.

[0043] Brief description of the figures

[0044] Figure 1. Global DNA methylation changes during metastatic progression in differentiated thyroid cancer (DTC). (A) Global DNA methylation levels in normal tissue (NT), non-metastatic primary DTC (non-mDTC), metastatic primary DTC (mDTC), lymph node metastases (LNM) and distant metastases (DM), measured as the mean of the beta-values of the CpGs located at different CGI-related regions according to the distribution shown in (B). The number of samples in each group is indicated in parentheses on the X- axis. The Kruskal-Wallis test and Dunn’s post hoc test were used to determine statistical significance.

[0045] Figure 2. DNA methylation dynamics during metastatic progression in papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC). (A) Number of hypomethylated (green) and hypermethylated (red) CpGs in non-metastatic (non-m) and metastatic (m) primary tumors, and distant metastases (DM) compared to normal tissue. CpGs are shaded with a color scale according to their specificity or overlap in the different risk categories. (B) Overlap of differentially methylated CpGs in non-m PTC and non-mFTC, mPTC and mFTC and in distant metastases (DM) derived from PTC (PTC-DM) and FTC (FTC-DM). (C) Relative distribution of the differentially hypomethylated (green) and hypermethylated (red) CpGs across different CGI-related regions (upper panel) and gene-related regions (lower panel).

[0046] Figure 3. Chromatin dynamics of the differentially methylated CpGs during metastatic progression in papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC). Enrichment of chromatin states (based on the Roadmap 18-state model and 18 cell types) of the hypomethylated and hypermethylated CpGs in non-metastatic (non-m) and metastatic (m) primary tumors and in distant metastases (DM) compared to normal tissue. Chromatin states are grouped as TSS-proximal regions (green), transcription-related states (cream color), enhancers (purple), ZFN genes and repeats (yellow), heterochromatin (dark grey), bivalent promoters and enhancers (blue), Polycomb-repressed regions (light grey), and quiescent region (white).

[0047] H1.hESC: human embryonic stem cell line; Skin.ker: foreskin keratinocyte; Skin.mel: foreskin melanocyte; Brain. hip: brain hippocampus middle; Brain. ctx: Brain dorsolateral prefrontal cortex; Pane. islet: pancreatic islets; GM12878: lymphoblastoid; HMEC: mammary epithelial; HSMM: skeletal muscle myoblasts; HUVEC: umbilical vein endotelial; K562: leukaemia; NHEK: epidermal keratinocyte; NHLF: lung fibroblast; TssA: active Transcription Start Site (TSS); TssFlnk: flanking TSS; TssFlnkU: flanking TSS upstream; TssFlnkD: flanking TSS downstream; Tx: strong transcription; TxWk: weak transcription; EnhG1: genic enhancer 1; EnhG2: genic enhancer 2; EnhA1: active enhancer 1; EnhA2: active enhancer 2; EnhWk: weak enhancer; ZNF / Rpts: ZNF genes and repeats; Het: heterochromatin; TssBiv: bivalent / poised TSS; EnhBiv: bivalent enhancer; ReprPC: repressed PolyComb; ReprPCWk: weak repressed PolyComb; Quies: Quiescent / low.

[0048] Figure 4. Gene ontology analysis of genes differentially methylated during metastatic progression of papillary (PTC) and follicular thyroid cancer (FTC). Dot chart showing the enriched GO terms, listed on the left according to the different GO categories, of the genes differentially methylated at their 5' regions in low-risk (LR) PTC and FTC, metastatic primary tumors (mPTC and mFTC), and distant metastases (DM), compared to normal tissue. Only GO terms shared between at least two groups of samples are represented. Other GO terms are listed in Supplemental Tables S9-S11).Colors indicate the different GO categories. Dot size is based on fold-change (FC).

[0049] Figure 5. Correlation between DNA methylation and gene expression. Pie charts showing (A) the percentage of genes differentially methylated (DMe) in non-metastatic papillary thyroid cancer (non-mPTC) compared to normal tissue from our discovery series (light yellow) and differentially expressed (DE) in low-risk PTC compared to normal tissue from TCGA dataset (dark yellow) and (B) the percentage of genes DMe in metastatic primary papillary cancer (mPTC) compared to normal tissue from our discovery series (light yellow) and DE in high-risk PTC compared to normal tissue from TCGA dataset (dark yellow). Stacked bar plots show the number of genes in each type of correlation between DNA methylation and gene expression.

[0050] Figure 6. Characterization of the 156 CpG-signature. (A) Unsupervised hierarchical cluster analysis of the 156 CpG-signature in our discovery series. Risk category, diagnosis of distant metastases (DM), mutations, and histology are shown in color-coded rows below the upper dendogram. (B) Global DNA methylation levels in normal tissue (NT), non-metastatic (non-m) and metastatic (m) primary differentiated thyroid cancer (DTC), lymph node metastases (LNM), and distant metastases (DM), measured as the mean of the beta-values of the 156-CpG signature. The number of samples in each group is indicated in brackets on the X- axis. The Kruskal-Wallis test and Dunn's post hoc test were used to determine statistical significance. (C) ROC curve for the capacity of the global DNA methylation levels of the156-CpG signature to predict the risk of developing DMs. (D) Distribution of the differentially hypo- and hypermethylated CpGs from the 156-CpG signature across different CGI-related regions, gene-related regions and enhancer-related regions. The number of hypo- and hypermethylated CpGs is indicated in brackets on the X- axis.

[0051] MUT: mutated; V600E: BRAFV600E mutation; other MUT: other BRAF mutations; AUC: area under the curve; EPIC: EPIC array used as background; Hypo: hypomethylated; Hyper: hypermethylated; Enh: enhancer.Figure 7. Characterization of the 4-CpG predictive model. (A) Mean beta-values of each of the four CpGs in normal tissue (NT), non-metastatic (non-m) and metasttic (m) primary differentiated thyroid cancer (mDTC), lymph node metastases (LMN) and distant metastases (DM). The Kruskal-Wallis test and Dunn's post hoc test were used to determine statistical significance. *P < 0.05; **P < 0.01; ***P < 0.001. The number of samples in each group is indicated in parentheses. (B) Unsupervised hierarchical cluster analysis of the four CpGs in our discovery series. (C) Results of the 4-CpG predictive model using beta-values from the EPIC array applied to samples from the discovery series and validation series I. The number of samples in each group is indicated in parentheses. (D) Results of the adjusted 4-CpG predictive model using DNA methylation values obtained by bisulfite pyrosequencing (BPS) applied to samples from the discovery series and validation series II. The number of samples in each group is indicated in parentheses. (E) ROC curve of the adjusted 4-CpG predictive model using DNA methylation values obtained by BPS applied to samples from validation series II.

[0052] AUC: area under the curve.

[0053] Figure 8. Model of differentiated thyroid cancer metastatic dissemination. Our results suggest a linear progression model of differentiated thyroid cancer (DTC), in which tumors evolve over time acquiring genetic and epigenetic alterations and giving rise to different subclones (represented here by different colors). Once a pro-metastatic subclone emerges and outgrows, becoming a major subclone (lightgolden), these cells can disseminate and seed distant metastases (DM). As a result, the paired primary tumor (mDTC) and DM present genetic and epigenetic similarities. The dissemination and seeding of DMs occur before the initial diagnosis (Dx) of the primary tumor; however, in many patients, the DM may be biologically and clinically dormant for years before progression (metachronous DM). Epigenetically, there is a progressive increase in alterations of DNA methylation from non-m DTC to DM, with a high prevalence of hypomethylation. Created with BioRender.com.

[0054] Hypo: hypomethylated; Hyper: hypermethylated.Figure 9. Global DNA methylation during metastatic progression of papillary (PTC) and follicular thyroid cancer (FTC). Global DNA methylation levels measured as the mean of the beta-values of the CpGs located at (A) CGIs, (B) CGI shores, and (C) nonCGI regions in normal tissue (NT), non-metastatic (non-m) and metastatic (m) primary tumors, lymph node metastases (LNM), and distant metastases (DM). The number of samples in each group is indicated in parentheses on the X-axis. The Kruskal-Wallis test and Dunn's post hoc test were used to determine statistical significance.

[0055] Figure 10. Global DNA methylation of repetitive and unique sequences during metastatic progression of papillary (PTC) and follicular thyroid cancer (FTC).

[0056] Global DNA methylation levels measured as the mean of the beta-values of the CpGs located at CGI shores (A and B) and non-CGI regions (C and D) depending on whether CpGs are located in repetitive or unique sequences. Methylation analysis was performed in normal tissue (NT), non-metastatic (non-m) and metastatic (m) primary tumors, lymph node metastases (LNM), and distant metastases (DM). The number of samples in each group is indicated in parentheses on the X- axis. The Kruskal-Wallis test and Dunn's post hoc test were used to determine statistical significance.

[0057] Figure 11. Study design. We profiled DNA methylation during metastatic progression in differentiated thyroid cancer (DTC) in the discovery series using EPIC arrays. We performed differential analyses, which allowed us to identify differentially methylated CpGs in the different tumor risk categories compared to normal tissue (NT). These CpGs were characterized through different analyses. In addition, we compared non-metastatic (non-m) and metastatic (m) primary tumors and identified a signature of 156 CpGs that distinguished non-mDTC from mDTC. The 156-CpG signature was reduced to four CpGs that were used to generate a predictive model for the risk of developing DM. Both the 156-CpG signature and the 4-CpG predictive model were validated in validation series I by EPIC array. In addition, the DNA methylation levels of the four CpGs were analyzed by bisulfite pyrosequencing (BPS) in a subset of samples from the discovery series, andthe BPS values were used to retrain the predictive model, generating an adjusted predictive model that was then validated in validation series II. Only primary DTCs were used to generate the predictive models, but NT and DM were also used to validate the model.

[0058] Figure 12. Alterations in DNA methylation during metastatic progression of differential thyroid cancer (DTC). (A) Number of hypomethylated (green) and hypermethylated (red) CpGs in non-metastatic (non-m) and metastatic (m) primary tumors, and distant metastases (DM). (B) Number of hypomethylated (green) and hypermethylated (red) CpGs in non-metastatic (non-m) and metastatic (m) papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC), lymph node metastases (LNM) and DM.

[0059] Figure 13. Correlation between the beta values of normal tissue (NT) and papillary thyroid cancer (PTC) samples from our discovery series (X-axis) and TCGA series (Y-axis). Of the differentially methylated (DMe) CpGs identified in our discovery series, 2001 in non-metastatic primary PTC (non-mPTC) and 2596 in metastatic primary PTC (mPTC) were available in TCGA database. (A) Beta values of the 2001 CpGs in NT (upper panel) and non-mPTC (lower panel). (B) Beta values of the 2596 CpGs in NT (upper panel) and mPTC (lower panel). Correlations were evaluated by Pearson’s correlation coefficient. Additionally, 909 out of the 2001 CpGs (45.43%) and 1419 out of the 2596 CpGs (54.7%) were differentially methylated in TCGA low-risk and high-risk PTC compared to NT.

[0060] Figure 14. Comparison of DNA methylation levels between tumor and normal tissue (NT). Density scatter plots showing the distribution pattern of DNA methylation (beta-values) for non-metastatic primary papillary thyroid cancer (non-mPTC), metastatic primary PTC (mPTC), distant metastases derived from PTC (PTC-DM), non-metastatic primary follicular thyroid cancer (non-mFTC), metastatic primary FTC (mFTC), and DMderived from FTC (FTC-DM), in comparison with NT. The number of samples is indicated in parentheses. Density is color-coded from low (blue) to high (red).

[0061] Figure 15. Overlap of differentially methylated CpGs during metastatic progression of papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC).

[0062] Venn diagrams representing the overlap of differentially methylated CpGs compared to normal tissue (NT). (A) Overlap between non-metastatic primary PTC (non-mPTC) and metastatic primary PTC (mPTC) and between mPTC and distant metastases (DM) derived from PTC (PTC-DM). (B) Overlap between non-metastatic primary FTC (non-mFTC) and metastatic primary FTC tumor (mFTC) and between mFTC and DM derived from FTC (FTC-DM). Green circles indicate hypomethylated CpGs; red circles indicate hypermethylated CpGs.

[0063] Figure 16. Chromatin dynamics of the CpGs that were differentially methylated during metastatic progression of papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC) in comparison with normal tissue (NT). Enrichment of chromatin states (based on the Roadmap 18-state model and 18 cell types) of the hypo-and hypermethylated CpGs located at 5' gene regions, gene bodies and intergenic regions in non-metastatic (non-m) and metastatic (m) primary PTC and FTC and distant metastases (DM) derived from PTC (PTC-DM), and FTC (FTC-DM). Chromatin states are grouped as TSS-proximal regions (green), transcription-related states (cream), enhancers (purple), ZFN genes and repeats (yellow), heterochromatin (dark grey), bivalent promoters and enhancers (blue), Polycomb-repressed regions (light grey), and quiescent region (white).

[0064] Hl.hESC: human embryonic stem cell line. Skin.ker: foreskin keratinocyte. Skin.mel: foreskin melanocyte. Brain. hip: brain hippocampus middle. Brain. ctx: Brain dorsolateral prefrontal cortex. Panc.islet: pancreatic islets. GM12878: lymphoblastoid. HMEC: mammary epithelial. HSMM: skeletal muscle myoblasts. HUVEC: umbilical vein endothelial. K562: leukaemia. NHEK: epidermal keratinocyte. NHLF: lung fibroblast.Figure 17. Clustering of metastases with the paired primary tumor. Unsupervised hierarchical cluster analysis of our 156-CpG signature showing that metastases did not form a distinct subcluster and most clustered next to their paired primary tumor, as highlighted by the second color-coded row. This is the same unsupervised hierarchical cluster analysis shown in Figure 6, based on the beta-values of the 156 CpG-signature; here each pair of primary tumor and metastasis is shown with a different color.

[0065] NT: normal tissue. DTC: differentiated thyroid cancer. non-mDTC: non-metastatic primary DTC. mDTC: metastatic primary DTC. LNM: lymph node metastases. DM: distant metastases. FTC: follicular thyroid cancer. PTC: papillary thyroid cancer.

[0066] Figure 18. Validation of the 156 CpG-signature in validation series I. Unsupervised hierarchical cluster analysis of our 156-CpG signature in the discovery series, in validation series I, and in three thyroid cancer invasive cell lines (TPC1, BCPAP, and HTh83). All but one mDTC in validation series I clustered with the mDTC of the discovery series. The three cell lines formed a subcluster within the mDTCs. Sample characteristics, including risk category, time of diagnosis of DM, and histology are represented in color-coded rows below the upper dendogram.

[0067] NT: normal tissue. DTC: differentiated thyroid cancer. non-mDTC: non-metastatic primary DTC. mDTC: metastatic primary DTC. MET: metastases (both lymph node and distant metastases). FTC: follicular thyroid cancer. PTC: papillary thyroid cancer.

[0068] Figure 19. Hypermethylation and downregulation of ELTD1 gene in metastatic primary thyroid tumors compared to non-metastatic primary thyroid tumors. (A) Beta-values measured by EPIC array of the CpGs covering the ELTD1 -associated 5’ region, including a CpG island (CGI; green square), for normal tissue (NT) and non-metastatic and metastatic primary papillary thyroid cancer (non-PTC and mPTC, respectively) and follicular thyroid cancer (non-FTC and mFTC, respectively). The two CpGs included in our 156-signature (cg04360793 and cg15084543) are indicated with a dotted blue border. (B) Lollipop diagram showing the methylation level of the CpG sites(circles) within the 5’ region of ELTD1 assessed by bisulfite sequencing in a subset of samples from the discovery series. Cg04360793 and cg15084543 are indicated with a dotted blue border. The plot was obtained using the Methylation plotter tool (Mallona I, Diez-Villanueva A, Peinado MA. Methylation plotter: a web tool for dynamic visualization of DNA methylation data. Source Code Biol Med. 2014;9:11). (C) Boxplot confirming the hypermethylation of cg04360793 and cg15084543 within the ELTD1 promoter in TCGA high-risk vs low-risk PTC. (D) Boxplot showing the downregulation of ELTD1 gene in TCGA high-risk vs low-risk PTC.

[0069] Figure 20. Correlation between EPIC array and bisulfite pyrosequencing (BPS) results of the methylation levels of the four CpGs. The comparison of the DNA methylation levels in 64 samples from the discovery series showed a strong correlation for all four CpGs.

[0070] NT: normal tissue. DTC: differentiated thyroid cancer. non-mDTC: non-metastatic primary DTC. mDTC: metastatic primary DTC. LNM: lymph node metastases. DM: distant metastases.

[0071] Figure 21. Validation of the prognostic value of the global DNA methylation of the 156-CpG signature measured by EPIC array. (A) Global DNA methylation levels, measured as the mean of the beta-values of the 156-CpG signature, in the samples used in the validation series (n=91). The number of samples in each group is indicated in brackets on the X-axis. The Kruskal-Wallis test and Dunn's post hoc test were used to determine statistical significance. *P < 0.05; **P < 0.01; ***P < 0.001. NT, normal tissue, non-metastatic primary differentiated thyroid cancer (non-mDTC), metastatic primary DTC (mDTC), lymph node metastases (LNM) and distant metastases (DM). (B) Receiver operating characteristic (ROC) curve analysis of the average methylation level of the 156-CpG signature. AUC, area under the curve.

[0072] Figure 22. Validation of the prognostic value of the 4 CpG-based predictive model measured by EPIC array. (A) Dot plot displaying the predicted risk of developing distantmetastases (DM) in the samples used in the validation series (n=91). The prediction of the risk of DM and the sample type (normal tissue (NT), non-metastatic primary differentiated thyroid cancer (non-mDTC), metastatic primary DTC (mDTC), lymph node metastases (LNM) and DM) of each sample are indicated. (B) Receiver operating characteristic (ROC) curves showing the individual discriminatory capacity of each of the 4 CpG biomarkers. (C) Receiver operating characteristic (ROC) curves showing the discriminatory capacity of the predictive model combining the 4 CpG biomarkers. AUC, area under the curve.

[0073] Figure 23. Validation of the prognostic value of the 4 CpG-based predictive model measured by bisulfite pyrosequencing. (A) Dot plot displaying the predicted risk of developing distant metastases (DM) in the samples used in the validation series (n=82). The prediction of the risk of DM and the sample type (normal tissue (NT), non-metastatic primary differentiated thyroid cancer (non-mDTC), metastatic primary DTC (mDTC), lymph node metastases (LNM) and DM) of each sample are indicated. (B) Receiver operating characteristic (ROC) curves showing the individual discriminatory capacity of each of the 4 CpG biomarkers. (C) Receiver operating characteristic (ROC) curves showing the discriminatory capacity of the predictive model combining the 4 CpG biomarkers. AUC, area under the curve.

[0074] Figure 24. Prognostic value of the 5 CpG-based predictive model measured by EPIC array. (A) Mean beta-values of each of the five CpG biomarkers in normal tissue (NT), non-metastatic (non-m) and metastatic (m) primary differentiated thyroid cancer (DTC), lymph node metastases (LMN) and distant metastases (DM). The Kruskal-Wallis test and Dunn's post hoc test were used to determine statistical significance. *P < 0.05; **P < 0.01; ***P < 0.001. (B) Unsupervised hierarchical cluster analysis of the five CpG biomarkers in our discovery series. (C) Dot plot displaying the predicted risk of developing distant metastases (DM) in the samples used in the discovery series (n=65primary tumors). The prediction of the risk of DM based and the sample type (NT, non-mDTC, mDTC, LNM and DM) of each sample are indicated.

[0075] Figure 25. Validation of the prognostic value of the 5 CpG-based predictive model measured by EPIC array. (A) Dot plot displaying the predicted risk of developing distant metastases (DM) in the samples used in the validation series (n=91). The prediction of the risk of DM and the sample type (normal tissue (NT), non-metastatic primary differentiated thyroid cancer (non-mDTC), metastatic primary DTC (mDTC), lymph node metastases (LNM) and DM) of each sample are indicated. (B) Receiver operating characteristic (ROC) curves showing the individual discriminatory capacity of each of the 5 CpG biomarkers. (C) Receiver operating characteristic (ROC) curves showing the discriminatory capacity of the predictive model combining the 5 CpG biomarkers. AUC, area under the curve.

[0076] Description of embodiments

[0077] In the present invention we have profiled DNA methylation during DTC metastatic progression and found a progressively increasing number of alterations - mostly hypomethylation events - from normal tissue through DTC to DM, both in PTC and FTC, although DNA methylation dynamics differed between the two histological subtypes. In addition, based on our findings, we have generated a novel 4-CpG predictive model that can easily be assessed in routine clinical practice to accurately identify patients at risk of developing DM.

[0078] Our previous study on Alu repeats pointed out the association of global Alu hypomethylation with progression and dedifferentiation (Klein Hesselink EN, et al. Increased Global DNA Hypomethylation in Distant Metastatic and Dedifferentiated Thyroid Cancer. J Clin Endocrinol Metab. 2018;103(2):397-406). Accordingly, here we confirm the relevance of global DNA hypomethylation in the progression of DTC, and in addition show that hypomethylation affects not only Alu repeats but also uniquesequences within non-CGI regions. In fact, some studies using genome-wide sequencing in different cancer types have demonstrated that DNA hypomethylation occurs in large DNA domains containing hyper-variably expressed genes (Berman BP, et al. Regions of focal DNA hypermethylation and long-range hypomethylation in colorectal cancer coincide with nuclear lamina-associated domains. Nat Genet.

[0079] 2011;44(1):40-46; Hansen KD, et al. Increased methylation variation in epigenetic domains across cancer types. Nat Genet. 2011;26;43(8):768-75). It is well known that global DNA hypomethylation can cause chromosomal instability, derepression of imprinted genes and retrotransposons, and aberrant gene expression (Ehrlich M. DNA hypomethylation in cancer cells. Epigenomics. 2009;1(2):239-259); however, the functional implications in thyroid cancer progression are unknown.

[0080] In contrast to global hypomethylation of Alu repeats that appears to be a late event in the molecular history of thyroid cancer progression (Klein Hesselink EN, et al. Increased Global DNA Hypomethylation in Distant Metastatic and Dedifferentiated Thyroid Cancer. J Clin Endocrinol Metab. 2018;103(2):397-406), we have found a global hypomethylation of non-CGI regions also in low-risk PTC, but not in low-risk FTC. In this sense, different DNA methylome studies have reported more hypomethylations than hypermethylations in PTC, but not in FTC (Zafon C, et al. DNA methylation in thyroid cancer. Endocr Relat Cancer. 2019;26(7): R415-R439). However, our analysis has revealed for the first time that this difference occurs only between low-risk PTC and low-risk FTC, whereas both mPTC and mFTC harbor more hypomethylations than hypermethylations, pointing out that the dynamics of DNA methylation are different for PTC and FTC. Moreover, while the overlap of DNA methylation alterations between low-risk PTC and low-risk FTC is low, it is high between mPTC and mFTC or between PTC-DM and FTC-DM, which suggests that the epigenetic alterations underlying PTC and FTC initiation are histology type-specific, whereas those underlying metastatic progression are largely shared. Since most PTC are classified as BRAF-like tumorscharacterized by high MEK-ERK signaling and FTC are classified as RAS-like tumors characterized by low MEK-ERK signaling and concurrent activation of PI3K (The Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell. 2014;159(3):676-690), these signaling pathways may regulate, at least in part, these epigenetic differences in low-risk tumors.

[0081] We have found that hypermethylated CpGs are associated with CGIs and CGI shores, reminiscent of the classical CGI methylator phenotype (CIMP) (Issa JP. CpG-island methylation in aging and cancer. Curr Top Microbiol Immunol. 2000;249:101-18) and many studies in different cancer types (Berman BP, et al. Regions of focal DNA hypermethylation and long-range hypomethylation in colorectal cancer coincide with nuclear lamina-associated domains. Nat Genet. 2011;44(1):40-46.; Hansen KD, et al. Increased methylation variation in epigenetic domains across cancer types. Nat Genet.

[0082] 2011;26;43(8):768-75; Kulis M, et al. Epigenomic analysis detects widespread genebody DNA hypomethylation in chronic lymphocytic leukemia. Nat Genet.

[0083] 2012;44(11):1236-42). These CpGs are also enriched in bivalent promoters and enhancers and in Polycomb-repressed regions, as occurs in other cancers (Easwaran H, et al. A DNA hypermethylation module for the stem / progenitor cell signature of cancer. Genome Res. 2012;22(5):837-849; Rodriguez J, et al. Bivalent domains enforce transcriptional memory of DNA methylated genes in cancer cells. Proc Natl Acad Sci U S A. 2008; 105(50): 19809-19814), pointing out their repressive effect on gene expression. In low-risk PTC, however, hypermethylated CpGs are mostly located in nonCGI regions not associated with bivalent regulatory elements nor Polycomb, indicating a different role of DNA methylation in the initiation of PTC compared to FTC. On the other hand, most hypomethylations are located in intergenic regions within non-CGI regions and enriched in enhancers, this enrichment being more remarkable in PTC. Altogether, these results show the potential role of DNA methylation alterations in the regulation of gene expression and thus in shaping the tumor phenotype. Yet, further analysis shouldbe warranted to distinguish passenger DNA methylation alterations from those driving tumor progression (Heery R, Schaefer MH. DNA methylation variation along the cancer epigenome and the identification of novel epigenetic driver events. Nucleic Acids Res.

[0084] 2021; 49(22): 12692-12705).

[0085] Cancer cells need to acquire certain properties and functions to metastasize such as migration capacity and invasiveness, cellular plasticity, or colonization of distant sites, among others (Turajlic S, Swanton C. Metastasis as an evolutionary process. Science.

[0086] 2016;352(6282):169-175). Accordingly, functional enrichment analyses have shown that DNA methylation changes along PTC and FTC progression are associated with genes related to cell adhesion, whose disruption is key to allow metastatic dissemination. Interestingly, both low-risk PTC and FTC display a few differentially methylated cadherin and protocadherin genes, while the proportion of epigenetically altered cadherin and protocadherin genes is greatly increased in mPTC and mFTC and even more in DM. In this regard, many studies have shown the epigenetic deregulation of these genes in different tumor types, resulting hyper- or hypomethylated and conferring pro- or anti-metastatic capacities depending on the specific cadherin / protocadherin gene (Vega-Benedetti AF, et al. Clustered protocadherins methylation alterations in cancer. Clin Epigenetics. 2019; 11 (1): 100; van Roy F. Beyond E-cadherin: roles of other cadherin superfamily members in cancer. Nat Rev Cancer. 2014;14(2):121-134). Another result of interest is the enrichment of differentially methylated genes in low-risk PTC in the extracellular region related to secretion. Therefore, a future specific analysis of some of the differentially methylated genes identified along PTC and FTC metastatic progression may contribute to better understanding of the molecular pathways underlying cell dissemination, helping to identify potential new therapeutic targets and leading to improve treatments.

[0087] Unlike most cancer types, thyroid carcinomas tend to remain subclinical with no or very slow progression, as evidenced by the fact that DTC, especially PTC, are commonlyfound at autopsy in people who have died without knowing they had a thyroid cancer (Furuya-Kanamori L, et al. Prevalence of Differentiated Thyroid Cancer in Autopsy Studies Over Six Decades: A Meta-Analysis. J Clin Oncol. 2016;34(30):3672-3679). In contrast, a small subset of DTCs develops DM, often many years after initial diagnosis, but the natural history of metastatic progression is unknown. Our results show a progressive gain of DNA methylation alterations from low-risk DTC to mDTC and DM and, more importantly, reveal that a high fraction of DNA methylation alterations in low-risk PTC or low-risk FTC are present in mPTC and PTC-DM or mFTC and FTC-DM, respectively, suggesting that these changes may be important not only for tumor initiation but also for tumor maintenance, and pointing out a linear progression model. At the same time, we have found common DNA methylation alterations in high-risk primary tumors and DM that are not present in low-risk primary tumors, such as the 156-CpG signature, which suggests that it may be involved in the acquisition of the metastatic phenotype and may also be used as predictors of DM. If DM were seeded by a minor subclone in the paired primary tumor, as posited by the classical progression model (Turajlic S, Swanton C. Metastasis as an evolutionary process. Science. 2016;352(6282):169-175), these alterations would be unlike to be detected by bulk DNA methylation analysis in the primary tumor; thus, the high epigenetic similarity we have found between primary tumors and paired DM suggests that DM derive from a major subclone in the primary tumor (Figure 8). In the same line, and according to previous reports (Song E, et al. Genetic profile of advanced thyroid cancers in relation to distant metastasis. Endocr Relat Cancer. 2020;27(5):285-293; Sohn SY, et al. Highly Concordant Key Genetic Alterations in Primary Tumors and Matched Distant Metastases in Differentiated Thyroid Cancer. Thyroid. 2016;26(5):672-68), we have found a high concordance between BRAF, RAS and TERTp mutations of primary tumors and paired DM. In this regard, prolonged dormancy of DM followed by progression many years after the initial diagnosis and treatment, which is common in DTC (Rajan N, et al. Progression and dormancy in metastatic thyroid cancer: concepts and clinical implications.Endocrine. 2020;70(1):24-35), may explain metachronous DM harboring genetic and epigenetic alterations similar to primary tumors.

[0088] With clinical application in mind, we generated, from the 156-CpG signature, a predictive model based on 4 CpGs that can predict the risk of developing DM with a 90% accuracy and established simple and quantitative assays based on BPS to measure their DNA methylation. In this sense, in connection to the 156-CpG signature and the predictive model based on 4 CpGs, apart from the data exemplified herein, please note the following. Using EPIC array, we validated that the average DNA methylation of the 156 CpGs was lower in mDTC and metastases compared to non-mDTC and normal tissues in a validation series of 91 samples (20 normal tissues, 17 non-mDTC, 29 mDTC and 25 metastases). This global DNA methylation value was able to discriminate between samples with low- and high-risk of DM with an area under the curve (AUC) of 0.8, sensitivity of 78.6% and specificity of 82.4% (Figure 21A-B). In addition, Using EPIC array, we validated the predictive model based on the 4 CpG biomarkers in an independent series of 91 samples (20 normal tissues, 17 non-mDTC, 29 mDTC and 25 metastases). Figure 22A shows that all normal tissues and most non-mDTC have a low-risk of DM, as expected, while most mDTC and metastases have a high-risk of DM. The analysis of each individual CpG biomarker already achieved a good ability to discriminate samples with low- and high-risk of DM; however, when the 4-CpG biomarkers were combined in the predictive model, the discriminatory performance was further improved, yielding higher AUC, specificity and sensitivity, and demonstrating an enhanced capacity to predict DM (Figure 22B-C). Furthermore, using bisulfite pyrosequencing, an easy-to-implement technique to determine DNA methylation, we validated the 4 CpG biomarkers in an independent series including 82 samples (19 normal tissues, 33 non-mDTC, 19 mDTC and 11 metastases). Figure 23A shows that all normal tissues and most non-mDTC have a low-risk of developing DM, as expected, while most mDTC and metastases have a high-risk of developing DM. The analysis ofeach individual CpG biomarker already achieved a good ability to discriminate between samples with low- and high-risk of DM; however, when the 4-CpG biomarkers were combined in the predictive model, the discriminatory performance was further improved, yielding higher AUG, specificity and sensitivity, and demonstrating an enhanced capacity to predict DM (Figure 23B-C).

[0089] Moreover, in addition to the 4-CpG signature, we also identified five CpG biomarkers that were differentially methylated in normal tissue and non-mDTC compared to mDTC, LNM and DM: cg25901805 (located in PLAGL1 gene, chr6: 144,286, 162-144,286, 162), cg10780897 (located in RP1 gene, chr8:55, 528, 099-55, 528, 099), cg07530798 (located in chr1:119, 535, 589-119, 535, 589), cg04690464 (located in LINC02691 gene, chr14:104, 711, 952-104, 711, 952), and cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259), the last of which is also included in the 4-CpG signature (Figure 24A). This group of five CpG biomarkers was still able to discriminate between non-mDTC and mDTC (Figure 24B). Using these five CpG biomarkers, we generated a predictive model to determine the risk of a primary DTC developing DM. When we applied this predictive model to the 65 DTC samples included in the discovery series, using the beta values from the EPIC arrays, it correctly classified 62 of them (Figure 24C). Furthermore, using EPIC array, we validated the 5 CpG biomarkers in an independent series including 91 samples (20 normal tissues, 17 non-mDTC, 29 mDTC and 25 metastases). Figure 25A shows that all normal tissues and most non-mDTC have a low-risk of developing DM, as expected, while most mDTC and metastases have a high-risk of developing DM. The analysis of each individual CpG biomarker already achieved a good ability to discriminate between samples with low-and high-risk of DM; however, when the 5-CpG biomarkers were combined in the predictive model, the discriminatory performance was further improved, yielding higher AUC, specificity and sensitivity, and demonstrating an enhanced capacity to predict DM (Figure 25B-C).1

[0090] Thus, the predictive model is an easy-to-implement prognostic tool that will allow physicians to identify those patients at higher risk of developing DM in early stages of the disease to ensure a closer follow-up to find DM and an early start of the most appropriate treatment, including multikinase inhibitors (Cabanillas ME, et al. Targeted Therapy for Advanced Thyroid Cancer: Kinase Inhibitors and Beyond. Endocr Rev.

[0091] 2019;40(6):1573-1604), to prevent progression and improve survival. On the other hand, it will also help to reduce overtreatment (Ullmann TM, et al. Current Controversies in Low-Risk Differentiated Thyroid Cancer: Reducing Overtreatment in an Era of Overdiagnosis. J Clin Endocrinol Metab. 2023; 108(2):271-280) of patients with a low risk of developing DM and its harmful side effects, thus improving their quality of life. It would be even more useful to apply the predictive model to preoperative samples obtained by fine-needle aspiration (FNA), which would enhance the risk stratification preoperatively allowing a better planning of the surgery extent. Although we have not assessed the performance of the epigenetic assays in FNA samples, the bisulfite pyrosequencing technique has been demonstrated to be compatible with this routine clinical sample type (Bisarro dos Reis M, et al. Prognostic Classifier Based on Genome-Wide DNA Methylation Profiling in Well-Differentiated Thyroid Tumors. J Clin Endocrinol Metab. 2017;102(11):4089-4099; Rodriguez-Rodero S, et al. Classification of follicular-patterned thyroid lesions using a minimal set of epigenetic biomarkers. Eur J Endocrinol.

[0092] 2022;187(3):335-347).

[0093] To our knowledge, this is the first comprehensive epigenetic analysis in the largest series of mDTCs and paired metastases reported to date. In conclusion, our comprehensive DNA methylation analysis unveils new insights into the dynamics and role of DNA methylation in thyroid cancer development and metastatic progression as well as into the natural history of DM and provides valuable information from a clinical perspective. In particular, we have generated an accurate and easy-to-implement prognostic tool topredict the risk of developing DM, which may improve the current clinical DTC patient risk classification and help move towards personalized medicine in thyroid cancer.

[0094] Therefore, a first aspect of the invention refers to an in vitro method of determining the methylation level of a set of biomarkers in a subject having thyroid cancer, the method comprising:

[0095] a. from an isolated or extracted genomic DNA obtained from a biological sample from the subject, determining the methylation level of at least 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150 or 156 biomarkers selected from Table 1A, Table 1B and / or Table 1C from the extracted genomic DNA,

[0096] where said determination of the methylation level is preferably used to determine the risk of developing metastasis, preferably distant metastasis (DM), in said subject.

[0097] The inventors have surprisingly found that the CpG biomarkers listed in Table 1A, Table 1B and Table 1C were differentially methylated between non-mDTC and mDTC. The listed CpGs can be associated with thyroid cancer, can be involved in the acquisition of the metastatic phenotype and can be used as predictors of developing distant metastasis.

[0098] The total set of biomarkers listed in Table 1A, Table 1B and Table 1C can be considered as the 156 CpG signature. The 156 CpG biomarkers differentially methylated between the mDTC and the non-mDTC (minimum of 20% of methylation difference between both groups) were determined through statistical analysis.

[0099] Table 1A. List of 4 CpG biomarkers associated with thyroid cancer and distant metastasis risk. The correspondence between the nomenclature used by Illumina array (EPIC array ID) and the loci in the human genome (GRCh37 / hg19) is provided.

[0100] CpG 1 CpG 2 CpG 3 CpG 4

[0101]

[0102] EPIC array ID cg01378044 cg18305394 cg11925561 cg14672100 Gene name TBX18 PCK1

[0103] Gene-related

[0104] Gene body 5’ region

[0105] region

[0106] Chromosome 6 20 5 5 Location (Feb. 85462137 56134788 3304588 6344259 2009 or or or or (GRCh37 / hg1 chr6:85,462,13 chr20:56, 134,78 chr5:3,304,58 chr5:6,344,25 9)) 7-85,462,137 8-56,134,788 8-3,304,588 9-6,344,259

[0107]

[0108] Table 1B. List of 5 CpG biomarkers associated with thyroid cancer and distant metastasis risk. The correspondence between the nomenclature used by the Illumina array (EPIC array ID) and the loci in the human genome (GRCh37 / hg19) is provided.

[0109] CpG 4 CpG 5 CpG 6 CpG 7 CpG 8 EPIC array eg 146721 cg25901805 eg 1078089 cg07530798 cg04690464 ID 00 7

[0110] Gene LINC02691 PLAGL1 RP1

[0111] name

[0112] Gene- related Gene body Gene body Gene body region

[0113] Chromoso

[0114] 5 6 8 1 14

[0115] me

[0116] 6344259 144286162 55528099 119535589 104711952 Location

[0117] or or or or or

[0118] (Feb. 2009

[0119] chr5:6,344 chr6:144,28 chr8:55,52 chr1:119,53 chr14:104,71 (GRCh37 / h

[0120] ,259- 6,162- 8,099- 5,589- 1,952- g19))

[0121]

[0122] 6,344,259 144,286,162 55,528,099 119,535,589 104,711,952

[0123] Table 1C. List of 148 CpG biomarkers associated with thyroid cancer and distant metastasis risk. The correspondence between the nomenclature used by the Illumina array (EPIC array ID) and the loci in the human genome (GRCh37 / hg19) is provided.

[0124] EPIC array ID Location (Feb. 2009 (GRCh37 / hg19))

[0125] cg26785499 chr1:2820433-2820433

[0126] cg01431482 chr1:2989085-2989085

[0127] cg26139866 chr1:4239037-4239037

[0128]

[0129] cg14035553 chr1:21029443-21029443cg09551984 chr1: 30229976-30229976 cg24364593 chr1: 34533065-34533065 cg06437931 chr1:38180356-38180356 cg10776919 chr1: 50889442-50889442 cg18279094 chr1:63790044-63790044 cg04360793 chr1:79472361 -79472361 cg15084543 chr1: 79472408-79472408 cg21618157 chr1:92202671-92202671 cg24398933 chr1:151492569-151492569 cg19387132 chr1:158254116-158254116 cg27534912 chr1:158258802-158258802 cg14275423 chr1: 158656484- 158656484 cg05304729 chr1: 158800024-158800024 cg25743622 chr1:223403446-223403446 cg26404722 chr2:9353201 -9353201 cg11448675 chr2:9375570-9375570 cg15366555 chr2:74601926-74601926 cg11412468 ch r2: 84433405-84433405 cg14231073 chr2:121347537-121347537 cg22540067 chr2:125305001-125305001 cg21445398 chr2: 125425654- 125425654 cg05053811 chr2:126401437-126401437 cg12906989 chr2: 137925983-137925983 cg12341429 chr2:161887154-161887154 cg18621766 chr2: 194689205-194689205 cg04314280 chr2:218687427-218687427 cg15446845 chr3: 5470909-5470909 cg22438338 chr3:14634313-14634313 cg18780425 chr3:69250024-69250024 cg11963062 chr3: 100633394- 100633394 cg16768018 chr3: 147108843-147108843 cg04871469 chr4:2257119-2257119 cg25175009 chr4:63265254-63265254 cg18549952 chr4: 129762662- 129762662 cg04783185 chr4: 190521538- 190521538 cg01653189 chr5:2569982-2569982 cg07826516 chr5:3106230-3106230 cg11328079 chr5: 3504656-3504656 cg26976756 chr5:3510004-3510004 cg18916149 chr5: 3532269-3532269 cg20009810 chr5: 3537215-3537215 cg09696091 chr5: 3564465-3564465 cg13453847 chr5: 5034274-5034274 cg03290223 chr5:5063319-5063319 cg06587623 chr5:5079159-5079159 cg20997704 chr5:13143814-13143814 cg00086247 chr5:29568333-29568333 cg17738861 chr5:32233621-32233621 cg04206637 chr5:63202277-63202277 cg07497602 chr5: 105881003-105881003 cg01664151 chr5: 132945292- 132945292 cg20876760 chr6: 9399379-9399379

[0130]

[0131] cg00237606 chr6:29455256-29455256cg22424903 chr6:35454434-35454434 cg24749265 chr6:41216090-41216090 cg25075684 chr7:2933297-2933297 cg08241694 chr7:50633896-50633896 cg16455187 chr7:50635841-50635841 cg06915574 chr7:76639592-76639592 cg01611504 chr7:108525503-108525503 cg10589286 chr7: 154916724- 154916724 cg15983026 chr7: 155325552-155325552 cg03974619 chr7: 156851553-156851553 cg06437651 chr7: 157466201 - 157466201 cg15819924 chr7: 157886456- 157886456 cg06593940 chr7: 157950633- 157950633 cg05647867 chr7: 158244818-158244818 cg07178994 chr8:816998-816998 cg07340894 chr8: 8238774-8238774 cg14793063 chr8:30557532-30557532 cg15325961 chr8:49183143-49183143 cg17406248 chr9:97807693-97807693 cg13811240 chr9:97851374-97851374 cg04027576 chr9:117996948-117996948 eg 19537600 chr9: 119836588-119836588 cg16681908 chr9: 140782525- 140782525 cg03478630 chr10: 10436568-10436568 cg11601663 chr10: 15666820-15666820 cg09918581 chr10:68812190-68812190 cg15438734 chr10:77811027-77811027 cg20069378 chr10:98117247-98117247 cg27567836 chr10:124091022-124091022 cg11174855 chr10:134598352-134598352 cg08264839 chr11:3191337-3191337 cg05133853 chr11:5265008-5265008 cg10092198 chr11:6129548-6129548 cg16573636 chr11:19917499-19917499 cg18157275 chr11: 56345140-56345140 cg04562441 chr11: 56380921-56380921 cg22151881 chr11:68082621-68082621 cg05355757 chr12:16762831-16762831 cg14096615 chr12: 16762843- 16762843 cg19104050 chr12:108701064-108701064 cg26998150 chr12:130734158-130734158 cg15858166 ch r12: 130953551 - 130953551 cg06626892 chr13:78995230-78995230 cg24690071 chr13:100635352-100635352 cg23117999 chr14:96508302-96508302 cg23450358 chr14:97581404-97581404 cg20925031 chr15:29988135-29988135 cg09262300 chr15:86519955-86519955 cg09783969 chr15:100886002-100886002 cg01252526 chr16:711033-711033 cg03287339 chr16:711135-711135 cg27039467 chr16:1730561-1730561

[0132]

[0133] cg13985307 chr16:25435877-25435877cg19530424 chr16:51795292-51795292

[0134] cg04770195 chr16:53114571-53114571

[0135] cg04980805 chr16:82377691-82377691

[0136] cg16122778 chr16:86671056-86671056

[0137] cg16596367 chr16:89070757-89070757

[0138] cg21851534 chr17:3907994-3907994

[0139] cg14601733 chr17:6280370-6280370

[0140] cg03916694 chr17: 18150707-18150707

[0141] cg06603761 chr17:32691699-32691699

[0142] cg23368579 chr18:519265-519265

[0143] cg19291618 chr18:4599319-4599319

[0144] cg15832062 chr19:3282509-3282509

[0145] cg19058865 chr19:4327350-4327350

[0146] cg19839825 chr19:7616011-7616011

[0147] cg20185461 chr19:18899082-18899082

[0148] cg19386774 chr19:22467847-22467847

[0149] cg08584061 chr19:22718926-22718926

[0150] cg17699869 chr19:29421611-29421611

[0151] cg20804418 chr19:32614622-32614622

[0152] cg22836826 chr19:42194016-42194016

[0153] cg11083325 chr19:49244668-49244668

[0154] cg16509509 chr19:51151383-51151383

[0155] cg26324258 chr20:40612816-40612816

[0156] cg21939911 chr20:44644971-44644971

[0157] cg25535351 chr20:50956552-50956552

[0158] cg02762107 chr20:51335826-51335826

[0159] cg16493240 chr20:58156194-58156194

[0160] cg16326851 chr20:60094855-60094855

[0161] cg05198244 chr21:26732918-26732918

[0162] cg02442623 chr21:30673051-30673051

[0163] cg14861920 chr22:35019460-35019460

[0164] cg01586506 chr22:38379506-38379506

[0165] cg13210663 chr22:47790645-47790645

[0166] cg04768709 chr22:48080590-48080590

[0167] cg16689883 chr22:48441270-48441270

[0168] cg27334056 chr22:48479650-48479650

[0169] cg17659636 chr22:48581422-48581422

[0170]

[0171] cg19491717 chr22:48843289-48843289

[0172] From hereinafter, we shall refer to each methylation site, also named CpG biomarker, in the present application with the nomenclature used by the Illumina array (EPIC array ID). The correlation with its exact locus is shown in Table 1A, Table 1B and Table 1C. For instance, the methylation site chr6:85,462, 137-85,462, 137 corresponds to cg01378044 according to EPIC array ID nomenclature.A “reference genome” or “reference assembly” is a digital nucleic acid sequence database, assembled as a representation of the set of genes in one idealized individual organism of a species. It is a set of the assembled sequences of DNA from a number of individual donors, therefore they do not represent the set of genes of any single individual organism, but it provides a haploid mosaic of different DNA sequences from each donor. GRCh37 / hg19 refers to the genome assembly Genome Reference Consortium Human Build 37 (GRCh37) from Homo sapiens (human) and was submitted by the Genome Reference Consortium on 2009 / 02 / 27. A synonym is hg19. The assembly type is “haploid-with-alt-loci”, the assembly level is “chromosome” and the genome representation is full. GRCh37 corresponds to the accession PRJNA31257 of the BioProject. GRCh37 / hg19 corresponds to GenBank assembly accession: GCA_000001405.1, also to RefSeq assembly accession: GCF_000001405.13. The RefSeq assembly and GenBank assembly are identical (National Centre Biotechnology Information. “GRCh37”. Accessed January 25, 2024. https: / / www.ncbi.nlm.nih.gov / assembly / GCF_000001405.13 / ?shouldredirect=false# / def ). The gene position can be accessed through Genome Browser Gateway by selecting Human Assembly Feb. 2009 (GRCh37 / hg19) (University of California Santa Cruz (UCSC). “Genome Browser Gateway”. Accessed January 25, 2024 htps: / / genome.ucsc.edu / cgi-bin / hgGateway?redirect=manual&source=genome. ucsc.edu).

[0173] ROC analysis is used in clinical epidemiology to quantify how accurately medical diagnostic tests can discriminate between two patient states, such as “diseased” and “non-diseased”. The area under a receiver operating characteristic (ROC) curve, also AUC, is a single scalar value that measures the global performance of a binary classifier (Hanley JA, McNeil BJ. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology. 1982;143(1):29-36. doi:10.1148 / radiology.143.1.7063747). The AUC value is found within the range [0.5-1.0], where the minimum value indicates the performance of a random classifier and the maximum value would represent to a perfect classifier (e.g., with a classification error rate equivalent to zero). The AUC is a robust total measure to evaluate the performance of score classifiers since its calculation relies on the complete ROC curve and thus involves all possible classification thresholds. The AUC is normally calculated by addition of successive trapezoid areas below the ROC curve (Melo, F.. Area under the ROC Curve. In: Dubitzky, W., Wolkenhauer, O., Cho, KH., Yokota, H. (eds) Encyclopedia of Systems Biology. Springer, New York, NY. (2013). https: / / doi.org / 10.1007 / 978-1-4419-9863-7_209). The ROC curve is based on the notion of a "separator" scale, on which results for the diseased and non-diseased form a pair of overlapping distributions. The complete separation of the two underlying distributions implies a perfectly discriminating test while complete overlap implies no discrimination. The ROC curve shows the trade off between the true positive fraction (TPF) or sensitivity and false positive fraction (FPF) or 1 -specificity as one change the criterion for positivity. The ROC curve is based on the notion of a separator (or decision) variable. The frequencies of positive and negative results of the diagnostic test will vary if one changes the criterion or cut-off for positivity on the decision axis. The optimal cut-off value can be determined using ROC curve analysis. The Youden index is one of the methods that have been proposed to determine optimal cut off values. The Youden index uses the maximum of vertical distance of ROC curve from the point (x,y) on diagonal line (chance line). Actually, Youden index maximizes the difference between sensitivity and 1 -specificity: Youden Index = Sensitivity + Specificity -1. Hence, by maximizing Sensitivity + Specificity across various cut-off points, the optimal cut-off point is calculated (Hajian-Tilaki K. Receiver Operating Characteristic (ROC) Curve Analysis for Medical Diagnostic Test Evaluation. Caspian J Intern Med. 2013 Spring;4(2):627-35).

[0174] Preferably, the method of the first aspect of the invention comprises determining the methylation level of at least one of the biomarkers selected from Table 1A. In particular,the method comprises determining the methylation level of cg01378044 (TBX18 gene). In another embodiment, the method comprises determining the methylation level of cg18305394 (5’ region of PCK1 gene). In another embodiment, the method comprises determining the methylation level of intergenic region cg11925561 in chromosome 5. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100 in chromosome 5.

[0175] In another embodiment of the first aspect, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A. In particular, the method comprises determining the methylation level of cg01378044 (TBX18 gene) and cg18305394 (5’ region of PCK1 gene). In another embodiment, the method comprises determining the methylation level of cg01378044 (TBX18 gene) and intergenic region cg14672100 in chromosome 5. In a preferred embodiment, the method comprises determining the methylation level of cg01378044 (TBX18 gene) and intergenic region cg11925561 in chromosome 5. In another embodiment, the method comprises determining the methylation level of cg18305394 (5’ region of PCK1 gene) and intergenic region cg11925561 in chromosome 5. In another embodiment, the method comprises determining the methylation level of cg18305394 (5’ region of PCK1 gene) and intergenic region cg14672100 in chromosome 5. In another embodiment, the method comprises determining the methylation level of intergenic region cg11925561 in chromosome and intergenic region cg14672100 in chromosome 5.

[0176] In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A. In particular, the method comprises determining the methylation level of cg01378044 (TBX18 gene), cg18305394 (5’ region of PCK1 gene) and intergenic region cg11925561 in chromosome 5. In another embodiment, the method comprises determining the methylation level of cg01378044 (TBX18 gene), cg18305394 (5’ region of PCK1 gene) and intergenic region cg14672100 in chromosome 5. In another embodiment, the method comprises determining themethylation level of cg18305394 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 and intergenic region cg14672100 in chromosome 5. In another embodiment, the method comprises determining the methylation level of cg01378044 (TBX18 gene), intergenic region cg11925561 in chromosome 5 and intergenic region cg14672100 in chromosome 5.

[0177] In a more preferred embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1A. Preferably, the method comprises determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 (TBX18 gene), cg18305394 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 and intergenic region cg14672100 in chromosome 5.

[0178] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1C. In particular, the method comprises determining the methylation level of cg01378044 (TBX18 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg18305394 (5’ region of PCK1 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg11925561 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C.

[0179] In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1C. In particular, the method comprises determining the methylation level of cg01378044 (TBX18 gene) and cg18305394 (5’ region of PCK1gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg01378044 (TBX18 gene) and intergenic region cg14672100 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg01378044 (TBX18 gene) and intergenic region cg11925561 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg18305394 (5’ region of PCK1 gene) and intergenic region cg11925561 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg18305394 (5’ region of PCK1 gene) and intergenic region cg14672100 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg11925561 in chromosome 5 and intergenic region cg14672100 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C.

[0180] In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1C. In particular, the method comprises determining the methylation level of cg01378044 (TBX18 gene), cg18305394 (5’ region of PCK1 gene) and intergenic region cg11925561 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg01378044 (TBX18 gene), cg18305394 (5’ region of PCK1 gene) and intergenic region cg14672100 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg18305394 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 and intergenicregion cg14672100 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of cg01378044 (TBX18 gene), intergenic region cg11925561 in chromosome 5 and intergenic region cg14672100 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C.

[0181] In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1C. In particular, the method comprises determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 (TBX18 gene), cg18305394 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 and intergenic region cg14672100 in chromosome 5 in combination with at least one of the biomarkers selected from Table 1C.

[0182] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1B. In particular, the method comprises determining the methylation level of intergenic region cg14672100. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene). In another embodiment, the method comprises determining the methylation level of cg10780897 (RP1 gene). In another embodiment, the method comprises determining the methylation level of cg07530798. In another embodiment, the method comprises determining the methylation level of cg04690464 (UNC02691 gene). In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1 B. In particular, the method comprises determining the methylation level of intergenic region cg14672100 and cg25901805 (PLAGL1 gene). In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100 and cg10780897 (RP1 gene). In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100 and cg07530798. In another embodiment, the methodcomprises determining the methylation level of intergenic region cg14672100 and cg04690464 (LINC02691 gene). In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene) and cg10780897 (RP1 gene). In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene) and cg07530798. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene) and cg04690464 (LINC02691 gene). In another embodiment, the method comprises determining the methylation level of cg10780897 (RP1 gene) and cg07530798. In another embodiment, the method comprises determining the methylation level of cg10780897 (RP1 gene) and cg07530798. In another embodiment, the method comprises determining the methylation level of cg07530798 and cg04690464 (LINC02691 gene).

[0183] In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1B. In particular, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene) and cg10780897 (RP1 gene). In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene) and cg07530798. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene) and cg04690464 (LINC02691 gene). In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg10780897 (RP1 gene) and cg07530798. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg10780897 (RP1 gene) and cg04690464 (LINC02691 gene). In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene) and cg07530798. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene),cg10780897 (RP1 gene) and cg04690464 (UNC02691 gene). In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg07530798 and cg04690464 (UNC02691 gene). In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene), cg07530798 and cg04690464 (UNC02691 gene). In another embodiment, the method comprises determining the methylation level of cg10780897 (RP1 gene), cg07530798 and cg04690464 (LINC02691 gene).

[0184] In another embodiment, the method comprises determining the methylation level of at least four of the biomarkers selected from Table 1B. In particular, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene) and cg07530798. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene) and cg04690464 (LINC02691 gene). In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene), cg07530798 and cg04690464 (LINC02691 gene). In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene), cg07530798 and cg04690464 (UNC02691 gene). In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg10780897 (RP1 gene), cg07530798 and cg04690464 (LINC02691 gene).

[0185] In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1B. Preferably the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene), cg07530798 and cg04690464 (LINC02691 gene).

[0186] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least one ofthe biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1 A in combination with at least one of the biomarkers selected from Table 1 B.

[0187] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least two of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A in combination with at least two of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least two of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1 A in combination with at least two of the biomarkers selected from Table 1B.

[0188] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1 A in combination with at least three of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1 A in combination with at least three of the biomarkers selected from Table 1 B. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least three of the biomarkers selected from Table 1B. In another embodiment, the method comprisesdetermining the methylation level of at least all of the biomarkers selected from Table 1 A in combination with at least three of the biomarkers selected from Table 1 B.

[0189] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least four of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1 A in combination with at least four of the biomarkers selected from Table 1 B. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least four of the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1 A in combination with at least four of the biomarkers selected from Table 1B.

[0190] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least all the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A in combination with at least all the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least all the biomarkers selected from Table 1B. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1A in combination with at least all the biomarkers selected from Table 1B.

[0191] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1B in combination with at least one of the biomarkers selected from Table 1C. In particular, the method comprises determining the methylation level of intergenic region cg14672100 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprisesdetermining the methylation level of cg25901805 (PLAGL1 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg10780897 (RP1 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C.

[0192] In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1B in combination with at least one of the biomarkers selected from Table 1C. In particular, the method comprises determining the methylation level of intergenic region cg14672100 and cg25901805 (PLAGL1 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100 and cg10780897 (RP1 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100 and cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100 and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene) and cg10780897 (RP1 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene) and cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene) andcg04690464 (UNC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg10780897 (RP1 gene) and cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg10780897 (RP1 gene) and cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg07530798 and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C.

[0193] In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1B in combination with at least one of the biomarkers selected from Table 1C. In particular, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene) and cg10780897 (RP1 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene) and cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene) and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg10780897 (RP1 gene) and cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg10780897 (RP1 gene) and cg04690464 (UNC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene),cg10780897 (RP1 gene) and cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene) and cg04690464 (UNC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg07530798 and cg04690464 (UNC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene), cg07530798 and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg10780897 (RP1 gene), cg07530798 and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C.

[0194] In another embodiment, the method comprises determining the methylation level of at least four of the biomarkers selected from Table 1B in combination with at least one of the biomarkers selected from Table 1C. In particular, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene) and cg07530798 in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene) and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene), cg07530798 and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene), cg07530798 and cg04690464(UNC02691 gene) in combination with at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of intergenic region cg14672100, cg10780897 (RP1 gene), cg07530798 and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C.

[0195] In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1 B in combination with at least one of the biomarkers selected from Table 1C. Preferably the method comprises determining the methylation level of intergenic region cg14672100, cg25901805 (PLAGL1 gene), cg10780897 (RP1 gene), cg07530798 and cg04690464 (LINC02691 gene) in combination with at least one of the biomarkers selected from Table 1C.

[0196] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1 B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least one of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1 A in combination with at least one of the biomarkers selected from Table 1 B and at least one of the biomarkers selected from Table 1C.

[0197] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least two of the biomarkers selected from Table 1B and at least one of the biomarkers selected fromTable 1C. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A in combination with at least two of the biomarkers selected from Table 1 B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least two of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1 A in combination with at least two of the biomarkers selected from Table 1 B and at least one of the biomarkers selected from Table 1C.

[0198] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1 A in combination with at least three of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A in combination with at least three of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least three of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1A in combination with at least three of the biomarkers selected from Table 1 B and at least one of the biomarkers selected from Table 1C.

[0199] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least four of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylationlevel of at least two of the biomarkers selected from Table 1A in combination with at least four of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least four of the biomarkers selected from Table 1 B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1A in combination with at least four of the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C.

[0200] In another embodiment, the method comprises determining the methylation level of at least one of the biomarkers selected from Table 1A in combination with at least all the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least two of the biomarkers selected from Table 1A in combination with at least all the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least three of the biomarkers selected from Table 1A in combination with at least all the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of at least all of the biomarkers selected from Table 1 A in combination with at least all the biomarkers selected from Table 1B and at least one of the biomarkers selected from Table 1C.

[0201] In another embodiment, the method comprises determining the methylation level as in any of the previous embodiments in combination with at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, and / or 148 biomarkers selected from Table 1C.In another embodiment, the method comprises determining the methylation level of all biomarkers selected from Table 1A and all biomarkers selected from Table 1B in combination with at least one biomarker selected from Table 1C. In another embodiment, the method comprises determining the methylation level of all biomarkers selected from Table 1A and all biomarkers selected from Table 1B in combination with at least 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, and / or 148 biomarkers selected from Table 1C. In another embodiment, the method comprises determining the methylation level of all biomarkers selected from Table 1A, Table 1B and Table 1C.

[0202] Examples of methods for determining DNA methylation include, but are not limited to, bisulfite-based methods (e.g. Sanger bisulfite sequencing), methylation-sensitive restriction enzyme-based methods, microarrays-based methods, next-generation sequencing (NGS) methods, methylation-specific melting curve analysis, mass spectrometry, immunoprecipitation methods, pyrosequencing, digital polymerase chain reaction (digital PCR), quantitative polymerase chain reaction (qPCR) or any combination thereof.

[0203] As used herein, the terms “digital PCR”, “single molecule PCR” and “single molecule amplification” refer to PCR and other nucleic acid amplification methods that are configured to provide amplification product or signal from a single starting molecule. Typically, samples are divided, e.g., by serial dilution or by partition into small enough portions (e.g., in microchambers or in emulsions) such that each portion or dilution has, on average, no more than a single copy of the target nucleic acid. Methods of single molecule PCR are described, e.g., in U. S. Pat. No. 6,143,496, which relates to a method comprising dividing a sample into multiple chambers such that at least one chamber has at least one target, and amplifying the target to determine how many chambers had a target molecule; U. S. Pat. No. 6,391,559; which relates to an assembly for containing and portioning fluid; and U. S. Pat. No. 7,459,315, whichrelates to a method of dividing a sample into an assembly with sample chambers where the samples are partitioned by surface affinity to the chambers, then sealing the chambers with a curable “displacing fluid.” See also U. S. Pat. No. 6,440,706 and U. S. Pat. No. 6,753,147, and Vogelstein, et al., Proc. Natl. Acad. Sci. USA Vol. 96, pp. 9236-9241, August 1999. See also US 20080254474, describing a combination of digital PCR combined with methylation detection. Preferably, the methylation levels are determined by a quantitative assay based on digital PCR.

[0204] More preferably, the methylation levels are determined by a quantitative assay based on bisulfite pyrosequencing. DNA is treated with bisulfite treatment, which converts unmethylated cytosine to uracil (which is subsequently converted to thymine during PCR amplification) while methylated cytosine remains unchanged. Next, bisulfite-treated DNA is amplified by polymerase chain reaction (PCR) and the resulting amplicons are sequenced by pyrosequencing, a sequence-by-synthesis method that measures the release of pyrophosphate as a byproduct of nucleotide incorporation by DNA polymerase, using an enzyme cascade that ends with detection of luciferase chemiluminescent signals.

[0205] It is herein noted that DNA methylation is the attachment of a methyl group at the C5 -position of the nucleotide base cytosine. In some instances, methylation of cytosine occurs in the CpG dinucleotides motif. In other instances, cytosine methylation occurs in, for example CHG and CHH motifs, where H is adenine, cytosine or thymine. In some instances, a CpG island is present in the 5' region of about one half of all human genes. CpG islands are typically, but not always, between about 0.2 to about 1 kb in length. Cytosine methylation further comprises 5 -methyl cytosine (5-mCyt) and 5-hydroxymethylcytosine. CGI shores are defined as regions up to 2 kb from the CGI boundary.

[0206] The CpG (cytosine-phosphate-guanine) or CG motif (which is the nomenclature used herein) refers to regions of a DNA molecule where a cytosine nucleotide occurs next toa guanine nucleotide in the linear strand. In some instances, a cytosine in a CpG dinucleotide is methylated to form 5 -methyl cytosine.

[0207] In some instances, one or more DNA regions are hypermethylated. In such cases, hypermethylation refers to an increase in the methylation of a region relative to a reference region. In some cases, hypermethylation is observed in one or more cancer types, and is useful, for example, as a diagnostic marker and / or a prognostic marker. In some instances, one or more DNA regions are hypomethylated. In some cases, hypomethylation refers to a loss in the methylation of a region relative to a reference region. In some cases, hypomethylation is observed in one or more cancer types, and is useful, for example, as a diagnostic marker and / or a prognostic marker.

[0208] In some embodiments, disclosed herein are CpG or CG methylation markers for prognosis of a cancer in a subject.

[0209] In particular, a second aspect of the invention refers to an in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:

[0210] a. determining the methylation level of one or more biomarkers as defined in the first aspect of the invention or in any of its preferred embodiments; and

[0211] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by a classification method that identifies the subject as a subject at risk of developing metastasis based on the biomarkers identified in step (a), preferably by a predictive model which correlates the methylation level of one or more of the biomarkers identified in step (a) with representative methylation level profiles from samples obtained from subjects having thyroid cancer previously identified as suffering from metastasis, preferably of distant metastasis(DM), said predictive model having been generated by training a computer with at least a plurality of methylation level profiles from previously identified subjects having thyroid cancer and suffering from metastasis, preferably of distant metastasis (DM), by machine learning on said plurality of methylation level profiles so as to obtain representative methylation level profiles associated with metastasis, preferably of distant metastasis (DM).

[0212] In a preferred embodiment of the method of the second aspect, the method comprises:

[0213] a. from an extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 or chr6:85,462,137-85,462, 137 (TBX18 gene), cg18305394 or chr20:56, 134,788-56, 134,788 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6,344,259- 6,344,259; and

[0214] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by a predictive model which correlates the methylation level of at least the biomarkers identified in step (a) with representative methylation level profiles from samples obtained from subjects having thyroid cancer previously identified as suffering from metastasis, preferably of distant metastasis (DM), said predictive model having been generated by training a computer with at least a plurality of methylation level profiles from previously identified subjects having thyroid cancer and suffering from metastasis, preferably of distant metastasis (DM), by machine learning on said plurality of methylationlevel profiles so as to obtain representative methylation level profiles associated with metastasis, preferably of distant metastasis (DM).

[0215] In another preferred embodiment of the method of the second aspect, the method comprises:

[0216] a. from an extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the methylation level of at least all of the biomarkers selected from the group consisting of cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259), cg25901805 (located in PLAGL1 gene, chr6:144, 286, 162-144, 286, 162), cg10780897 (located in RP1 gene, chr8:55, 528, 099-55, 528, 099), cg07530798 (located in chr1:119,535,589- 119,535,589) and cg04690464 (located in LINC02691 gene, chr14:104,711,952-104,711,952); and

[0217] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by a predictive model which correlates the methylation level of at least the biomarkers identified in step (a) with representative methylation level profiles from samples obtained from subjects having thyroid cancer previously identified as suffering from metastasis, preferably of distant metastasis (DM), said predictive model having been generated by training a computer with at least a plurality of methylation level profiles from previously identified subjects having thyroid cancer and suffering from metastasis, preferably of distant metastasis (DM), by machine learning on said plurality of methylation level profiles so as to obtain representative methylation level profiles associated with metastasis, preferably of distant metastasis (DM).In other embodiments, the methylation pattern or level can be determined by, but is not limited to, bisulfite-based methods, methylation-sensitive restriction enzyme-based methods, microarrays-based methods, next-generation sequencing (NGS), methylation-specific melting curve analysis, mass spectrometry, immunoprecipitation methods, pyrosequencing, digital polymerase chain reaction (digital PCR), quantitative polymerase chain reaction (qPCR) or any combination thereof. In a preferred embodiment, the methylation level is determined by a quantitative assay based on bisulfite pyrosequencing or digital PCR to measure the DNA methylation.

[0218] In an embodiment, the second step of the second aspect of the invention, step b), is performed by a machine learning method selected from a regression method, a classification method or a combination thereof.

[0219] It will be appreciated that the term "machine learning" generally refers to algorithms that give a computer the ability to learn without being explicitly programmed, including algorithms that learn from and make predictions about data. Machine learning algorithms employed by the embodiments disclosed herein may include, but are not limited to, random forest (" RF"), least absolute shrinkage and selection operator (" LASSO") logistic regression, regularized logistic regression, XGBoost, decision tree learning, artificial neural networks (" ANN"), deep neural networks (" DNN"), support vector machines, rulebased machine learning, and / or others.

[0220] For clarity, algorithms such as linear regression or logistic regression can be used as part of a machine learning process. However, it will be understood that using linear regression or another algorithm as part of a machine learning process is distinct from performing a statistical analysis such as regression with a spreadsheet program. Whereas statistical modeling relies on finding relationships between variables (e.g., mathematical equations) to predict an outcome, a machine learning process may continually update model parameters and adjust a classifier as new data becomes available, without relying on explicit or rules-based programming.In a particular embodiment, the second step of the second aspect of the invention is performed by a classification method, which results in identifying the subject as a subject at risk of metastasis.

[0221] In one embodiment, step (b) is carried out by a classification method; preferably selected from logistic regression, random forest, gradient boosting (GB), adaptive boosting (AB), extreme Gradient Boosting (XGB) k-nearest neighbors (kNN), artificial neural network (ANN), support vector machine (SVM), and combinations thereof.

[0222] In an embodiment, the predictive model is generated by training the computer with a plurality of methylation profiles from previously identified samples from subjects having thyroid cancer and suffering from DM by machine learning on said plurality of methylation profiles so as to obtain representative multivariable data sets associated with this group of patients; wherein the training comprises the following steps: (i) training data, from a plurality of methylation profiles, is randomly stratified into:

[0223] - a training or calibration dataset (for example in a percentage of about 75%), and - a validation dataset (for example in a percentage of about 25%);

[0224] (ii) the predictive model is seeded on the calibration dataset (particularly is developed by applying a machine learning method selected from a regression method, a classification method or a combination thereof on the calibration dataset);

[0225] (iii) the predictive model is optimized by an internal cross validation; preferably by the leave one out cross validation, wherein the model is trained n times, and each time one variable is left out from the training dataset, which is of n size minus one variable; and (iv) the predictive model is further validated by predicting new samples using the validation or test dataset, which is the variable left out from the training dataset each time; and the overall performance is the average of all individual performances (n). In another embodiment, the predictive model is optimized by internal cross validation using the k-fold cross validation.In some embodiments, the second step is performed by a classification method wherein the patients are assigned a probability of belonging to given category such as patients having thyroid cancer and at risk of suffering from DM. In some embodiments, the classification method is carried out by a method selected from gradient boosting, support vector machine (SVM), decision trees, K nearest neighbors, Naive Bayes or neural networks. In a preferred embodiment, the classification method is carried out by a Gradient Boosting. As used herein, Gradient Boosting is a machine learning algorithm that uses a gradient boosting framework. Gradient Boosting trees, a decision-tree-based ensemble model, differ fundamentally from conventional statistical techniques that aim to fit a single model using the entire dataset. Such ensemble approach improves performance by combining strengths of models that learn the data by recursive binary splits, such as trees, and of "boosting", an adaptive method for combining several simple (base) models. At each iteration of the gradient boosting algorithm, a subsample of the training data is selected at random (without replacement) from the entire training data set, and then a simple base learner is fitted on each subsample. The final boosted trees model is an additive tree model, constructed by sequentially fitting such base learners on different subsamples. This procedure incorporates randomization, which is known to substantially improve the predictor accuracy and also increase robustness. Additionally, boosted trees can fit complex nonlinear relationships, and automatically handle interaction effects between predictors as addition to other advantages of tree-based methods, such as handling features of different types and accommodating missing data. Hence, in many cases their predictive performance is superior to most traditional modelling methods. In a particular embodiment, the second step is performed by a regression method; preferably selected from multiple linear regression (MLR), principal component regression (PCR), partial least squares regression (PLSR), artificial neural network (ANN), support vector machine (SVM), random forest (RF), lasso regression, ridge regression and combinations thereof.In a particular embodiment, the second step is performed by a classification method, more in particular, by gradient boosting, which includes the value of one or more variables of the methylation profile collected in step (i) and which contribute to the identification of the patients having thyroid cancer and at risk of suffering from DM. In a particular embodiment, the second step is performed by a regression method which includes the value of one or more variables of the methylation profile collected in step (i) and which contribute to the identification of the patients having thyroid cancer and at risk of suffering from DM.

[0226] In another embodiment, the predictive model uses a reference dataset, wherein there are samples of thyroid cancer patients from risk 0 to 1 of suffering from DM, wherein from 0 to 0.5 is indicative of low risk of suffering from DM and from 0.5 to 1 is indicative of high risk of suffering from DM. In a preferred embodiment, the predictive model is a high risk predictive model, wherein a cut-off value is established above which it is indicative of a high risk of suffering from DM and below which it is indicative of a low risk of suffering from DM.

[0227] On a third aspect the invention is directed to a prognostic method. In particular, to an in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:

[0228] a. determining the methylation level of one or more biomarkers as defined in the first aspect of the invention or in any of its preferred embodiments; and

[0229] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by comparing said methylation level / s with a reference value, wherein a deviation in the methylation level of said at least one biomarker selected from the list shown in a) with respect to said reference value, is indicative that the subject is at risk of metastasis, preferably of distant metastasis (DM).Preferably the method of the third aspect comprises:

[0230] a. from an extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 or chr6:85,462,137-85,462, 137 (TBX18 gene), cg18305394 or chr20:56, 134,788-56, 134,788 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6,344,259- 6,344,259; and

[0231] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by (ii) comparing said methylation levels with one or more reference values, wherein a deviation in the methylation level of said biomarkers selected from the list shown in a) with respect to said one or more reference values, is indicative that the subject is at risk of metastasis, preferably of distant metastasis (DM).

[0232] In another preferred embodiment, the method of the third aspect comprises:

[0233] a. from an extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the methylation level of at least all of the biomarkers selected from the group consisting of cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259), cg25901805 (located in PLAGL1 gene, chr6:144, 286, 162-144, 286, 162), cg10780897 (located in RP1 gene, chr8:55, 528, 099-55, 528, 099), cg07530798 (located in chr1:119,535,589- 119,535,589) and cg04690464 (located in LINC02691 gene, chr14:104,711,952-104,711,952); andb. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by (ii) comparing said methylation levels with one or more reference values, wherein a deviation in the methylation level of said biomarkers selected from the list shown in a) with respect to said one or more reference values, is indicative that the subject is at risk of metastasis, preferably of distant metastasis (DM).

[0234] More preferably, the methylation level is determined by a quantitative assay based on bisulfite pyrosequencing to measure the DNA methylation.

[0235] In a first step, the third aspect of the invention comprises the determination of the methylation level of the indicated biomarkers in a biological sample obtained from the subject.

[0236] In some embodiments, the biological sample can come from tissue comprising fresh tissue, frozen tissue and / or formalin-fixed paraffin-embedded tissue (FFPE), and / or body fluids comprising blood, urine, stool, semen, tears, mucus, sweat, milk, cerebrospinal fluid and / or saliva / buccal swabs. In a preferred embodiment, the biological sample is a sample isolated from a subject, such as a biopsy. In some embodiments, the biological sample can contain nucleic acids derived from the thyroid tumor. In a preferred embodiment, the biological sample is a sample containing DNA derived from the thyroid tumor. The thyroid tumor can comprise both the tumor and the tumor microenvironment. In some embodiments, the biological sample is a sample containing DNA derived from normal tissue. In some instances, the normal tissue is non-tumor thyroid tissue adjacent to tumor tissue. In a preferred embodiment, the thyroid tumor is derived from differentiated thyroid cancer (DTC).

[0237] In the second step, the third aspect of the invention comprises comparing said level of methylation with a reference value.Reference value”, as used herein, refers to a laboratory value used as a reference for values / data obtained by laboratory examinations of subjects or samples collected from subjects. The reference value or reference level can be an absolute value; a relative value; a value that has an upper and / or lower limit; a range of values; an average value; a median value, a mean value, or a value as compared to a particular control or baseline value. A reference value can be based on an individual sample value, such as for example, a value obtained from a sample from the subject being tested, but at an earlier point in time or from a non-cancerous tissue. The reference value can be based on a large number of samples, such as from population of subjects of the chronological age matched group, or based on a pool of samples including or excluding the sample to be tested. Various considerations are taken into account when determining the reference value of the marker. Among such considerations are the age, weight, sex, general physical condition of the patient and the like. For example, equal amounts of a group of at least 2, at least 10, at least 100 to preferably more than 1000 subjects, preferably classified according to the foregoing considerations, for example according to various age categories, are taken as the reference group. In another embodiment, the quantity of the biomarker in a sample from a tested subject may be determined directly relative to the reference value (e.g., in terms of increase or decrease, or fold-increase or folddecrease). Advantageously, this may allow to compare the quantity of the biomarker in the sample from the subject with the reference value (in other words to measure the relative quantity of any one or more biomarkers in the sample from the subject vis-a-vis the reference value) without the need to first determine the respective absolute quantities of said biomarker.

[0238] Typically, reference values are the methylation level of the CpG or CG being compared in a reference sample. Thus, in an embodiment, the reference value is the mean methylation level in a pool of samples from healthy subjects, thyroid patients not suffering from DM or thyroid patients suffering from DM.In another embodiment, the reference value for the methylation level of the CpG or CG of interest is the mean level of said CpG or CG in a pool of samples from primary tumours, preferably obtained from subjects suffering from the same type of cancer as the patient object of the study. In a particular embodiment, the reference value is the methylation levels of the CpG or CG of interest in a pool obtained from primary tumor tissue obtained from patients. The methylation profile in the reference sample can preferably, be generated from a population of two or more individuals. The population, for example, can comprise 3, 4, 5, 10, 15, 20, 30, 40, 50 or more individuals. Furthermore, the methylation profile in the reference sample and in the sample of the individual that is going to be prognosticated according to the methods of the present invention can be generated from the same individual, provided that the profiles to be assayed and the reference profile are generated from biological samples taken at different times and are compared to one another. For example, a sample of an individual can be obtained at the beginning of a study period. A reference biomarker profile from this sample can then be compared with the biomarker profiles generated from subsequent samples of the same individual. In a preferred embodiment, the reference sample is a pool of samples from several individuals and corresponds to portions of tissue that are far from the tumor area and which have preferably been obtained in the same biopsy but which do not have any anatomopathological characteristic of tumor tissue.

[0239] Once this reference value is established, the methylation level of this marker in tumor tissue from subjects can be compared with this reference value, and thus be assigned a level of deviation with respect to a reference value. The “deviation” can be either an increase ora decrease in the methylation levels. For example, an increase in methylation level above the reference value of at least 1.1 -fold, 1.5-fold, 2-fold, 5-fold, 10- fold, 20-fold, 30-fold, 40-fold, 50-fold, 60-fold, 70-fold, 80-fold, 90-fold, 100-fold or even more compared with the reference value is considered as “increased” expression level. Similarly, the methylation level is considered increased in a sample of the subject understudy when the levels increase with respect to the reference sample by at least 5%, by at least 10%, by at least 15%, by at least 20%, by at least 25%, by at least 30%, by at least 35%, by at least 40%, by at least 45%, by at least 50%, by at least 55%, by at least 60%, by at least 65%, by at least 70%, by at least 75%, by at least 80%, by at least 85%, by at least 90%, by at least 95%, by at least 100%. On the other hand, a decrease in methylation levels below the reference value of at least 0.9-fold, 0.75- fold, 0.2-fold, 0.1-fold, 0.05-fold, 0.025-fold, 0.02-fold, 0.01 -fold, 0.005-fold or even less compared with reference value is considered as “decreased” methylation level. Similarly, the methylation is considered decreased when its levels decrease with respect to the reference sample by at least 5%, by at least 10%, by at least 15%, by at least 20%, by at least 25%, by at least 30%, by at least 35%, by at least 40%, by at least 45%, by at least 50%, by at least 55%, by at least 60%, by at least 65%, by at least 70%, by at least 75%, by at least 80%, by at least 85%, by at least 90%, by at least 95%, by at least 100% (i.e., absent). The comparison of the methylation levels of the CpG or CG of interest with the reference value allows differential prognosis.

[0240] The methods according to the second or third aspect of the invention allow the determination of the prognosis of a patient suffering from thyroid cancer.

[0241] As it is used herein, the term “prognosis” refers to the prediction in a subject having a thyroid cancer of the likelihood of cancer-attributable death or progression, including metastatic spread. Prognosis may also be referred to in terms of "aggressiveness" or "severity": an aggressive cancer is determined to have a high risk of negative outcome (i.e., negative or poor prognosis) and a non-aggressive cancer has a low risk of negative outcome (i.e., positive or favorable prognosis). An "aggressive" or "severe" tumour is a cell-proliferation disorder that has the biological capability to rapidly spread outside of its primary location or organ.Methods for selecting a therapy for a subject based on the prognostic or determination of risk of the methods according to the invention.

[0242] In another aspect (a fourth aspect), the invention relates to an in vitro method for selecting a subject having thyroid cancer as a candidate to receive an adequate therapy to treat said cancer, the method comprising:

[0243] (i) determining whether the patient is at risk of metastasis, preferably of distant metastasis (DM), following any of the methods according to the invention and (ii) selecting said patient as a candidate to receive an adequate therapy to treat the disease depending on the outcome of the method, in particular depending on whether the patient is at risk or not of developing or having metastasis, preferably distant metastasis (DM). In a first step, the in vitro method for selecting a subject having thyroid cancer as a candidate to receive an adequate therapy to treat said cancer, comprises determining whether the patient is at risk of metastasis or not by using any of the methods according to the invention.

[0244] In a second step, the method comprises selecting a patient to be treated with either a therapy adequate for the treatment of thyroid cancer when the patient is at risk of metastasis, preferably of distant metastasis (DM), or with a therapy adequate for the treatment of thyroid cancer when the patient is not at risk of metastasis, preferably of distant metastasis (DM).

[0245] In some embodiments, when the patient is at risk of metastasis, preferably of distant metastasis (DM), the patient is selected to be treated by a therapy selected from the group consisting of surgery, radiation therapy, chemotherapy, hormonal therapy, or targeted therapy or any combination thereof.

[0246] In some embodiments, the radiation therapy is systemic radiation therapy, in which a patient may swallow or receive an injection of a radioactive substance, such as radioactive iodine or a radioactive substance bound to a monoclonal antibody.Radioactive iodine (131I) is a type of systemic radiation therapy commonly used to help treat cancer, such as thyroid cancer. Thyroid cells naturally take up radioactive iodine. In some embodiments, the chemotherapy includes one or more drugs selected from the group consisting of dacarbazine, docetaxel, doxorubicin, epirubicin, gemcitabine, ifosfamide, paclitaxel, temozolomide, trabectedin and vinorelbine. Conventional chemotherapy has limited efficacy and significant toxicities. Doxorubicin remains the most effective conventional chemotherapeutic agent but has poor response rates. In some embodiments, the chemotherapeutic agent is doxorubicin, which can be administered at 60–75 mg / m2every 3–4 weeks, among other chemotherapy regimens. In some embodiments, therapy is a targeted therapy. The term “targeted therapy”, as used herein, refers to drugs which attack or target proteins that control growth, proliferation and / or survival within cancer cells. Targeted therapies may include, but are not limited to, multikinase inhibitors, also called kinase inhibitors or MKIs, and they normally block the enzymatic activity of kinases. Multi kinase inhibitors that have shown efficacy in the treatment of DTC can be, but are not limited to, lenvatinib, sorafenib, cabozantinib for radioactive iodine (RAI)-refractory DTC who have progressed during previous treatment with sorafenib or lenvatinib, selpercatinib and pralsetinib for DTC harbouring of RET fusion, larotrectinib and entrectinib for tumors with harbouring functional neurotrophic tyrosine receptor kinase (NTRK) fusions (Filetti S, Durante C, Hartl DM, et al. ESMO Clinical Practice Guideline update on the use of systemic therapy in advanced thyroid cancer. Ann Oncol. 2022;33(7):674-684. doi:10.1016 / j.annonc.2022.04.009). Other kinase inhibitors that have been investigated in the clinical trials for radioactive iodine (RAI)-refractory DTC are, but are not limited to, sulfatinib, apatinib, donafenib, pazopanib, drabrafenib, vemurafenib and drabrafenib in combination with trametinib (Zhang L, Feng Q, Wang J, Tan Z, Li Q, Ge M. Molecular basis and targeted therapy in thyroid cancer: Progress and opportunities. Biochim Biophys Acta Rev Cancer. 2023; 1878(4): 188928. doi:10.1016 / j.bbcan.2023.188928).Therapeutic methods coupled to the differential prognostic methods of the invention The invention also provides methods for the treatment of subjects which have been identified as at risk of metastasis, preferably of distant metastasis (DM), by any of the methods according to the invention, wherein if the subject has been identified as a patient at risk of metastasis, preferably of distant metastasis (DM), the subject is treated with a therapy adequate for the treatment of thyroid cancer. In some embodiments, the therapy is selected from the group consisting of surgery, radiation therapy, chemotherapy, hormonal therapy or a replacement therapy.

[0247] The term “treatment”, as used herein comprises any type of therapy, which aims at terminating, preventing, ameliorating and / or reducing the susceptibility to a clinical condition as described herein. In a preferred embodiment, the term treatment relates to prophylactic treatment (i.e. a therapy to reduce the susceptibility of a clinical condition, a disorder or condition as defined herein). Thus, “treatment,” “treating,” and the like, as used herein, refer to obtaining a desired pharmacologic and / or physiologic effect, covering any treatment of a pathological condition or disorder in a mammal, including a human. The effect may be prophylactic in terms of completely or partially preventing a disorder or symptom thereof and / or may be therapeutic in terms of a partial or complete cure for a disorder and / or adverse effect attributable to the disorder. That is, “treatment” includes (1) preventing the disorder from occurring or recurring in a subject, (2) inhibiting the disorder, such as arresting its development, (3) stopping or terminating the disorder or at least symptoms associated therewith, so that the host no longer suffers from the disorder or its symptoms, such as causing regression of the disorder or its symptoms, for example, by restoring or repairing a lost, missing or defective function, or stimulating an inefficient process, or (4) relieving, alleviating, or ameliorating the disorder, or symptoms associated therewith, where ameliorating is used in a broad sense to refer toat least a reduction in the magnitude of a parameter, such as inflammation, pain, and / or immune deficiency.

[0248] The therapeutic method according to the invention is applied to patients which have been diagnosed or determined to be at risk of metastasis, preferably of distant metastasis (DM), by using any of the methods according to the invention. In some embodiments, the method comprises a first step in which the method is applied to the patient and, a second step in which patients determined to be at risk of metastasis, preferably of distant metastasis (DM), are selected and a third step in which the patients are treated with a therapy adequate for the treatment of thyroid cancer.

[0249] Suitable surgical therapies, radiation therapies, chemotherapies, hormonal therapies or targeted therapies have been described in the context of the methods for selecting a therapy for a subject based on the prognostic methods according to the invention and are equally applicable to the therapeutic methods according to the invention.

[0250] Kit of the invention and uses thereof

[0251] In another aspect, the invention relates to a kit, package or device that contains reagents adequate for implementing any of the methods of the invention. It will be understood that, depending on the nature of the method, the reagents adequate for its implementation will vary.

[0252] In the context of the present invention, “kit” is understood as a product containing the different reagents required for carrying out the methods of the invention packaged such that it allows being transported and stored. The materials suitable for the packaging of the components of the kit include glass, plastic (polyethylene, polypropylene, polycarbonate, and the like), bottles, vials, paper, sachets, and the like. Where there are more than one component in a kit they may be packaged together if suitable or the kit will generally contain a second, third or other additional container into which theadditional components may be separately placed. However, in some embodiments, certain combinations of components may be packaged together comprised in one container means. A kit can also include a means for containing any reagent containers in close confinement for commercial sale. Such containers may include injection or blow-molded plastic containers into which the desired vials are retained. One or more compositions of a kit can be lyophilized. In some embodiments, all compositions of a kit of the disclosure will be lyophilized. In some embodiments, a kit of the disclosure with one or more lyophilized agents will be supplied with a re-constitution buffer. Reagents and components of kits may be comprised in one or more suitable container means. A container means may generally comprise at least one vial, test tube, flask, bottle, syringe or other container means, into which a component may be placed, and preferably, suitably aliquoted.

[0253] Furthermore, the kits of the invention can contain instructions for the simultaneous, sequential, or separate use of the different components that are in the kit. Said instructions can be in the form of printed material or in the form of an electronic medium capable of storing instructions such that they can be read by a subject, such as electronic storage media (magnetic disks, tapes, and the like), optical media (CD-ROM, DVD), and the like. The media may additionally or alternatively contain Internet addresses providing said instructions.

[0254] Computer Systems and Devices suitable for carrying out the methods of the invention Methods of the invention can be performed using software, hardware, firmware, hardwiring, or combinations of any of these.Accordingly, in another aspect, the invention relates to a computer-implemented method, wherein the method is any of the methods according to the invention. In another aspect, the invention relates to a computer containing instructions for carrying any of said methods.

[0255] In particular, the present invention relates to a computer implemented method for determining the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer, wherein the method comprises the steps:

[0256] a. receiving the methylation level as described in the first aspect of the invention, said methylation level being determined according to the method of any of previous embodiments, such as bisulfite pyrosequencing;

[0257] b. integrating the methylation levels of step a) above into a predictive model configured to determine the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer; and

[0258] c. generating with the predictive model of step (b) an output indicative of whether the subject is at risk of having metastasis, preferably of distant metastasis (DM).

[0259] In one embodiment, the invention relates to a computer implemented method for determining the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer, wherein the method comprises the steps:

[0260] (i) receiving the methylation level as described in the first aspect of the invention, said methylation level being determined according to the method of any of previous embodiments, such as bisulfite pyrosequencing;

[0261] (ii) correlating the methylation level of one or more of the biomarkers identified in step (i) with representative methylation level profiles from samples obtained from subjects having thyroid cancer previously identified as suffering from metastasis, preferably of distant metastasis (DM), using a predictive model;(iii) associating the subject having thyroid cancer with a methylation level profile according to the correlation of step (ii); and

[0262] (iv) determining the risk of developing metastasis, preferably of distant metastasis of the subject having thyroid cancer as the risk of the profile associated in step (iii);

[0263] wherein preferably the predictive model has been generated by training a computer with at least a plurality of methylation level profiles from previously identified subjects having thyroid cancer and suffering from metastasis, preferably of distant metastasis (DM), by machine learning on said plurality of methylation level profiles so as to obtain representative methylation level profiles associated with metastasis, preferably of distant metastasis (DM).

[0264] In particular, the present invention refers to a computer implemented method for determining the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer, wherein the method comprises all previously mentioned steps, wherein step (i) comprises receiving the methylation levels of one or more biomarkers selected from the group consisting of cg01378044 or chr6:85,462,137-85,462 (TBX18 gene), cg18305394 or chr20:56, 134,788-56, 134,788 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259. In another embodiment, the present invention also refers to a computer implemented method for determining the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer, wherein the method comprises all previously mentioned steps, wherein step (i) comprises receiving the methylation levels of the combination of biomarkers selected from the group consisting of cg01378044 or chr6:85,462,137-85,462 (TBX18 gene), cg18305394 or chr20:56, 134,788-56, 134,788 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259.In a preferred embodiment, the computer implemented method for determining the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer, wherein the method comprises all previously mentioned steps, wherein step (i) comprises receiving the methylation levels of one or more biomarkers selected from the group consisting of cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259), cg25901805 (located in PLAGL1 gene, chr6:144,286,162-144,286,162), cg10780897 (located in RP1 gene, chr8:55,528,099-55,528,099), cg07530798 (located in chr1:119,535,589-119,535,589) and cg04690464 (located in LINC02691 gene, chr14:104,711,952-104,711,952). In another embodiment, the present invention also refers to a computer implemented method for determining the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer, wherein the method comprises all previously mentioned steps, wherein step (i) comprises receiving the methylation levels of the combination of biomarkers selected from the group consisting of cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259), cg25901805 (located in PLAGL1 gene, chr6:144,286,162-144,286,162), cg10780897 (located in RP1 gene, chr8:55,528,099-55,528,099), cg07530798 (located in chr1:119,535,589-119,535,589) and cg04690464 (located in LINC02691 gene, chr14:104,711,952-104,711,952).

[0265] In another embodiment, the invention relates to a computer implemented method for classifying a subject having thyroid cancer as a subject at risk of developing DM, wherein the method comprises the steps:

[0266] (i) receiving the methylation level as described in the first aspect of the invention, said methylation level being determined according to the method of any of previous embodiments, such as bisulfite pyrosequencing;

[0267] (ii) correlating the methylation level of one or more of the biomarkers identified in step (i) with representative methylation level profiles from samples obtained fromsubjects having thyroid cancer previously identified as suffering from metastasis, preferably of distant metastasis (DM), using a predictive model;

[0268] (iii) associating the subject having thyroid cancer with a methylation level profile according to the correlation of step (ii); and

[0269] (iv) classifying the subject having thyroid cancer as a subject with a risk of developing DM as the risk of the profile associated in step (iii);

[0270] wherein the predictive model has been generated by training a computer with a plurality of methylation level profiles from previously identified subjects having thyroid cancer and suffering from metastasis, preferably of distant metastasis (DM), by machine learning on said plurality of methylation level profiles so as to obtain representative methylation level profiles associated with metastasis, preferably of distant metastasis (DM).

[0271] In particular, the present invention refers to a computer implemented method for classifying a subject having thyroid cancer as a subject at risk of developing DM, wherein the method classifies said subject on the basis of the quantitative methylation levels of one or more biomarkers selected from the group consisting of cg01378044 or chr6:85,462,137-85,462 TBX18 gene), cg18305394 or chr20:56, 134,788-56, 134,788 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259. In another embodiment, the present invention also refers to a computer implemented method for classifying a subject having thyroid cancer as a subject at risk of developing DM, wherein the method classifies said subject on the basis of the quantitative methylation levels of the combination of biomarkers selected from the group consisting of cg01378044 or chr6:85,462,137-85,462 (TBX18 gene), cg18305394 or chr20:56, 134,788-56, 134,788 (5’ region of PCK1 gene), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259.In another embodiment, the computer implemented method for classifying a subject having thyroid cancer as a subject at risk of developing DM, wherein the method classifies said subject on the basis of the quantitative methylation levels of one or more biomarkers selected from the group consisting of cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259), cg25901805 (located in PLAGL1 gene, chr6:144, 286, 162-144, 286, 162), cg10780897 (located in RP1 gene, chr8:55, 528, 099-55, 528, 099), cg07530798 (located in chr1:119,535,589-119,535,589) and cg04690464 (located in LINC02691 gene, chr14:104, 711,952-104,711,952).

[0272] It is further noted that preferably, for any of the methods according to the invention the subject has differentiated thyroid cancer (DTC).

[0273] It is further noted that preferably, for any of the methods according to the invention, the subject has papillary thyroid cancer (PTC) or follicular thyroid cancer (FTC).

[0274] Lastly, the present invention further refers to a computer program comprising processor readable instructions which, when the program is executed by a computer, cause the computer to carry out steps of the computer implemented method of the present invention.

[0275] The invention is described below by way of the following examples, which are merely illustrative and do not limit the scope of the invention.

[0276] Examples

[0277] Materials and methods

[0278] Patients, samples and study design

[0279] This was a retrospective study of normal and tumor tissue obtained from non-mDTC and mDTC patients treated at five European centers. Non-mDTC patients were defined asthose with non-aggressive PTC or FTC variants, with no DM at diagnosis, and who were disease-free after initial treatment and remained disease-free during a follow-up of at least five years. mDTC patients were any PTC or FTC patient with either synchronous or metachronous (occurring at least six months after diagnosis) DM. Primary tumors <1 cm in size were excluded. Normal tissue was defined as non-tumor thyroid tissue adjacent to tumor tissue, as indicated by the pathologist. Analyses were performed in a discovery series of 97 samples and in two validation series - one of 13 and one of 45 samples.

[0280] The discovery series for the DNA methylome analysis included 97 formalin-fixed, paraffin-embedded (FFPE) samples obtained from 68 DTC patients who underwent thyroidectomy: 65 primary tumors (30 non-mDTC and 35 mDTC); 15 paired normal tissues; 11 paired DM (mostly from bone or lung, three metachronous); and six paired LNM. DNA methylation was profiled by the Infinium MethylationEPIC array (Illumina).

[0281] Baseline patient and tumor characteristics are shown in Table 2.

[0282] Table 2. Baseline and tumor characteristics of the discovery series.

[0283] Normal Low-risk Lymph node Distant mDTC

[0284] tissue DTC metastases* metastases Patients n n = 15 n = 30 n = 35 n = 6 n = 11 Age (y) mean±SD 52.3±10.7 46±15.3 58.8±14.6 53.9±19.1 63.2±7.4

[0285]

[0286] Female sex n (%) 13 (86.7) 24 (80) 28 (80) 4 (66.7) 8 (72.7) Histology n (%)

[0287] PTC 11 (73.3)B20 (66.7) 21 (60) 5 (83.3)B6 (54.5)BFTC 4 (26.7)B10 (33.3) 14 (40) 1 (16.7)B5 (45.5)BTNM n (%)

[0288] T stage

[0289] T1 8 (26.7) 6 (17.1)

[0290] T2 17 (56.7) 7 (20)

[0291] T3 5 (16.6) 10 (28.6)

[0292] T4 0 (0) 11 (31.4)

[0293] NA 0 (0) 1 (2.9)

[0294] N stage

[0295] Nx-NO 22 (73.3) "4 (40)

[0296] N1 8 (26.3) 20 (57.1)

[0297] NA 0 (0) 1 (2.9)

[0298] DM site n (%)

[0299] Bone 10 (28.6) 8 (72.7)

[0300]

[0301] Lung 9 (25.7) 3 (27.3)Soft tissues 1 (2.8)

[0302] Multiple sites 10 (28.6)

[0303]

[0304] Unknown 5 (14.3)

[0305] DM ocurrence

[0306] Synchronous 32 (91.4) 4 (66.7)c10 (90.1) Metachronous 3 (8.6) 2 (33.3)c0 Unknown 0 0 1 (9.1) Mutations n (%)

[0307] BRAF (unique) 8 (26.6) 1 (2.9) 0 0 RAS (unique) 3 (10) 4 (11.4) 0 1 (9.1) TERT (unique) 0 (0) 10 (28.6) 0 1 (9.1) All wt 15 (50) 7 (20) 0 1 (9.1) Multiple MUT 0 9 (25.7) 0 3 (27.3) At least 1 MUT 2 (6.7) 1 (2.9) 2 (33.3) 3 (27.3)

[0308]

[0309] Unknown 2 (6.7) 3 (8.5) 4 (66.7) 2 (18.1)

[0310] ADerived from mDTC patients.BData corresponding to the histology of the paired primary tumor.cData corresponding to the time of diagnosis of the paired distant metastases.DA unique mutation was identified.EThere is one BRAF mutation other than V600E.FThe three genes were wt.GMore than one mutation was identified, specifically TERT+BRAF or TERT+RAS.HOne mutation was identified and at least data for one mutation was missing. ’All data missing, or at least data for one mutation was missing and the rest of genes were wt. DTC: differentiated thyroid cancer; mDTC: high-risk DTC that develop DM; DM: distant metastasis; MUT: mutation.

[0311] The validation series I included 13 FFPE mDTC (nine from patients with metachronous DM). These samples were analyzed by the Infinium MethylationEPIC array to validate the 156-CpG signature identified in the discovery series. The validation series II included 45 FFPE samples (ten normal tissues, 13 non-mDTC, 11 mDTC, 5 LNM and 6 DM). These samples were analyzed by bisulfite pyrosequencing to validate the adjusted 4-CpG predictive model. Baseline and tumor characteristics of the validation series are shown in Table 3.

[0312] Table 3. Baseline patient and tumor characteristics of the validation series I and II Validation I (n =..... „.

[0313] Io) ' Validation II (n= 45)

[0314] NT (n = non-mDTC mDTC (n MET (n Parameter mDTC (n = 13)

[0315]

[0316] 10) (n = 13) = 11 ) = 11 )Age (years) mean 52.9 ±.?. 66.6 ± 61.2 ± ± SD 57.7 ± 14.6 16.9

[0317]

[0318] 16.9 14Female sex n (%) 8 (61.5) 6 (60) 10 (76.9) 14 (73.7) 8 (72.7) Histology n (%)

[0319] PTC 7 (53.8) 5 (50)a12 (92.3) 8 (72.7) 6 (54.5) FTC 6 (46.2) 5 (50)a1 (7.7) 3 (27.3) 5 (45.5) DM ocurrence

[0320] Synchronous 4 (30.8) 7 (63.6) Metachronous 9 (69.2) 4 (36.4) Unknown 0 (0) 0 (0)

[0321] aData corresponding to the histology of the paired

[0322] primary tumor.

[0323] MET: metastases from mDTC, both lymph node

[0324] and distant metastases

[0325] Cell lines

[0326] All cell lines were maintained at 37°C and 5% CO2in humidified atmosphere. The thyroid

[0327] cancer cell lines TPC1 (PTC-derived) and HTh83 (ATC-derived) were kindly provided by Dr. P. Santisteban, and the BCPAP cell line (PDTC-derived) by Dr. A. Velazquez. TPC1 and HTh83 cells were grown in DMEM (high glucose) and BCPAP in RPMI 1640 (Gibco), both supplemented with 10% of fetal bovine serum, 20pM L-glutamine and 10pM sodium pyruvate. All cell lines were authenticated with the AmpFLSTR® Identifier® Plus PCR Amplification Kit (Applied Biosystems). All cells were tested and found to be Mycoplasma-free during experiments.

[0328] DNA extraction

[0329] For each FFPE sample, 10 non-stained, 10-pm-thick slices (containing at least 80% tumor cells to minimize the effect of contamination by normal cells) were cut from the paraffin blocks and deparaffinized using xylene. Genomic DNA was extracted with E. Z. N. A. FFPE DNA-Kit (Omega Bio-tek) and quantified using a fluorometric method (Qubit, Thermo Fisher Scientific). For cell lines, the AllPrep DNA / RNA / miRNA Universal kit (QIAGEN) was used.Mutation analyses

[0330] BRAF mutations at codon 600 in exon 15, and H-, N-, and K-RAS mutations at codons 12 and 13 in exon 2 and codon 61 in exon 3 were sequenced by Sanger method as follows:

[0331] PTC tissues were screened only for BRAF mutation, whereas follicular variants of PTC with no BRAF mutation and FTC samples were screened for RAS mutations. Matched metastases were analyzed following the same approach.

[0332] TERTp mutations were analyzed by digital PCR analysis using the QX200 Droplet Digital PCR (ddPCR) system (Bio-Rad), according to the manufacturer’s protocol. Specifically, two ddPCR-specific assays, dHsaEXD72405942 and dHsaEXD46675715, were used to detect TERTp C228T and C250T mutations, respectively. Mutant populations were identified and fractional abundance calculations of mutant to total molecules were generated for each sample along with a 95% confidence interval (Cl) using Poisson statistics from the QuantaSoft.

[0333] Infinium MethylationEPIC array

[0334] DNA methylation was analyzed using the Infinium MethylationEPIC v1.0 array. FFPE-derived DNA samples were checked for their eligibility as indicated by the qPCR-based Infinium HD FFPE QC Assay (Illumina). All samples qualified for restoration, with ACq values <5 relative to intact genomic DNA control. The amount of 300 ng of DNA from FFPE samples and 600 ng of DNA from cell lines was used for bisulfite conversion with the EZ-96 DNA Methylation kit (Zymo Research), following the manufacturer’s recommendations for Infinium methylation profiling assays. Bisulfite-converted FFPE DNA was used for restoration with the Infinium FFPE Restoration kit (Illumina) using the Infinium HD FFPE Restore Protocol supplied by the manufacturer. The resulting materialwas used as input for hybridization on methylation arrays, using the manufacturer’s recommended incubation, centrifugation and hybridization oven setup and following the automated processing protocol on a Tecan EVO 150 robot for denaturing, hybridization and staining steps. Arrays were scanned on a HiScan (Illumina) with HiScan Control Software. Data quality was first evaluated with Genome Studio Methylation module.

[0335] Analysis of DNA methylation data

[0336] Data from Infinium Human MethylationEPIC array were analyzed using the RnBeads version 1.2.2 R package (Assenov Y, et al. Comprehensive analysis of DNA methylation data with RnBeads. Nat Methods. 2014; 11 (11): 1138–1140). For each CpG site, methylation levels were quantified by beta values ranging from 0 to 1, with 0 indicating completely unmethylated and 1 indicating completely methylated. To rule out technical and biological biases, we removed a subset of probes that did not pass filtering criteria: (i) 139,816 probes overlapping with known single-nucleotide polymorphisms (SNPs); (ii) 50,620 probes with low signal (detection P > 0.05, determined by the Greedycut algorithm); (iii) 1,720 non-CpG targeting probes; and (iv) 16,073 probes mapping in ChrX and ChrY. We also excluded two samples with <99% of probes with detection P > 0.05 as determined by the Greedycut algorithm. Beta values were normalized using the BMIQ (Beta Mixture Quantile dilation) method to adjust the beta values of Type II probes into a statistical distribution characteristic of Type I probes (Teschendorff AE, et al. A betamixture quantile normalization method for correcting probe design bias in Illumina Infinium 450 k DNA methylation data. Bioinformatics. 2013;29(2): 189-96).

[0337] Global DNA methylation was calculated as the mean beta value across the CpGs covered by EPIC array after quality control for each sample. When indicated, global DNA methylation was calculated for specific genomic compartments of interest: CGIs, CGI shores and non-CGI regions.Differential DNA methylation at the site level was analyzed through hierarchical linear models as implemented in the limma package (Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3: Article3) and using M-values - defined as log2(β-value) / (1 - β value) - which exhibit a distribution that is more consistent with limma's statistical model assumptions than the beta values (Du P, et al. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinformatics. 2010;11(1):587). P-values were corrected for multiple testing using a false discovery rate (FDR). Probes with adjusted P < 0.05 and delta-beta (Δβ) > 0.2 or < -0.2 between groups of samples were considered hyper- or hypomethylated, respectively.

[0338] Genomic annotation of CpGs

[0339] Annotation was done using data from the Illumina MethylationEPIC Manifest. For the location relative to a gene, we used the following categories: 5' region (TSS 1500: 201-1500 bp upstream of the transcriptional start site [TSS]); TSS 200 (1-200 bp upstream of the TSS, 5' UTR and first exon); gene body (from the first intron to the last exon and 3' UTR); and intergenic regions. When both 5’ region and gene body were listed as relative positions to a gene, priority was given to 5’ region. Some CpGs located in regions overlapping with more than one gene or at the same distance from different genes were assigned multiple gene annotations. For the location relative to a CGI, we used three categories: CGI; CGI shore (0-2 kb from the CGI edge); and non-CGI (>2 kb from the CGI edge). Repetitive sequences for the hg38 version of the human genome were downloaded using the R package Annotationhub (version 3.2.2) (Morgan Martin, Lori Shepherd. AnnotationHub: Client to access AnnotationHub resources. R package version 3.2.2. https: / / bioconductor.org / packages / AnnotationHub. 2022).We annotated potential enhancers based on published open chromatin data from three sources: the ENCODE Consortium (ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489(7414):57-74) downloaded through the University of California, Santa Cruz (UCSC) genome browser (http: / / hgdownload.soe.ucsc.edu / gbdb / hg38 / encode3 / ccre / ); enhancer data from the FANTOM5 project (http: / / fantom.gsc.riken.jp / 5 / datafiles / latest / extra / Enhancers / ) (Andersson R, et al. An atlas of active enhancers across human cell types and tissues. Nature. 2014;507(7493):455-461); and predicted enhancer-gene maps based on the activity-by-contact (ABC) model in thyroid tissue (Nasser J, et al. Genome-wide enhancer maps link risk variants to disease genes. Nature. 2021;593(7858):238-243).

[0340] Chromatin state enrichment

[0341] To study chromatin state enrichment of the CpGs that were differentially methylated during DTC progression, we downloaded data from 18 cell types (since there is no model for thyroid samples) using the 18-state model of the NIH Roadmap Epigenomics Consortium (Kundaje A, et al. Integrative analysis of 111 reference human epigenomes. Nature. 2015;518(7539):317-330).

[0342] Regions with chromatin states 1-4 (active TSS, flanking TSS, flanking TSS upstream, and flanking TSS downstream) were considered " TSS-proximal regions"; states 5-6 (strong and weak transcription) "transcription-related states"; states 7-11 (genic, active and weak enhancers) “enhancers”; state 12 " ZFN genes and repeats"; state 13 “heterochromatin”; states 14-15 (bivalent / poised TSS and bivalent enhancer) "bivalent promoters and enhancers"; states 16-17 (repressed Polycomb or weak repressed Polycomb) " Polycomb-repressed regions"; and state 18 (quiescent / low) "quiescent region". We estimated the ratio of observed to expected differentially methylated CpGs in the different groups of tumors (non-mPTC, mPTC, PTC-DM, non-mFTC, mFTC, andFTC-DM) to normal tissue, using the CpGs included in the Infinium MethylationEPIC array as background.

[0343] TCGA dataset on DNA methylation in thyroid tumors

[0344] We used DNA methylation data from TCGA thyroid dataset (The Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell. 2014;159(3):676-690) generated by Infinium HumanMethylation450 BeadChip (covering 450,000 CpGs), including 56 normal, 110 low-risk PTC (tumor size of T1 / T2, NO or N, and resection of R0 / R1; excluding tall cell or columnar subtypes) and 17 high-risk PTC (tumor size of T4 / T4a / T4b or with DM) fresh-frozen tissues, as classified by the American Thyroid Association’s (ATA) risk stratification, and excluding microcarcinomas and intermediate-risk PTC. Baseline patient and tumor characteristics for this dataset are shown in Table 4.

[0345] Table 4. Baseline patient and tumor characteristics of The Cancer Genome Atlas (TCGA) dataset used in this study.

[0346] The Cancer Genome Atlas - TCGA (n = )a

[0347]

[0348] Parameter NT (n = 56) Low-risk (n = 110) High-risk (n = 17)

[0349] Age (years) mean ± SD 45.66 ± 16.87 46.94 ± 16.29 62.24 ± 11

[0350] Female sex n (%) 36 (64.29) 86 (78.18) 9 (52.94)

[0351] Histology n (%)

[0352] Classical PTC - 71 (64.55) 13 (76.47) FVPTC - 39 (35.45) 4 (23.53)

[0353] Mutation n (%)

[0354] BRAFv600E- 49 (44.55) 13 (76.47) RAS - 22 (20) 3 (17.67) TERTp 6 (5.45) 10 (58.8)aRisk defined according to the American Thyroid Association’s (ATA) classification

[0355] (TCGA, Cell 2014).

[0356] Low-risk PTC: tumor size of T1 / T2, NO or N, and resection of R0 / R1; excluding tall cell

[0357] or columnar subtypes)

[0358] High-risk PTC: tumor size of T4 / T4a / T4b or with DM.Microcarcinomas and intermediate-risk PTC were excluded. Samples with availability of

[0359] DNA methylation data.

[0360] Methylation data (expressed in beta values) were downloaded through the XENA browser (Goldman MJ, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. 2020;38(6):675-678). Differential DNA methylation analysis was performed following the same procedure as with data from the discovery series.

[0361] TCGA dataset on gene expression in thyroid tumors

[0362] We used gene expression data from TCGA thyroid dataset (The Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell. 2014;159(3):676-690) generated by RNA-Seq, including the same samples analyzed for DNA methylation (Table 4). Raw count expression data were downloaded from the Firebrowse portal. Differential expression analysis of TCGA raw counts was performed using DESeq2 (version 1.34.0) R package (Love MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.

[0363] 2014;15(12):550), removing genes with low expression across all samples (<10 counts). Differentially expressed genes were defined as those with an adjusted P < 0.05 (Benjamini & Hochberg correction) and a |logFC| > 1.

[0364] Cluster analysis

[0365] Heatmaps were created using the pheatmap R package (version 1.0.12). The clustering distance was determined using the Minkowski method, and the clustering scale with the Ward.D2 method.Functional enrichment

[0366] Functional enrichment was analyzed with the GO analysis tool from the GO Consortium website (geneontology.org) (Ashburner M, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25-29). We used the BP, MF and CC terms annotated in the GO database. Statistical significance was defined as FDR < 0.05. This analysis was applied to the genes associated with CpGs located in the 5’ region (according to Illumina’s EPIC) that were differentially methylated in non-mPTC, non-mFTC, mPTC, mFTC, PTC-DM and FTC-DM vs normal tissue. To reinforce the results, GO analysis was also applied to the genes differentially expressed in low-risk and high-risk PTC vs normal tissue in TCGA dataset that had CpGs located in their 5’ regions that were differentially methylated in low-risk PTC and mPTC vs normal tissue in our discovery series. Finally, GO analysis was also applied to the genes associated with the 156-CpG signature identified in our discovery series.

[0367] Bisulfite sequencing and pyrosequencing

[0368] DNA bisulfite treatment was performed using the EZ DNA Methylation Kit (Zymo Research) according to the manufacturer’s instructions, and bisulfite-treated DNA was amplified by nested PCR using IMMOLASE hot-start DNA Polymerase (Bioline). Each CpG was assessed twice in each sample and amplicons were pooled to ensure a representative methylation profile. DNA methylation of individual CpGs within the ELTD 1 -associated CGI was analyzed by Sanger sequencing, while methylation of the CpGs in the 4-CpG predictive model was analyzed by pyrosequencing. Direct sequencing was performed using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) by Eurofins Genomics. The degree of methylation was calculated by comparing the peak of the cytosine residues with the peak of the thymine residues, as previously described (Melki J, et al. Cancer-specific region of hypermethylation identified within the HIC1 putative tumour suppressor gene in acute myeloid leukaemia.Leukemia. 1999;13(6):877-836883). Pyrosequencing assays were performed with the PyroMark Q24 sequencer (Qiagen, Hilden, Germany), as previously described (Tost J, Gut IG. DNA methylation analysis by pyrosequencing. Nat Protoc. 2007;2(9):2265-75) and DNA methylation levels were quantified with the PyroMark Q24 software (Qiagen). DNA methylation ranged from 0 to 1, with 0 indicating completely unmethylated and 1 indicating completely methylated. The primers and sequences are shown in Table 5.

[0369] Table 5. Primers and sequences used for the bisulfite sequencing

[0370] and pyrosequencing assays.

[0371] Annealing Amplicon Primer nameaSequence

[0372] temperature (°C) length (bp)

[0373]

[0374] cl Trid l_. A AAGGGA GG ELTD1_Fwext _

[0375]

[0376] Tri D. ATCAATAATCACTAAACA ELTD1_Rvext..

[0377] r r

[0378] . TTGAGTATTTATTTAGGA

[0379] ELTD1 Fw int

[0380] G 50 199 ELTD1_Rv int CACACACTCACACAAAAA

[0381] ■- * TTTATGATGTAATAAGTT

[0382] cg01378044 Fw ext

[0383] AATT 48 119 cg01378044 Rv ext CAACCTTAATTCCTACCA

[0384] 5' BIO- cg01378044 Fw int TTTATGATGTAATAAGTT

[0385] AATT 48 103 cg01378044 Rv int CAA„TTTAAAAAAAATCAAT

[0386] Ax- /

[0387] Pyrosequencing CA TTTAAAAAAAA TCAA T

[0388] primer AC

[0389] cg18305394 Fw ext G. A.TAA. GATATAATTGGTTG

[0390] A50 125 cg18305394 Rv ext C. T „ATAACAAATACAAATC

[0391] AC 1

[0392] 5' BIO- cg18305394 Fw int GAAAGATATAATTGGTTG

[0393] AATA 52 111 cg18305394 Rv - in *t T. T. A „AATCACACTATATTA

[0394] AAv

[0395] Pyrosequencing TTAAA TCACAC TA TA TTA

[0396] primer AAC

[0397] cg14672100 Fw ext AAGGGTTAGTAGAAGTTAGT

[0398] cg14672100 Rv ext C. A. AAAATAAACCCTACAA50 120

[0399] r r5' BIO- cg14672100 Fw int AGGTAGTAGAAGTTAGT

[0400] AGT 52 84 cg d14xe6-7r2od1A0A0 D Rv - in *t ATC. T. A

[0401] I AArATAACTAAAAATAA w Pyrosequencing A C TA A TA A C TA AAA A TA A p

[0402]

[0403] rimer TAAC

[0404] cg11925561 Fw ext TTTATTTGGTGTGGGGTT

[0405] * AAAACAAACAACAAAACC 50 129 cg11925561 Rv extA

[0406] cg 4141A 92O5C5CA641 r Frw ■ intTGTTGG I I I I I uuA I I GG I I

[0407] 5' BIO- 52 98 cg11925561 Rv int AAAACAAACAACAAAACC

[0408] A

[0409] Pyrosequencing G l GG 1111 IGGATT GG TT

[0410] primer TT

[0411] aThe assays are based on

[0412] nested-PCRs.

[0413] ext: external PCR; int:

[0414] internal PCR

[0415] Predictive model for risk of DM

[0416] We identified a CpG-based predictive model for the risk of developing DM. Specifically, the 156 CpGs differentially methylated between non-m DTC and mDTC were reduced to 34 CpGs by applying the SimpleLogistic algorithm from WEKA (Frank E, et al. The WEKA Workbench. Online Appendix for “Data Mining: Practical Machine Learning Tools and Techniques”. Morgan Kaufmann Publishers; 2016) to the beta values obtained with EPIC arrays. To avoid overfitting, we selected the ten CpGs with the highest size effect, as the quality of the classification was not affected. These ten CpGs were still highly correlated. We then applied a forward stepwise logistic regression model to select the best fit by Akaike information criterion (AIC) and identified four CpGs as the best variables for the predictive model. The performance of the model was estimated by leave-one-out cross-validation. The predictive model was then adjusted to bisulfite pyrosequencing values with a generalized linear model for binomial distributions (glm function in R-base).Statistics

[0417] All analyses were performed using R v4.1.3. Values are expressed as means ± standard deviation (SD). All correlations were evaluated by the Pearson correlation coefficient (r). Normality of data was checked with the Shapiro- Wilk test, and data were analyzed with the Mann-Whitney U or the Kruskal-Wallis test, as appropriate, followed by Dunn's post-hoc test when needed, using the R package dunn.test (version 1.3.5). Significance was set at P < 0.05. ROC curves were created using the R packages OptimalCutPoints (version 1.1-5) (Lopez-Raton M, et al. OptimalCutpoints: An R Package for Selecting Optimal Cutpoints in Diagnostic Tests. J StatSoftw. 2014;61(8)) and verification (version 1.42) (NCAR - Research Applications Laboratory, verification: Weather Forecast Verification Utilities, https: / / cran.r-project.org / web / packages / verification / index.html.

[0418] 2015).

[0419] Results

[0420] Characterization of the discovery series

[0421] The discovery series consisted of a total of 97 samples: 15 normal thyroid tissue, 65 primary 5 DTC tumors (30 non-mDTC and 35 mDTC), 6 LMN, and 11 DM. Among the 65 primary tumors, 38 (58.5%) harbored at least one mutation in BRAF, RAS or the promoter of TERT (TERTp) (Table 6).

[0422] Table 6. Mutations detected in primary tumors and metastases of differentiated thyroid cancer (DTC) classified according to tumor histology and risk

[0423] Lymph node Distant non-mDTC mDTC * *t.

[0424] metastases metastases Mutations PTC n FTC PTC n FTC n PTC FTC PTC FTC

[0425]

[0426] n (%) = 20 = 10 = 21 = 14 n = 5 n = 1 n = 6 n = 5

[0427] BRAF. 11 Jnn 13 14M(55)(10° (61.9) (100) COO) (50) (100)

[0428]

[0429] 4

[0430] BRAF(V60 7 (35) 0 0 0 0 0 0 OE)(19 0)

[0431] otherA / rr\ 0 0 0 0 BRAF MUT1 (5>0 1 (4 8) 0

[0432] 3

[0433] Unknown 1 (5) 0 05o30

[0434]

[0435] (100) (50)

[0436] RAS

[0437] wt 14 9 12 812 ( 0 0

[0438] (70) (90) 57.2) (57.1)<1®7(40) MUT 4 5

[0439] 3<15> (10) (19.0) ( 0 0 2

[0440] 35.7) <33 3(40) Unknown 3 (15) 0 5 1 5 1 3 1 (23.8) (7.2) (100) (100) (50) (20) TERT

[0441] 1R2wt1010 4 3 1 (90) (1°° (47.6) (28.6) (80)0(33'3(20)

[0442] 11 10

[0443] MUT 1 (5)a0 (52.4) (71.4) 1 1 3 4 b c (20) (100) (50) (80)

[0444] 1 Unknown 1 (5) 0 0 (0) 0 (0)10 (16.7 0 (20)

[0445] aMutant allelic

[0446] frequency < 5%

[0447] bOne tumor with mutant

[0448] allelic frequency < 5%

[0449] cTwo tumors with mutant

[0450] amallelic frequency < 5%

[0451] BRAFV600E mutations were detected in 11 tumors, RAS mutations in 13, and TERTp mutations in 22. Concomitant RAS and TERTp mutations were detected in five tumors, and concomitant BRAF and TERTp mutations in four. TERTp mutations have been associated with poor outcomes in thyroid cancer patients (Liu R, Xing M. TERT promoter mutations in thyroid cancer. Endocr Relat Cancer. 2016;23(3): R143-R155; Park J, et al. TERT Promoter Mutations and the 8th Edition TNM Classification in Predicting the Survival of Thyroid Cancer Patients. Cancers (Basel). 2021;13(4):648), and as would be expected, we identified this mutation only in mDTC (with the sole exception of one nonmetastatic primary PTC [non-mPTC] with a TERTp mutant allele frequency of 2.9%). Interestingly, LNM and DM harbored the same mutations as their paired primary tumors.

[0452] Global DNA methylation changes along PTC and FTC metastatic progression

[0453] DNA methylation levels have been shown to be dependent on genomic context (Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13(7):484-92.; Baubec T, Schubeler D. Genomic patterns and context specific interpretation of DNA methylation. Curr Opin Genet Dev. 2014;25:85-92) and, although the EPIC array analysis showed no significant differences in global DNA methylation levels between normal tissue and tumors (p = 0.09), we found that CGIs were globally hypermethylated in mDTC compared to normal tissue and non-mDTC, while non-CGI regions were globally hypomethylated in mDTC compared to normal tissue (Figure 1). CGI shore did not show significant differences. When PTC and FTC were analyzed separately, this methylation pattern was maintained, except in non-mPTC, which also showed global hypomethylation of non-CGI regions (Figure 9). Global methylation of LNM and DM was similar to that of metastatic primary tumors, but differences were not statistically significant in most cases, probably due to the relatively small sample size of LMN and DM. Since EPIC array mostly covers unique sequences, these results revealed that global hypomethylation did not exclusively affect Alu elements, as we had previously reported (Klein Hesselink EN, et al. Increased Global DNA Hypomethylation in Distant Metastatic and Dedifferentiated Thyroid Cancer. J Clin Endocrinol Metab. 2018; 103(2): 397-406), but also affected unique sequences (Table 7, Figure 10).

[0454] Table 7. Distribution of differentially methylated CpGs according to their location in unique or repetitive sequences

[0455]

[0456]

[0457] Differentially methylated CpGs in tumor vs normal tissue

[0458]

[0459] Hypormothyl.ition

[0460] Hypomethylations CGI CGI-shore non-CGI Repetiti Repetiti Repeti Repetit Repetiti Unique ve Unique ve Unique tive Unique ive Unique ve no 78

[0461] n- 88 11 79 21 98 1. 85 14.8 21 mP 47.4.6 45.0 12.0.8 2 53.9 8.1 44 0 11.2 TC 8 % 63 % 81 % 19 % 82 % 1 % 7 % 8 % 40 % 93 %

[0462] 88

[0463] 91 8. 75 24 96 3. 94 5..6 11 mP 19.7 18 3 32.3 10.7 43.4 1 6 33.7 1 3 11 0 14.4 TC 95 % 0 % 50 % 66 % 2 % 6 % 8 % 9 % 25 % 5 %

[0464] 72

[0465] PT 91 8. 19 71 28 98 1. 84 7 16 19.3 27 C- 77.3 74 7 06.4 76.6 30.2 5 8 39.0 5.0 79 0 75.7 DM 27 % 0 % 6 % 49 % 50 % 7 % 48 % 4 % 5 % 78 % no 81

[0466] n- 93 6. 77 22 98 1. 89 10.5 18 mF 14.1 11 9 13.5 38.5 38.5 5 56.8 6.2 18 0 42.5 TC 95 % 0 % 30 % 6 % 2 % 6 % 5 % 4 % 78 % 6 %

[0467] 72

[0468] 93 6. 71 28 97 2. 83 4 17.9 27 mF 37.2 27 8 93.7 36.3 15.9 3 1 21.0 4.0 94 0 34.1 TC 11 % 1 % 83 % 97 % 32 % 3 % 60 % 3 % 02 % 92 %

[0469] 76

[0470] FTC 90 9. 72 27 98 1. 86 4 13.6 23 73.9 73 1 58.7 21.3 23.2 4 8 27.7 1.3 80 0 24.4 DM 78 % 8 % 17 % 88 % 92 % 4 % 41 % 9 % 62 % 63 %

[0471] 2

[0472] EPI 78 21 78 21 96 5 3. 81 8 18 71 28 C 67.2 18.7 67.2 18.7 15.8 1 1 12.3 8.6 39.8 15.1 arr 72 1 86 9 72 1 86 9 63 4 0 6 56 2 7 8 52 7 47 3

[0473]

[0474] ay 04 % 76 % 04 % 76 % 23 % 3 % 67 % 0 % 14 % 03 %

[0475] In summary, metastatic primary PTC (mPTC), metastatic primary FTC (mFTC), and non-mPTC - but not non-metastatic primary FTC (non-mFTC) -are all characterized by global hypomethylation of non-CGI regions, suggesting that this is an early event in PTC - but not FTC - carcinogenesis. In contrast, only mPTC and mFTC are characterized by global hypermethylation of CGIs, suggesting an association with progression and metastases.

[0476] Progressive increase of DNA methylation alterations along PTC and FTC metastatic progression

[0477] We then analyzed single CpGs to better understand the dynamics of DNA methylation during DTC metastatic progression (Figure 11) and observed a progressive increase of DNA methylation alterations from non-mDTC to DM, with a predominance ofhypomethylation (Figure 12A), in both PTC and FTC (Figure 2A). Compared to normal tissue, in non-mPTC, mPTC, and DM derived from PTC (PTC-DM), we identified 6,343, 6,493 and 35,191 differentially methylated CpGs, respectively. In non-mFTC, mFTC and DM from FTC (FTC-DM), we identified 3,323, 17,069 and 16,125 differentially methylated CpGs, respectively. We also analyzed LNM derived from mPTC and found an increase of DNA methylation alterations compared to primary tumor but less than in DM (Figure 12B). We found a high correlation between these results and PTC data from The Cancer Genome Atlas (TCGA), highlighting the robustness of our DNA methylome analysis (Figure 13).

[0478] In line with our findings on global DNA methylation, 91% of CpGs in non-mPTC and 52% of CpGs in non-mFTC were hypomethylated in comparison with normal tissue, while 67% of CpGs in mPTC and 77% of CpGs in mFTC were hypomethylated in comparison with normal tissue (Figure 12B). We also observed a high percentage of hypomethylated CpGs in DM, especially in PTC-DM. Differences in CpG methylation were greater between DM and normal tissue than between primary tumors (both metastatic and non-metastatic) and normal tissue (Figure 14).

[0479] More importantly, a subset of the methylation alterations in non-mPTC and non-mFTC (32.2% and 46.9%, respectively) were maintained in mPTC and PTC-DM or, mFTC and FTC-DM, respectively, suggesting a linear progression of these alterations. In addition, a significant proportion of the de novo alterations found in mPTC (49.6%) and mFTC (34.7%) 200 were maintained in DM (Figure 2A, Figure 15), including in three metachronous DM, suggesting that the DM were seeded from a major subclone of the primary tumor before the initial diagnosis of DTC.

[0480] We also compared methylation alterations between PTC and FTC at different points of metastatic progression. We found a low overlap of methylation alterations between non-mPTC and non-mFTC, while the overlap between mPTC and mFTC and between FTC-DM and PTC-DM was higher (Figure 2B).Taken together, these results showed a progressive gain of methylation alterations -mostly hypomethylation - during metastatic progression in PTC and FTC, leading us to postulate histology-specific epigenetic alterations in the initiation of PTC and FTC but common epigenetic alterations in the metastatic progression of both histologies.

[0481] Genomic distribution and functional characterization of DNA methylation alterations in DTC metastatic progression

[0482] To better understand the functional implications of methylation changes during metastatic progression, we analyzed the genome distribution of the differentially methylated CpGs (Figure 2C, Table 8). Tumor-associated hypomethylation events in PTC and FTC showed a high enrichment of CpGs far from CGIs at intergenic regions (all P < 2.2e-16, odds ratio [OR] > 3 and OR > 1.1, respectively), and depletion of CpGs within CGIs at 5’ regions (all P < 2.2e-16, OR < 0.07 and OR < 0.5, respectively).

[0483] Table 8. Genomic distribution of differentially methylated CpGs

[0484] Hypomethylated CpGs in tumors vs NT

[0485] Odd ratios (ORs) respect to EPIC array

[0486] Nr CpGs CGI Shore non-CGI

[0487] non-mPTC 5800 0.5049428 5.10

[0488] mPTC 4316 0.01909322 0.5077165 5.15825

[0489]

[0490] PTC-DM 26718 0.05547585 0.4965457 4.808174

[0491] non-mFTC 1716 0.07288878 0.7678447 3.132611

[0492] mFTC 13080 0.03888646 0.6600674 3.841563

[0493]

[0494] FTC-DM 8005 0.03717858 0.6851253 3.721811

[0495] P-values (ORs) respect to EPIC array

[0496] Nr CpGs CGI Shore non-CGI

[0497] non-mPTC 5800 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0498] mPTC 4316 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0499]

[0500] PTC-DM 26718 < 2.2e-16 < 2.2e-16 < 2.2e-16non-mFTC 1716 < 2.2e-16 9.06E-02 < 2.2e-16 mFTC 13080 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0501]

[0502] FTC-DM 8005 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0503] Hypermethylated CpGs in tumors vs NT

[0504]

[0505] Odd ratios (ORs) respect to EPIC array Nr CpGs CGI Shore non-CGI non-mPTC 541 0.3897293 0.4672445 2.748282 mPTC 2175 1.062641 0.9010225 1.021828

[0506]

[0507] PTC-DM 8467 1.970526 1.506648 0.449285

[0508] non-mFTC 1605 1.175975 1.438078 0.700104 mFTC 3982 2.325567 1.453578 0.3957734

[0509]

[0510] FTC-DM 8116 1.683548 1.622295 0.4863325

[0511] P-values (ORs) respect to EPIC array Nr CpGs CGI Shore non-CGI non-mPTC 541 3.61E-09 2.19E-05 < 2.2e-16 mPTC 2175 0.2561 0.07293 0.6414

[0512]

[0513] PTC-DM 8467 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0514] non-mFTC 1605 0.007481 1.85E-06 1.64E-09 mFTC 3982 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0515]

[0516] FTC-DM 8116 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0517] Hypomethylated CpGs in tumors vs NT

[0518] Odd ratios (ORs) respect to EPIC array Nr CpGs 5' region Body Intergenic non-mPTC 5800 0.51 1.26163 1.42 mPTC 4316 0.5017078 0.8997507 2.021478

[0519]

[0520] PTC-DM 26718 0.4499802 0.7577117 2.528985

[0521] non-mFTC 1716 0.5426821 1.483267 1.131368 mFTC 13080 0.5094272 0.8507028 2.107453

[0522]

[0523] FTC-DM 8005 0.4967928 0.9533145 1.926668

[0524] P-values (ORs) respect to EPIC array Nr CpGs 5' region Body Intergenic non-mPTC 5800 < 2.2e-16 < 2.2e-16 < 2.2e-16 mPTC 4316 < 2.2e-16 0.0009112 < 2.2e-16

[0525]

[0526] PTC-DM 26718 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0527] non-mFTC | 1716 | < 2.2e-16 | 6.26E-13 | 0.02019 |mFTC 13080 < 2.2e-16 < 2.2e-16 < 2.2e-16

[0528]

[0529] FTC-DM 8005 < 2.2e-16 0.04054 < 2.2e-16

[0530] Hypermethylated CpGs in tumors vs NT

[0531] Odd ratios (ORs) respect to EPIC array

[0532] Nr CpGs

[0533] non-mPTC 541 0.433187 2,2.153723 0.8506042

[0534] mPTC 2175 0.7728413 1.33102 0.9311054

[0535]

[0536] PTC-DM 8467 1.114824 0.8909831 1.010149

[0537] non-mFTC 1605 1.005474 0.8267779 1.225832

[0538] mFTC 3982 1.046344 0.8864464 1.09095

[0539]

[0540] FTC-DM 8116 1.059214 0.847558 1.13027

[0541] P-values (ORs) respect to EPIC array

[0542] Nr CpGs 5' region Body Intergenic

[0543] non-mPTC 541 4.74E-13 < 2.2e-16 0.1127

[0544] mPTC 2175 4.02E-05 4.90E-08 0.1499

[0545]

[0546] PTC-DM 8467 1.89E-03 4.05E-04 0.6778

[0547] non-mFTC 1605 0.916 0.0002861 0.0001811

[0548] mFTC 3982 0.1745 0.000281 0.01292

[0549]

[0550] FTC-DM 8116 0.01427 1.40E-09 6.02E-04

[0551] Data was checked for normality using Shappiro. p-values were obtained using Fisher's t-test. Significant ORs (p-value < 0.05) are indicated in green if associated with depletion and red if associated with enrichment.

[0552] Interestingly, hypomethylated CpGs in non-mPTC and non-mFTC showed an enrichment of CpGs at gene bodies (all P < 2.2e-16, OR > 1.3). In FTC, hypermethylated CpGs were enriched in CGIs (P < 0.0075, OR > 1.2) and depleted in non-CGI regions (P < 1.64e-09, OR < 0.7). This pattern was also observed in PTC-DM (P < 2.2e-16, OR > 1.5 and OR < 5, respectively). In contrast, in non-mPTC, hypermethylated CpGs were depleted in CGIs (P < 3.6e-09, OR < 0.4) and enriched in non-CGI regions (P < 2.2e-16, OR > 2.7). However, there was no clear pattern of hypermethylation events in regard to the relative position of the nearest gene. The differences in distribution between hypomethylated and hypermethylated CpGs may be due to different regulatory mechanisms. Interestingly, there were differences in the genomic locations ofhypermethylated CpGs between non-mPTC and non-mFTC, again suggesting that many of the epigenetic alterations involved in PTC and FTC initiation may be histology-specific. We then analyzed chromatin state enrichment (Figure 3, Figure 16). Hypomethylated CpGs were enriched in enhancers, and CpGs in 5' regions were also enriched in TSS-proximal states, suggesting a potential regulatory role. However, there was a progressive decrease in the enrichment of these states during metastatic progression, with a greater decrease seen in PTC. Hypermethylated CpGs were enriched in bivalent promoters and enhancers, characterized by the coexistence of active and repressive histone marks, as well as in Polycomb repressive domains. We also observed enrichment in TSS-proximal states, especially for CpGs in 5' regions, and to a lesser extent, in enhancers and zinc finger protein genes, especially for CpGs in gene bodies. This pattern was less evident in non-mPTC, although PTC displayed an increasing enrichment of these states during progression. Neither hypomethylated nor hypermethylated CpGs were enriched in transcription-related states or constitutive heterochromatin or quiescent regions (associated with very low levels of histone marks). These results provide evidence for the differential role of hypomethylation and hypermethylation events during metastatic progression in PTC and FTC.

[0553] To further investigate the potential functional implications of these alterations in DNA methylation, we analyzed gene ontology (GO), focusing on the differentially methylated CpGs located in 5’ regions and their associated genes.

[0554] In all tumor groups, we found high enrichment of GO biological process (BP) terms associated with cell adhesion, which has a key role in metastatic dissemination (Figure 4, Tables 9-11).

[0555] Table 9. Gene ontology (GO) enrichment analysis of the genes associated with CpGs differentially methylated at their 5' region. GO term: biological process

[0556]

[0557] PTCgenesEPI Cou Ex Fold

[0558] C - nt pe Enric P- GO term: biological REFLIST (92 cte hmen val FD

[0559] process (17968) 3) d t ue R Reduced Term n 0.0

[0560] o regulation of plasma 00 0.0 cellular component Fimembrane organization 0.8 10 31 organization or biogenesis rn (GO:1903729) 17 7 7 8.02 8 8 (membrane)

[0561] P 0.0

[0562] T negative regulation of actin 00 0.0 cellular component C filament depolymerization 2.1 15 40 organization or biogenesis V (G0:0030835) 41 10 1 4.75 2 8 (cytoskeleton) s negative regulation of actin 2.5 0.0 cellular component N filament polymerization 2.8 2e- 13 organization or biogenesis T (G0:0030837) 56 13 8 4.52 05 1 (cytoskeleton)

[0563] 5.9 0.0 cellular component negative regulation of protein 3.6 7e- 22 organization or biogenesis polymerization (G0:0032272) 71 14 5 3.84 05 7 (cytoskeleton) homophilic cell adhesion via 6.1 3.2

[0564] plasma membrane adhesion 8.4 6e- e- molecules (G0:0007156) 164 31 2 3.68 09 05 cell adhesion

[0565] 0.0

[0566] regulation of Rho protein 00 0.0

[0567] signal transduction 4.1 20 47

[0568] (G0:0035023) 81 14 6 3.36 4 5 signal transduction cellular component negative regulation of protein- 1.6 0.0 organization or biogenesis containing complex assembly 6.7 2e- 09 (protein-containing (G0:0031333) 131 21 3 3.12 05 01 complex)

[0569] negative regulation of 7.5 0.0 cellular component supramolecular fiber 7.4 9e- 05 organization or biogenesis organization (GO:1902904) 145 23 5 3.09 06 63 (supramolecular fiber) negative regulation of 4.6 0.0 cellular component cytoskeleton organization 7.2 9e- 19 organization or biogenesis (G0:0051494) 142 21 9 2.88 05 2 (cytoskeleton)

[0570] 0.0

[0571] cell-cell adhesion via plasma- 1.3 00

[0572] membrane adhesion 13. 1e- 25

[0573] molecules (G0:0098742) 257 37 2 2.80 07 5 cell adhesion

[0574] 6.5 0.0

[0575] negative regulation of cell-cell 9.6 2e- 24

[0576] adhesion (G0:0022408) 187 25 1 2.60 05 2 cell adhesion

[0577] 6.6 0.0

[0578] negative regulation of cell 14. 8e- 05

[0579] adhesion (G0:0007162) 277 35 23 2.46 06 47 cell adhesion

[0580] 0.0

[0581] 00 0.0 cellular component regulation of protein 9.7 15 40 organization or biogenesis polymerization (G0:0032271) 190 24 6 2.46 3 3 (cytoskeleton)

[0582] 3.2 2.5

[0583] cell-cell adhesion 26. 1e- e- (G0:0098609) 518 63 61 2.37 09 05 cell adhesion regulation of small GTPase 1.1 0.0

[0584] mediated signal transduction 14. 7e- 07

[0585] (G0:0051056) 290 35 9 2.35 05 59 signal transduction 6.8 0.0 cellular component regulation of actin filament 3e- 24 organization or biogenesis organization (GO:0110053) 253 30 13 2.31 05 7 (cytoskeleton)

[0586] 0.0

[0587] regulation of l-kappaB 00 0.0

[0588] kinase / NF-kappaB signaling 12. 15 41

[0589] (G0:0043122) 240 28 33 2.27 1 2 signal transduction cellular component regulation of protein- 2.2 0.0 organization or biogenesis containing complex assembly 19. 9e- 02 (protein-containing (G0:0043254) 389 45 98 2.25 06 75 complex) regulation of supramolecular 2.3 0.0 cellular component fiber organization 18. 2e- 12 organization or biogenesis (G0:1902903) 351 39 03 2.16 05 4 (supramolecular fiber)

[0590] 4.1 6.5

[0591] 47. 8e- e-

[0592]

[0593] cell adhesion (G0:0007155) 921 99 31 2.09 11 07 cell adhesion0.0

[0594] 00 0.0 cellular component regulation of actin filament- 19. 20 47 organization or biogenesis based process (G0:0032970) 379 38 47 1.95 8 6 (cytoskeleton)

[0595] 4.3 0.0

[0596] positive regulation of kinase 24. 8e- 19 regulation of molecular activity (G0:0033674) 470 47 14 1.95 05 5 function

[0597] 0.0

[0598] 00 0.0

[0599] regulation of cell-cell 23. 14 40

[0600] adhesion (G0:0022407) 465 45 89 1.88 7 8 cell adhesion

[0601] 3.1 0.0

[0602] regulation of cell adhesion 38. 8e- 03

[0603] (G0:0030155) 755 71 78 1.83 06 31 cell adhesion

[0604] 5.4 0.0

[0605] regulation of cell migration 45. 5e- 21

[0606] (G0:0030334) 880 75 2 1.66 05 2 cell motility regulation of cellular 4.3

[0607] component biogenesis 45. 9e- 0.0 cellular component (G0:0044087) 892 76 82 1.66 05 19 organization or biogenesis 0.0

[0608] 00 0.0

[0609] regulation of locomotion 50. 17 42

[0610] (G0:0040012) 981 79 39 1.57 1 9 locomotion

[0611] 0.0

[0612] 00 0.0

[0613] intracellular signal 73. 15 39

[0614] transduction (G0:0035556) 1426 107 25 1.46 5 7 signal transduction 0.0

[0615] negative regulation of 00 0.0

[0616] response to stimulus 78. 12 35

[0617] (G0:0048585) 1533 114 75 1.45 4 7 response to stimulus regulation of intracellular 9.1

[0618] signal transduction 83. 3e- 0.0

[0619] (GO:1902531) 1631 121 78 1.44 05 29 signal transduction 7.3 0.0

[0620] cell surface receptor signaling 98. 7e- 24

[0621] pathway (G0:0007166) 1922 139 73 1.41 05 4 signal transduction 11 3.9 0.0

[0622] regulation of catalytic activity 5.3 1e- 18 regulation of molecular (G0:0050790) 2246 160 8 1.39 05 4 function

[0623] regulation of cellular 11 7.2 0.0

[0624] component organization 4.9 8e- 24 cellular component (G0:0051128) 2238 158 6 1.37 05 6 organization or biogenesis 16 2.5 0.0

[0625] regulation of signaling 4.4 9e- 02

[0626] (G0:0023051) 3201 223 3 1.36 06 88 signal transduction regulation of cell 16 3.9 0.0

[0627] communication 3.8 e- 03

[0628] (G0:0010646) 3190 221 7 1.35 06 57 cell communication 14 2.6 0.0

[0629] regulation of signal 5.0 4e- 12

[0630] transduction (G0:0009966) 2824 195 7 1.34 05 9 signal transduction 0.0

[0631] 12 00 0.0

[0632] regulation of developmental 0.6 15 40

[0633] process (G0:0050793) 2349 162 7 1.34 4 1 developmental process regulation of multicellular 13 9.1 0.0

[0634] organismal process 3.9 6e- 28 regulation of multicellular (G0:0051239) 2608 179 7 1.34 05 5 organismal process 12 0.0 0.0

[0635] cellular response to chemical 7.8 00 44

[0636] stimulus (G0:0070887) 2489 170 6 1.33 18 5 response to stimulus 0.0

[0637] 22 1.0 00

[0638] negative regulation of cellular 6.0 1e- 26 regulation of molecular process (G0:0048523) 4401 300 8 1.33 07 3 function

[0639] 19 6.0 0.0

[0640] regulation of response to 0.9 3e- 05

[0641] stimulus (G0:0048583) 3717 250 4 1.31 06 22 response to stimulus 0.0

[0642] 14 00 0.0

[0643] regulation of molecular 9.3 18 44 regulation of molecular

[0644]

[0645] function (G0:0065009) 2907 194 3 1.30 2 2 function24 2.4 0.0

[0646] anatomical structure 8.6 4e- 00

[0647]

[0648] development (G0:0048856) 4840 322 3 1.30 07 38 developmental process

[0649] Coun

[0650] genesEPIC - t Exp Fold P- REFLIST (4402 ecte Enrich val FD Reduced GO term: biological process (17968) ) d ment ue R Term no homophilic cell adhesion via plasma 6.3 9.9

[0651] n- membrane adhesion molecules 7e- 2e- cell m (G0:0007156) 164 37 6.86 5.39 15 11 adhesion FT cell-cell adhesion via plasma- 1.4 7.4

[0652] C membrane adhesion molecules 10.7 4e- 5e- cell vs (G0:0098742) 257 42 6 3.9 12 09 adhesion N 2.3 1.8

[0653] T 21.6 7e- 5e- cell cell-cell adhesion (G0:0098609) 518 64 8 2.95 13 09 adhesion 2.3 9.1

[0654] 38.5 5e- 6e- cell cell adhesion (G0:0007155) 921 89 5 2.31 12 09 adhesion 1.7 0.0 developm nervous system development 86.5 e- 004 ental (G0:0007399) 2069 137 9 1.58 07 42 process 2.8 developm regulation of developmental process 98.3 9e- 0.0 ental (G0:0050793) 2349 140 1 1.42 05 451 process 1.0 0.0 developm 151. 3e- 003 ental system development (G0:0048731) 3622 214 59 1.41 07 2 process 7.0 0.0 developm multicellular organism development 167. 4e- 015 ental (G0:0007275) 3991 227 03 1.36 07 7 process 7.5 0.0 developm anatomical structure development 202. 9e- 014 ental

[0655]

[0656] (G0:0048856) 4840 266 56 1.31 07 8 process

[0657] genesEPI Cou Ex Fold

[0658] C - nt pe Enric P- GO term: biological REFLIST (113 cte hmen val FD

[0659] process (17968) 2) d t ue R Reduced Term mP 2.

[0660] TO homophilic cell adhesion via 1.6 59

[0661] vs plasma membrane adhesion 10. 7e- e- NT molecules (G0:0007156) 164 40 33 3.87 11 07 cell adhesion

[0662] 1.

[0663] cell-cell adhesion via 2.2 76

[0664] plasma-membrane adhesion 16. 6e- e- molecules (G0:0098742) 257 51 19 3.15 11 07 cell adhesion

[0665] 0.0 cellular component negative regulation of 00 0. organization or biogenesis protein-containing complex 8.2 10 04 (protein-containing assembly (GO:0031333) 131 22 5 2.67 1 48 complex)

[0666] 4.

[0667] 1.4 41

[0668] cell-cell adhesion 32. 2e- e- (G0:0098609) 518 74 63 2.27 09 06 cell adhesion

[0669] 3.

[0670] 1.8 56

[0671] 58. 3e- e- cell adhesion (G0:0007155) 921 106 02 1.83 08 05 cell adhesion

[0672] 1.

[0673] 13 4.7 23

[0674] nervous system 0.3 4e- e- development (G0:0007399) 2069 200 5 1.53 09 05 developmental process positive regulation of 5.5 0.

[0675] multicellular organismal 90. 6e- 02 regulation of multicellular

[0676]

[0677] process (GO:0051240) 1439 131 66 1.44 05 79 organismal process3.

[0678] 22 9.6 74

[0679] system development 8.1 e- e- (G0:0048731) 3622 317 9 1.39 10 06 developmental process 1.

[0680] 25 5.3 19

[0681] multicellular organism 1.4 3e- e- development (G0:0007275) 3991 339 4 1.35 09 05 developmental process anatomical structure 13 7.8 0.

[0682] morphogenesis 4.6 2e- 03

[0683]

[0684] (G0:0009653) 2137 181 3 1.34 05 69 developmental process

[0685] Coun

[0686] genesEPIC - t Exp Fold P- REFLIST (4402 ecte Enrich val FD Reduced

[0687]

[0688] GO term: biological process (17968) ) d ment ue R Term m homophilic cell adhesion via plasma 5.1 5.6

[0689] FT membrane adhesion molecules e- 8e- cell C (G0:0007156) 164 62 24.3 2.55 09 06 adhesion vs cell-cell adhesion via plasma- 1.2 3.7 response N membrane adhesion molecules 38.0 e- 4e- to

[0690] T (G0:0098742) 257 89 8 2.34 10 07 stimulus detection of chemical stimulus 2.8 2.2 response involved in sensory perception of 51.5 9e- 5e- to smell (G0:0050911) 348 111 6 2.15 11 07 stimulus 4.8 2.5 response sensory perception of smell 55.4 5e- 2e- (G0:0007608) 374 116 1 2.09 11 07 stimulus detection of chemical stimulus 1.1 2.1 response involved in sensory perception 58.3 e- 3e- (G0:0050907) 394 115 7 1.97 09 06 stimulus 3.2 8.3 response detection of stimulus involved in 3e- 9e- to sensory perception (G0:0050906) 461 131 68.3 1.92 10 07 stimulus 1.4 2.3 response sensory perception of chemical 66.6 7e- e- to stimulus (G0:0007606) 450 126 7 1.89 09 06

[0691] 1.4 1.2 response detection of chemical stimulus 63.8 8e- 8e- (G0:0009593) 431 118 5 1.85 08 05 stimulus 1.3 8.0

[0692] 76.7 5e- 9e- cell cell-cell adhesion (G0:0098609) 518 131 4 1.71 07 05 adhesion 2.5 1.9 response 85.6 7e- 1e- to detection of stimulus (GO:0051606) 578 146 3 1.7 08 05 stimulus 2.4 1.9 response 129. 4e- e- sensory perception (G0:0007600) 876 202 78 1.56 08 05 stimulus 2.7 0.0

[0693] 136. 4e- 001 cell cell adhesion (G0:0007155) 921 204 45 1.5 07 38 adhesion 1.5 8.9 signal G protein-coupled receptor signaling 4e- e- transduct pathway (G0:0007186) 1080 234 160 1.46 07 05 ion 1.6 1.3

[0694] nervous system process 195. 9e- 8e- system (G0:0050877) 1321 283 71 1.45 08 05 process 8.0 7.8

[0695] 281. 3e- 1e- system

[0696]

[0697] system process (G0:0003008) 1901 385 64 1.37 09 06 process

[0698] Coun genesEPIC - t Exp Fold P- REFLIST (4402 ecte Enrich val FD Reduced GO term: biological process (17968) ) d ment ue R Term p detection of chemical stimulus 2.5 4.0 response T involved in sensory perception of 85.2 8e- 1e- ■■■

[0699]

[0700] C- smell (G0:0050911) 348 223 6 2.62 27 23 stimulusD 2.0 1.5 response M sensory perception of smell 91.6 1e- 6e- to V (G0:0007608) 374 227 3 2.48 25 21 stimulus s detection of chemical stimulus 8.5 3.3 response N involved in sensory perception 96.5 1e- 1e- T (G0:0050907) 394 229 3 2.37 24 20 stimulus 2.0 1.0 response detection of stimulus involved in 112. 9e- 9e- sensory perception (G0:0050906) 461 256 94 2.27 24 20 stimulus 8.5 2.6 response detection of chemical stimulus 105. 4e- 6e- to (G0:0009593) 431 239 59 2.26 23 19 stimulus 1.0 2.7 response sensory perception of chemical 110. 7e- 8e- to stimulus (G0:0007606) 450 242 25 2.2 21 18

[0701] 1.4 3.2 response 141. 5e- 2e- to detection of stimulus (GO:0051606) 578 287 6 2.03 21 18 stimulus homophilic cell adhesion via plasma 6.3 0.0 membrane adhesion molecules 40.1 8e- 001 cell (G0:0007156) 164 81 8 2.02 07 99 adhesion cell-cell adhesion via plasma- 2.2 9.9 membrane adhesion molecules 62.9 4e- 8e- cell (G0:0098742) 257 124 6 1.97 09 07 adhesion 0.0 developm embryonic skeletal system 31.1 001 0.0 ental development (G0:0048706) 127 58 1 1.86 77 275 process 9.1 1.4 response 214. 7e- 3e- sensory perception (G0:0007600) 876 378 61 1.76 20 16 stimulus 7.9 1.5 signal G protein-coupled receptor signaling 264. 3e- 4e- transducti pathway (G0:0007186) 1080 448 59 1.69 21 17 on 3.7 1.2

[0702] 126. 1e- 8e- cell cell-cell adhesion (G0:0098609) 518 202 91 1.59 08 05 adhesion 5.7 6.4 nervous system process 323. 8e- 3e- system (G0:0050877) 1321 500 63 1.54 17 14 process 4.1 0.0 synaptic trans-synaptic signaling 100. 2e- 087 transmissi (G0:0099537) 411 150 69 1.49 05 9 on 9.2 synaptic chemical synaptic transmission 96.7 2e- 0.0 transmissi (G0:0007268) 395 143 7 1.48 05 173 on 9.2 synaptic anterograde trans-synaptic signaling 96.7 2e- 0.0 transmissi (G0:0098916) 395 143 7 1.48 05 171 on 1.3 1.6

[0703] 465. 8e- 6e- system system process (G0:0003008) 1901 674 73 1.45 17 14 process 0.0 synaptic 107. 001 0.0 transmissi synaptic signaling (G0:0099536) 438 153 31 1.43 91 294 on 9.2 3.0

[0704] 225. 3e- 6e- cell cell adhesion (G0:0007155) 921 319 64 1.41 08 05 adhesion 0.0 developm skeletal system development 120. 001 0.0 ental (G0:0001501) 491 169 29 1.4 55 252 process 0.0 developm 121. 001 0.0 ental chemotaxis (G0:0006935) 494 170 03 1.4 6 255 process 0.0 developm 126. 001 0.0 ental locomotion (G0:0040011) 515 177 17 1.4 21 204 process 0.0 developm 121. 001 0.0 ental

[0705]

[0706] taxis (G0:0042330) 496 170 52 1.4 99 301 processTable 10. Gene ontology (GO) enrichment analysis of the genes associated with CpGs that are differentially methylated at their 5' region. GO term: molecular function.

[0707]

[0708] genesEPIC - Exp Fold P- GO term: molecular REFLIST Count ecte Enrich val FD Reduced funcion (17968) (923) d ment ue R Term

[0709] non- 3.36 0.00 cell adhesion mPT cadherin binding 16.1 e- 054 molecule C vs (G0:0045296) 315 41 8 2.53 07 9 binding NT 6.17

[0710] calcium ion binding e- 3.02 calcium ion (G0:0005509) 697 76 35.8 2.12 09 e-05 binding 2.65

[0711] phospholipid binding 22.8 e- 0.02 (G0:0005543) 445 46 6 2.01 05 16 lipid binding 5.00 cell adhesion cell adhesion molecule 27.1 E- 0.03 molecule binding (G0:0050839) 529 51 7 1.88 05 5 binding 1.03

[0712] 40.7 e- 0.00

[0713]

[0714] lipid binding (G0:0008289) 794 76 9 1.86 06 126 lipid binding genesEPIC - Exp Fold P- GO term: molecular REFLIST Count ecte Enrich val FD Reduced funcion (17968) (1132) d ment ue R Term mPT 6.49 0.00

[0715] C vs calcium ion binding 43.9 e- 031 calcium ion

[0716]

[0717] NT (G0:0005509) 697 85 1 1.94 08 8: —

[0718] genesEPIC - Exp Fold P- GO term: molecular REFLIST Count ecte Enrich val FD Reduced funcion (17968) (4402) d ment ue R Term

[0719] PTC- 1.82

[0720] DM odorant binding 23.0 e- 1.28 odorant (G0:0005549) 94 67 3 2.91 10 e-07 binding 1.95 signaling olfactory receptor activity 84.7 e- 9.54 receptor (G0:0004984) 346 223 7 2.63 27 e-24 activity 2.37 signaling G protein-coupled receptor 184. e- 5.81 receptor activity (G0:0004930) 752 354 23 1.92 23 e-20 activity transmembrane signaling 1.59 signaling receptor activity 286. e- 1.56 receptor (G0:0004888) 1170 485 64 1.69 22 e-19 activity 3.98 signaling molecular transducer 337. e- 6.51 receptor activity (G0:0060089) 1379 553 84 1.64 23 e-20 activity 3.98 signaling signaling receptor activity 337. e- 4.88 receptor (G0:0038023) 1379 553 84 1.64 23 e-20 activity 1.12 0.00

[0721] calcium ion binding 170. e- 068 calcium ion

[0722]

[0723] (G0:0005509) 697 246 76 1.44 06 9 binding

[0724]

[0725] Cou

[0726] genesEPIC nt Ex Fold P- -REFLIST (440 pec Enrich val FD Reduced

[0727]

[0728] GO term: molecular funcion (17968) 2) ted ment ue R Termnon- 2.2 1.1 ■■■■■I mFTC 29. 6e- 1e- MHHHB

[0729]

[0730] vs NT calcium ion binding (G0:0005509) 697 67 17 2.3 09 05 binding

[0731] Cou

[0732] genesEPIC nt Ex Fold P- -REFLIST (440 pec Enrich val FD Reduced GO term: molecular funcion (17968) 2) ted ment ue R Term mFTC 1.8 2.2

[0733] vs NT 13. 4e- 5e- odorant odorant binding (G0:0005549) 94 43 93 3.09 08 05 binding 2.5 1.2 signaling olfactory receptor activity 51. 3e- 4e- receptor (G0:0004984) 346 111 26 2.17 11 07 activity 6.9 1.6 signaling G protein-coupled receptor activity 111 1e- 9e- receptor (G0:0004930) 752 187.41 1.68 10 06 activity 8.9 1.4 signaling transmembrane signaling receptor 173 7e- 7e- receptor activity (G0:0004888) 1170 258.34 1.49 09 05 activity 4.7 4.6 signaling molecular transducer activity 204 3e- 4e- receptor (G0:0060089) 1379 290.3 1.42 08 05 activity 4.7 3.8 signaling signaling receptor activity 204 3e- 6e- receptor

[0734]

[0735] (G0:0038023) 1379 290.3 1.42 08 05 activity Cou

[0736] genesEPIC nt Ex Fold P- - REFLIST (440 pec Enrich val FD Reduced GO term: molecular funcion (17968) 2) ted ment ue R Term FTC- 0.0 signaling DM vs neurotransmitter receptor activity 17. 001 0.0 receptor NT (G0:0030594) 104 37 3 2.14 55 421 activity 9.4 00

[0737] DNA-binding transcription activator 75. 9e- 003 DNA activity (G0:0001216) 455 126 69 1.66 07 88 binding 0.0 ansmem 54. 001 00 brane gated channel activity (G0:0022836) 327 87 4 1.6 58 408 transport 0.0 ansmem 71. 001 00 brane ion channel activity (G0:0005216) 428 108 2 1.52 64 383 transport 1.8 00 calcium 115 7e- 006 ion calcium ion binding (G0:0005509) 697 175.95 1.51 06 54 binding RNA polymerase II cis-regulatory 1.6 80 region sequence-specific DNA binding 189 4e- 2e- DNA (G0:0000978) 1139 283.47 1.49 09 06 binding 0.0 ansmem 78. 001 00 brane channel activity (GO:0015267) 470 116 19 1.48 95 434 transport 2.2 54 cis-regulatory region sequence- 192 e- e- DNA specific DNA binding (G0:0000987) 1159 286.8 1.48 09 06 binding 0.0 ansmem passive transmembrane transporter 78. 001 00 brane activity (G0:0022803) 471 116 35 1.48 99 425 transport RNA polymerase II transcription 3.5 43 regulatory region sequence-specific 223 8e- 9e- DNA DNA binding (G0:0000977) 1341 321.08 1.44 09 06 binding 1.4 1 4

[0738] DNA-binding transcription factor 224 5e- 3e- DNA activity (G0:0003700) 1352 319.91 1.42 08 05 binding 3.4 28 transcription regulatory region nucleic 239 5e- 2e- DNA acid binding (G0:0001067) 1438 333.21 1.39 08 05 binding 4.1 29 transcription cis-regulatory region 238 4e- e- DNA

[0739]

[0740] binding (G0:0000976) 1436 332.88 1.39 08 05 binding4.6 2.5 sequence-specific double-stranded 247 1e- 1e- DNA DNA binding (GO: 1990837) 1490 342.86 1.38 08 05 binding 4.2 2.6 sequence-specific DNA binding 264 5e- e- DNA (G0:0043565) 1592 362.83 1.37 08 05 binding 4.0 0.0

[0741] double-stranded DNA binding 2e- 001 DNA (G0:0003690) 1575 351 262 1.34 07 79 binding 2.2 0.0 transcription regulator activity 303 9e- 001 DNA

[0742]

[0743] (G0:0140110) 1822 400.09 1.32 07 12 binding

[0744] Table 11. Gene ontology (GO) enrichment analysis of the genes associated with CpGs that are differentially methylated at their 5' region. GO term: cellular component.

[0745]

[0746] PTC genesEPIC - Coun Exp Fold P- REFLIST t ecte Enrich val FD Reduced GO term: cellular component (17968) (923) d ment ue R Tenn no o.o BSHH n- 001 0.0 cvtoskeleto m podosome (G0:0002102) 31 9 1.59 5.65 06 238 n PT 1.2 C 20.9 8e- 0.0 cell vs focal adhesion (G0:0005925) 408 44 6 2.1 05 043 junction NT 2.2 0.0 cell-substrate junction 21.3 6e- 065 cell (G0:0030055) 415 44 2 2.06 05 3 junction 0.0 membrane microdomain 16.2 003 0.0 (G0:0098857) 317 33 8 2.03 63 386 membrane 0.0 16.2 003 0.0 membrane raft (G0:0045121) 317 33 8 2.03 63 367 membrane 8.3 5.6 65.0 9e- 5e- cell anchoring junction (G0:0070161) 1266 116 3 1.78 09 06 junction 3.4 1.7 103. 6e- 5e- cell cell junction (G0:0030054) 2022 162 87 1.56 08 05 junction 0.0 plasma membrane region 61.2 001 0.0 (G0:0098590) 1193 93 8 1.52 26 255 membrane 0.0 67.6 001 0.0 cell neuron projection (G0:0043005) 1317 100 5 1.48 9 214 projection 0.0 integral component of plasma 82.3 001 0.0 membrane (G0:0005887) 1603 118 4 1.43 65 196 membrane 0.0 intrinsic component of plasma 001 0.0 membrane (GO:0031226) 1680 123 86.3 1.43 42 24 membrane 2.6 0.0 116. 6e- 067 cell cell projection (G0:0042995) 2263 162 25 1.39 05 2 projection 0.0 extracellular exosome 100. 001 0.0 extracellula (G0:0070062) 1962 140 79 1.39 32 243 r region 0.0 101. 001 0.0 extracellula extracellular vesicle (GO:1903561) 1983 141 86 1.38 42 221 r region 0.0 extracellular organelle 101. 001 0.0 extracellula

[0747]

[0748] (G0:0043230) 1984 141 92 1.38 43 206 r region0.0 extracellular membrane-bounded 101. 001 0.0 extracellula organelle (G0:0065010) 1984 141 92 1.38 43 193 r region 8.5 1.7

[0749] 291. 2e- 26- cell periphery (GO:0071944) 5671 398 31 1.37 13 09 periphery 0.0

[0750] plasma membrane bounded cell 110. 001 0.0 cell projection (GO:0120025) 2156 151 75 1.36 61 204 projection 5.5 5.5

[0751] 267 e- 5e-

[0752]

[0753] plasma membrane (G0:0005886) 5208 358 53 1.34 10 07 membrane Coun

[0754] genesEPIC - t Exp Fold P- REFLIST (1132 ecte Enrich val FD Reduced GO term: cellular component (17968) ) d ment ue R Tenn m 0.0

[0755] PT integral component of synaptic 002 0.0

[0756] C membrane (G0:0099699) 150 23 9.45 2.43 59 402 membrane vs 0.0

[0757] NT intrinsic component of synaptic 10.2 004 0.0 membrane (G0:0099240) 162 24 1 2.35 14 465 membrane 0.0

[0758] postsynaptic density 19.4 003 0.0 cell (G0:0014069) 309 38 7 1.95 08 415 junction 0.0

[0759] postsynaptic specialization 20.6 002 0.0 cell (G0:0099572) 328 40 6 1.94 01 339 junction 0.0

[0760] asymmetric synapse 19.7 003 0.0 cell (G0:0032279) 314 38 8 1.92 51 444 junction 0.0

[0761] glutamatergic synapse 19.8 003 0.0 cell (G0:0098978) 315 38 5 1.91 63 431 junction 5.0

[0762] 37.4 0E- 0.0 cell postsynapse (G0:0098794) 594 65 2 1.74 05 101 junction 5.5 3.7 intrinsic component of plasma 105. 5e- 4e- membrane (G0:0031226) 1680 174 84 1.64 10 07 membrane 2.5 1.2 integral component of plasma 100. 1e- 7e- membrane (G0:0005887) 1603 165 99 1.63 09 06 membrane 4.4

[0763] plasma membrane region 75.1 e- 0.0 (G0:0098590) 1193 113 6 1.5 05 111 membrane 7.3 7.4

[0764] 328. 2e- e- plasma membrane (G0:0005886) 5208 450 11 1.37 14 11 membrane 3.6 7.3

[0765] 357. 5e- 8e- cell periphery (GO:0071944) 5671 483 28 1.35 14 11 periphery 0.0

[0766] 127. 001 0.0 cell

[0767]

[0768] cell junction (G0:0030054) 2022 171 39 1.34 42 261 junction Coun

[0769] genesEPIC - t Exp Fold P- REFLIST (4402 ecte Enrich val FD Reduced GO term: cellular component (17968) ) d ment ue R Term PT 0.0

[0770] C- integral component of postsynaptic 12.4 009 0.0

[0771] D density membrane (G0:0099061 ) 51 28 9 2.24 59 318 membrane M 0.0

[0772] intrinsic component of postsynaptic 13.4 008 0.0 density membrane (G0:0099146) 55 30 7 2.23 45 294 membrane 1.2 0.0 integral component of postsynaptic 27.9 3e- 007 membrane (G0:0099055) 114 58 3 2.08 05 99 membrane 8.0 0.0 intrinsic component of postsynaptic 5e- 005

[0773]

[0774] membrane (G0:0098936) 120 61 29.4 2.07 06 81 membra neintegral component of postsynaptic 0.0 specialization membrane 17.6 009 0.0 (G0:0099060) 72 36 4 2.04 03 304 membrane intrinsic component of postsynaptic 0.0

[0775] specialization membrane 18.6 005 0.0 (G0:0098948) 76 38 2 2.04 41 215 membrane 2.5 0.0

[0776] intrinsic component of synaptic 39.6 4e- 002 membrane (G0:0099240) 162 78 9 1.97 06 14 membrane 1.1 0.0

[0777] integral component of synaptic 36.7 7e- 007 membrane (G0:0099699) 150 71 5 1.93 05 9 membrane 0.0

[0778] intrinsic component of presynaptic 19.8 016 0.0 membrane (G0:0098889) 81 38 4 1.91 1 485 membrane 0.0

[0779] postsynaptic density membrane 21.5 011 0.0 (G0:0098839) 88 41 6 1.9 4 366 membrane 0.0 0.0 presynaptic membrane 35.2 001 066 (G0:0042734) 144 64 8 1.81 45 4 membrane 2.8 0.0

[0780] 3e- 015

[0781] synaptic membrane (G0:0097060) 360 136 88.2 1.54 05 9 membrane 0.0

[0782] postsynaptic membrane 62.4 006 0.0 (G0:0045211) 255 95 7 1.52 73 247 membrane 1.0 6.8

[0783] intrinsic component of plasma 411. 1e- 1e- membrane (G0:0031226) 1680 611 59 1.48 17 15 membrane 1.3 6.6

[0784] integral component of plasma 392. 1e- 1e- membrane (G0:0005887) 1603 581 72 1.48 16 14 membrane 0.0

[0785] neuron to neuron synapse 82.8 010 0.0 (G0:0098984) 338 118 1 1.43 5 342 membrane 0.0

[0786] 150. 003 cell axon (G0:0030424) 616 202 91 1.34 33 14 projection 0.0

[0787] 146. 006 0.0 cell dendritic tree (G0:0097447) 597 194 26 1.33 32 241 projection 0.0

[0788] 145. 007 0.0 cell dendrite (GO: 0030425) 595 193 77 1.32 26 262 projection 1.2 1.2

[0789] 127 7e- 8e- plasma membrane (G0:0005886) 5208 1682 5.91 1.32 31 28 membrane 9.6 1.9

[0790] 138 6e- 5e- cell periphery (GO:0071944) 5671 1818 9.34 1.31 34 30 periphery 0.0 cellular somatodendritic compartment 199. 002 0.0 anatomical

[0791]

[0792] (G0:0036477) 815 259 67 1.3 5 107 entity

[0793]

[0794] genesEPIC - Exp Fold P- GO term: cellular REFLIST Count ecte Enrich valu Reduced

[0795]

[0796] component (17968) (4402) d ment e FDR Teim no intrinsic component of plasma 70.3 3.3e 6.68

[0797] n- membrane (GO:0031226) 1680 123 1 1.75 -09 e-06 membrane m integral component of plasma 67.0 1.07 7.22

[0798] FT membrane (G0:0005887) 1603 117 9 1.74 e-08 e-06 membrane C plasma membrane 217. 6.2e 3.13

[0799] vs (G0:0005886) 5208 289 97 1.33 -08 e-05 membrane NT 237. 1.02 1.03

[0800]

[0801] cell periphery (GO:0071944) 5671 314 34 1.32 e-08 e-05 peripherygenesEPIC - Exp Fold P- GO term: cellular REFLIST Count ecte Enrich valu Reduced component (17968) (4402) d ment e FDR Term m 0.00

[0802] FT 39.1 045 0.03 cell C distal axon (GO:0150034) 264 65 1 1.66 7 85 projection vs intrinsic component of plasma 248. 1.36 9.15

[0803] NT membrane (GO:0031226) 1680 359 9 1.44 e-10 e-08 membrane integral component of plasma 237. 1.07 5.43 membrane (G0:0005887) 1603 340 49 1.43 e-09 e-07 membrane 0.00 cellular somatodendritic compartment 120. 035 0.03 anatomical (G0:0036477) 815 164 74 1.36 6 13 entity plasma membrane region 176. 9.83

[0804] (G0:0098590) 1193 233 75 1.32 e-05 11 membrane 0.00

[0805] anchoring junction 187. 014 0.01

[0806]

[0807] (G0:0070161) 1266 244 56 1.3 8 58 cell junction

[0808] genesEPIC - Exp Fold P- REFLIST Count ecte Enrich valu Reduced GO term: biological process (17968) (4402) d ment e FDR Term FT 0.00

[0809] C- chromaffin granule 098 0.03

[0810] D (G0:0042583) 14 10 2.33 4.29 8 57

[0811] M ionotropic glutamate receptor 5.04 0.00

[0812] vs complex (G0:0008328) 36 20 5.99 3.34 e-05 51 membrane NT 0.00 neurotransmitter receptor 027 0.01 complex (G0:0098878) 41 20 6.82 2.93 9 66 membrane intrinsic component of 0.00

[0813] postsynaptic membrane 19.9 036 0.01 (G0:0098936) 120 40 6 2 3 83 membrane 28.1 4.44 0.00 cell dendritic spine (G0:0043197) 169 55 1 1.96 e-05 527 projection integral component of 0.00

[0814] postsynaptic membrane 18.9 092 0.03 (G0:0099055) 114 37 6 1.95 8 41 membrane 28.2 4.68 0.00 cell neuron spine (G0:0044309) 170 55 8 1.94 e-05 498 projection extrinsic component of 0.00

[0815] plasma membrane 27.7 028 0.01 (G0:0019897) 167 51 8 1.84 1 62 periphery intrinsic component of 0.00

[0816] synaptic membrane 26.9 047 0.02 (G0:0099240) 162 49 5 1.82 5 13 membrane integral component of 0.00

[0817] synaptic membrane 24.9 090 0.03 (G0:0099699) 150 45 5 1.8 4 38 membrane 0.00

[0818] cation channel complex 34.9 015 0.01 (G0:0034703) 210 62 3 1.77 4 07 membrane 0.00

[0819] postsynaptic membrane 42.4 040 0.01 (G0:0045211) 255 70 2 1.65 1 88 membrane 0.00

[0820] ion channel complex 46.7 017

[0821] (G0:0034702) 281 77 4 1.65 4 11 membrane 189. 2.44 1.64 chromatin (G0:0000785) 1141 305 81 1.61 e-13 e-10 chromosome 0.00

[0822] neuron to neuron synapse 56.2 011 0.00 (G0:0098984) 338 90 3 1.6 5 933 cell junction 0.00

[0823] postsynaptic specialization 54.5 016 0.01 (G0:0099572) 328 87 6 1.59 4 07 cell junction 0.00

[0824] 43.9 086 cell

[0825]

[0826] distal axon (GO:0150034) 264 70 2 1.59 6 33 projection0.00 postsynaptic density 042 0.01 (G0:0014069) 309 81 51.4 1.58 4 95 cell junction 0.00

[0827] glutamatergic synapse 047 0.02 (G0:0098978) 315 82 52.4 1.56 6 09 cell junction 0.00

[0828] asymmetric synapse 52.2 061 0.02 (G0:0032279) 314 81 3 1.55 2 58 cell junction 0.00

[0829] extrinsic component of 083 0.03 membrane (GO:0019898) 303 78 50.4 1.55 8 39 periphery 0.00

[0830] synaptic membrane 59.8 049 0.02 (G0:0097060) 360 91 9 1.52 7 14 membrane 102. 5.99 0.00 cell axon (G0:0030424) 616 155 47 1.51 e-06 101 projection 98.8 2.26 0.00 postsynapse (G0:0098794) 594 147 1 1.49 e-05 304 cell junction 98.9 3.61 0.00 cell dendrite (G0:0030425) 595 146 8 1.48 e-05 455 projection plasma membrane protein 93.8 7.3e 0.00

[0831] complex (G0:0098797) 564 138 2 1.47 -05 671 membrane 99.3 4.63 0.00 cell dendritic tree (G0:0097447) 597 146 1 1.47 e-05 52 projection cell-cell junction 79.6 0.00 0.01 (G0:0005911) 479 116 8 1.46 035 86 cell junction 0.00 cellular somatodendritic compartment 135. 3.07 056 anatomical (G0:0036477) 815 197 58 1.45 e-06 4 entity 0.00 protein- transcription regulator 80.0 084 0.03 containing complex (G0:0005667) 481 114 1 1.42 2 27 complex integral component of plasma 266. 5.79 1.95 membrane (G0:0005887) 1603 378 66 1.42 e-10 e-07 membrane intrinsic component of plasma 279. 1.94 7.85 membrane (GO:0031226) 1680 396 47 1.42 e-10 e-08 membrane anchoring junction 210. 1.07 3.08 (G0:0070161) 1266 296 6 1.41 e-07 e-05 cell junction actin cytoskeleton 81.1 0.00 0.04 (G0:0015629) 488 114 8 1.4 137 53 cytoskeleton cellular 90.9 0.00 0.03 anatomical cell body (G0:0044297) 547 126 9 1.38 113 95 entity 213. 4.05 9.1e

[0832] synapse (G0:0045202) 1282 295 26 1.38 e-07 -05 cell junction neuron projection 219. 2.6e 6.56 cell (G0:0043005) 1317 303 08 1.38 -07 e-05 projection 336. 1.04 5.25

[0833] cell junction (G0:0030054) 2022 464 36 1.38 e-10 e-08 cell junction plasma membrane region 198. 7.85 0.00 (G0:0098590) 1193 268 46 1.35 e-06 122 membrane 0.00

[0834] 278. 5.65 011

[0835]

[0836] chromosome (G0:0005694) 1675 369 64 1.32 e-07 4 chromosome

[0837] Accordingly, GO molecular function (MF) terms involving calcium dependent processes and cellular component (CC) terms related to membrane and cell junction were also significantly enriched. These terms mostly involve cadherin and protocadherin genes. In addition, certain GO terms were associated with specific tumor groups. For example, BP terms associated with the organization or biogenesis of cellular components (mainly thecytoskeleton), and CC terms related to the extracellular region (mainly to exosomes), were highly enriched in non-mPTC.

[0838] We also performed an integrative analysis of our DNA methylation data in non-mPTC and mPTC with TCGA gene expression data on low- and high-risk PTC (Figure 5). Compared to normal tissue, a total of 651 and 944 differentially methylated CpGs in our non-mPTC and mPTC, respectively, were associated with 342 and 567 differentially expressed genes in TCGA low-risk and high-risk PTC, respectively (Table 12).

[0839] Table 12. Differentially methylated CpGs at 5' region in our discovery series associated with differentially expressed genes in The Cancer Genome Atlas (TCGA).

[0840] Nr unique differentially methylated CpGs at 5' 299 382

[0841] Nr unique associated genes 164 258

[0842] non-mPTC Low-risk PTC mPTC vs High-risk PTC Comparison vs NT vs NT NT vs NT Series Discovery TCGA Discovery TCGA eg id gene cgid gene cg26267310 PLEKHA4 cg10234511 CASZ1 cg26992963 SPSB1 cg09727611 AKAP2 cg21236655 TNC cg08994526 AKAP2 cg06001396 LGALS3 cg06001396 LGALS3 cg06182018 LIFR cg09749435 C1orf64 cg16254374 BNC1 cg05213092 AKAP2 cg16594139 PLEKHA4 cg24059418 SHANK2 cg10275315 BNC1 cg16790847 ZIC4 cg07686035 SLC26A4 cg10127822 AKAP2 cg03957481 KLHDC7B cg11523424 SHANK2 cg09749435 C1orf64 cg25923918 CSGALNACT1 cg23876631 PNPLA5 cg05196406 ALDH4A1 cg03392126 DLG2 cg13398586 C1orf64 cg07925867 CDH16 cg10236239 SULT1C4 cg00962755 PHLDA3 cg07667018 TPO cg18244358 BAIAP3 cg19239270 MPPED2 cg27373426 C1orf64 cg01668287 SLC25A42 cg04777726 PLEKHA4 cg16510654 C1orf64 cg13107973 KCNK2 cg01141459 HORMAD2 cg27169020 BNC1 cg15923139 SORBS2 cg16000321 CRYAB cg18279094 FOXD3 cg15868302 FOXD2 cg10714061 IP6K3 cg25923918 CSGALNACT1 cg17179314 CHRNB4 cg14669515 SPSB1 cg09065714 UCN

[0843]

[0844] cg15988320 PCDHGA11 cg27373426 C1orf64cg15988320 PCDHGC5 cg21970261 TFF3 cg15988320 PCDHGC4 cg12732899 AKAP2 cg26341879 PLA2G12B cg24813021 TNS3 cg02526950 MET cg10514711 SHANK2 cg08660634 CDH6 cg16433265 SORBS2 cg13500471 TMC6 cg15474728 ABR cg01724702 SIRPG cg23189410 ZIC4 cg00322915 MGAT3 cg23189410 ZIC1 cg05064785 SLC34A2 cg00962755 PHLDA3 cg20157095 ALDH3B1 cg20591728 HOXB6 cg25677394 HRH1 cg13095301 GRAP2 cg14254452 BBC3 cg23876631 PNPLA5 cg11941811 TMEM108 cg27169020 BNC1 cg14575854 GRB7 cg26867637 ABR cg22217144 CLU cg24663003 SRCRB4D cg17924081 ETV4 cg05528187 AKAP2 cg27397879 FAM20A cg19908768 SULT1C4 cg07003993 SYT12 cg05666820 FZD9 cg17276289 CLDN16 cg00043564 EMID2 cg00852595 ANKS1B cg23136645 TPO cg08995609 RIN1 cg03249747 DNAH8 cg13953433 SLMO1 cg07149030 ABR cg16854606 DAND5 cg25290307 EMID2 cg27092248 HRH1 cg23101261 ACSL6 cg18911683 CRHR2 cg22813391 CDH4 cg01246266 TMC6 cg07925867 CDH16 cg01891789 SEMA3F cg02495208 GRAP2 cg10402521 STARD13 cg22169053 SPOCD1 cg16983159 TMEM173 cg18115040 HOXD10 cg04607238 TMPRSS4 cg18437144 CASZ1 cg00182273 CA7 cg02942546 L3MBTL4 cg08378932 TM4SF4 cg20170028 SHANK2 cg05673137 UNC5CL cg13594711 SRCRB4D cg25655096 LPAR5 cg02829152 IL5RA cg24562785 SIRPG cg00796973 RAP1GAP cg06082750 MGAT3 cg03073060 TNS3 cg14275207 STRA6 cg05036173 SULT1C4 cg15607221 NRIP1 cg15868302 FOXD2 cg18186672 DYSF cg08887581 C1orf64 cg19445684 FGF1 cg15433210 BTBD11 cg15236354 GREM2 cg09583199 CDHR3 cg03762081 KLK10 cg13833437 UCN cg00457165 CDSN cg13026754 C14orf86 cg04253697 BEAN cg23449696 ZIC1 cg22111043 MYO1G cg02973171 FAM181B cg26091557 GRHL3 cg14509403 HORMAD2 cg08686879 CSF2 cg27409154 IGFBP6 cg23634952 ENG cg02531794 SRCRB4Dcg11984463 TTC9B cg19476303 ABR cg06196379 TREM1 cg17068383 MASP1 cg26774156 STRA6 cg04689058 CDH16 cg12323347 DUSP5 cg07964848 MKX cg06321045 CDSN cg03957481 KLHDC7B cg12653876 LGALS3 cg24690071 ZIC2 cg12589564 PRDM1 cg11981672 LINGO1 cg04472592 KRT80 cg22129545 UCN cg20951054 C14orf176 cg06654109 SLC18A3 cg09432154 GPR87 cg06654109 CHAT cg14722693 CSGALNACT1 cg17482089 PCDHA6 cg06555246 CTXN1 cg17482089 PCDHA2 cg15231202 APOD cg17482089 PCDHA1 cg23720929 APOD cg17482089 PCDHA7 cg14936961 MSLN cg17482089 PCDHA5 cg09575189 SYT8 cg17482089 PCDHA11 cg11702866 HRH1 cg17482089 PCDHA3 cg14894368 NAPSA cg17482089 PCDHA4 cg16680624 NOD1 cg17482089 PCDHA8 cg06646796 DHRS2 cg01335367 MGC 14436 cg05924780 LPAR5 cg12981354 GNAL cg23704150 CAPN8 cg03656583 SGK1 cg00633815 RXRG cg00548268 NPTX2 cg13818688 PAPPA2 cg04761722 GRIN2A cg10773212 FAM178B cg11132045 PC

[0845] cg18057851 NRCAM cg24635109 IPCEF1 cg18640660 STRA6 cg15366423 MAPK4 cg20709371 SERPINA3 cg25371081 IQSEC3 cg04677882 NPW cg04522339 PAX5 cg26338105 EPS8 cg16919569 NBLA00301 cg20036207 ICAM4 cg16919569 HAND2 cg20036207 ICAM1 cg03118854 PNPLA5 cg00797102 KCNJ2 cg02843201 VWA3A cg06457736 HRH1 cg11076487 SORBS2 cg13959356 BID cg12197399 SORBS2 cg23698062 RPSAP52 cg00545229 C1 QTNF7 cg23698062 HMGA2 cg25234117 PLCH1 cg13235717 SFN cg21752601 UCN cg24426405 AGR2 cg22249529 KLHDC7B cg08622198 CHRM3 cg14993439 UHRF1 cg21837069 VWA1 cg12179939 RAP1GAP cg01507563 RGS5 cg25959149 BANF2 cg08443014 LYPD6 cg01326836 HOXA11 cg00714130 SOX5 cg16191009 CPNE9 cg07435254 SNX25 cg22273285 HOXB6 cg14476479 GPT cg01227537 ZIC1 cg16587265 GPT cg02510853 PKMYT1 cg27398499 TACSTD2 cg14456683 ZIC1cg05954120 TMEM79 cg17280346 ZIC1 cg11319300 NTM cg25212131 SLC25A2 cg11569121 RAB27A cg03876116 ADAM33 cg16505204 THRSP cg22689630 KLK13 cg13778865 IPCEF1 cg13578841 CGNL1 cg19709585 LAMB3 cg16986720 DLG2 cg07435386 KLHDC8A cg02712145 CRHR2 cg08887622 SEMA3F cg00028336 SORBS2 cg26176246 SLC34A2 cg18628164 MAP2 cg01406381 SLC1A5 cg04072430 S100A3 cg17240725 WBSCR26 cg05099221 PRR15L cg08328750 KRT15 cg23268208 HORMAD2 cg19836199 RUNX1 cg21273765 HIF3A cg15714343 MFAP4 cg08618219 SEMA6D cg05893475 RAET1E cg05113927 UCN cg08099797 LIPH cg06182018 LIFR cg06754890 UNC5CL cg08844745 FBXO17 cg14562054 RIN1 cg13080565 GRIK2 cg12521187 LINGO1 cg01876130 FER1L5 eg 17874802 LPAR5 cg09206268 MATN4 cg06912965 RUNX1 cg14904619 GPRIN1 cg21544658 SYT8 cg20442078 UCN cg12746703 C6 cg17546247 ZIC4 cg01799015 PALM cg17546247 ZIC1 cg25883083 LPAR5 cg03415532 B3GALT2 cg25452264 TPD52L1 cg21183011 SERPIND1 cg15761609 FAM20A cg17632937 HORMAD2 cg13103915 CSF2 cg13107973 KCNK2 cg24331162 SYT8 cg07421287 KCNA3 cg13420744 DHRS3 cg21810373 SHANK2 cg18717554 STRA6 cg19591003 C14orf86 cg21201572 AGR2 cg25415941 DLG2 cg27235662 CLDN16 cg03143395 SORBS2 cg21950166 SFN cg19448837 MAPK4 cg04042305 RAB27A cg04053572 RGS4 cg00280812 KRT80 cg09905635 TNS3 cg25742326 PLEKHN1 cg16000321 CRYAB cg23505644 CDSN cg25711779 EFEMP1 cg23220439 RIN1 cg08239282 FHL2 cg01136183 TM7SF4 cg22274074 HOXA11 cg23168040 KLK10 cg15745560 NUPR1 cg24306984 C10orf67 cg10046892 PAQR6 cg21538450 KLHDC8A cg26179497 ACSL6 cg17605371 ETNK2 cg10521706 OSR1 cg10981439 TREM1 cg04135144 TEDDM1 cg05082205 CLDN16 cg12426652 CLDN14 cg10352418 CA7 cg19308222 EREG cg08443845 RUNX1 cg10509254 NBLA00301cg02196805 CSF2 cg10509254 HAND2 cg04701290 BID cg02582940 FER1L5 cg24237576 DNASE1L2 cg13899108 PDE4C cg08964563 MET cg23547429 SLC43A3 cg20457051 GPR115 cg00635674 GRAMD2 cg05812089 LOC100126784 cg01295399 SLC18A3 cg25286393 NAPSA cg01295399 CHAT cg14841694 SEMA3F cg15808668 IP6K3 cg07196577 MYBPHL cg13305951 RNF112 cg04242992 RAB27A cg05565668 SLC25A15 cg21875969 DOK7 cg04939673 C4orf38 cg06073141 CLDN16 cg00576132 GNAL cg21328082 TREM1 cg04498349 GRIN2A cg06438797 CA7 cg07766263 SLC7A14 cg08739221 EPS8 eg 10978403 FLRT1 cg01767040 CHRM3 cg04305670 MATN4 cg05239841 ARNTL cg18392348 TFF3 cg05798299 RAB27A cg03190099 FAHD2A cg02902865 DHRS3 cg09486642 TRPV6 cg05346878 TACSTD2 cg02801193 SULT1C4 cg21301496 RNASE11 cg21843594 HORMAD2 cg15974698 TGM1 cg21636683 SLC25A2 cg24428058 CLU cg18608703 CFTR cg01441698 PLA2G2E cg26267081 C14orf86 cg04071779 NOD1 cg22826333 CWH43 cg17324980 ICAM1 cg02412684 FGFR2 cg17324980 ICAM4 cg04398720 TCL1A cg09135399 APOBEC3H cg12419052 ZIC1 cg13888509 STRA6 cg16845108 CECR2 cg16510709 LPAR5 cg00447208 TMC6 cg12942329 TMEM79 cg15656257 ABCB9 cg16657886 SYT8 cg12125606 CXCL13 cg11865061 MSLN cg01319940 RXRG cg15442792 MUC21 cg25928446 CADM3 cg13286116 ARNTL cg23744654 PLD5 cg07381930 GOLGA7B cg07175433 KCNIP4 cg07538160 MUC21 cg23291270 EPS8 cg01420517 ETNK2 cg16960593 LPA cg00176348 RNASE11 cg26104206 PHYHIP cg04654868 TPD52L1 cg22717825 PLA2G2E cg02579816 CST5 cg21501207 SH2D1 B cg19856444 SLC39A12 cg02196805 CSF2 cg21536783 TACSTD2 cg25886621 PLEKHG7 cg07189401 LOC100126784 cg05893475 RAET1E cg26733338 SYT1 cg17874802 LPAR5 cg14042308 RIN1 cg05999628 CCL11 cg23266464 ARNTL cg09928596 ACTBL2 cg26037504 RDH5 cg04531728 FAM135Bcg02294791 SYT8 cg10713220 ADAH

[0846] cg01716680 GJA4 cg21275690 HCFC1R1 cg00721951 SEMA3F cg09006659 CLU cg09022733 SERPINA1 eg 13745346 CBFA2T3 cg25881591 NTM cg01246266 TMC6 cg07978005 PTPRE eg 12966875 SLPI

[0847] cg10709671 CTXN1 cg23501051 L0C100126784 cg06720467 SFN cg10780897 RP1

[0848] cg14246614 PLA2G2E eg 16274229 LRFN5 cg03190661 TMC6 cg25203007 GALE

[0849] cg13333107 EPS8 cg21538450 KLHDC8A cg09723979 LOC100126784 cg13471990 ENTPD1 cg13171103 LPAR5 cg19856444 SLC39A12 cg09404790 ARNTL eg 19907725 S100A2 cg11154784 TGM1 cg17240725 WBSCR26 eg 19035595 ENG cg24507742 SLC1A5 cg26369418 CSF2 cg15442792 MUC21 cg03547181 ADORA1 cg05791813 PHACTR3 cg10918527 NOD1 cg11706865 CHST8

[0850] eg 18052274 LOC100126784 cg12653876 LGALS3 cg23133335 LOXL2 eg 14696870 FCER1A cg20826957 ETV5 cg08999895 SLC22A18AS cg15908125 SYTL1 cg17707274 MMP7

[0851] eg 14085262 CHI3L1 cg20716202 SLC22A18AS cg07999953 GJB3 cg22401524 SLC01A2 cg24565274 LOC100126784 cg09446222 CAPSL cg23291270 EPS8 cg16680624 NOD1

[0852] cg01853665 RNASE11 cg07538160 MUC21 cg06588028 LPAR5 cg10109486 MYT1L cg20943461 TMC6 cg08900770 IL16

[0853] cg24610399 ENG cg19047119 PDE4D cg24735489 CDSN cg04943586 RFX4 cg24790676 KLK10 cg03861143 SYT8 cg07429623 KLHDC8A cg00176348 RNASE11 cg01077915 SYT12 cg04452020 PDE4D cg20140201 C2orf54 cg25881591 NTM cg02325250 CSF2 cg17871234 MICAL2 cg09928596 ACTBL2 cg18165389 SLC6A1 cg21587879 DOK7 cg10981439 TREM1 cg23739048 RXRG cg00214970 SPHKAP cg06100542 RXRG cg08443845 RUNX1 cg03290977 C1 QTNF7 cg05368724 FAM178B cg03611487 LOC100126784 cg06130787 KLK10

[0854] eg 17243643 RDH5 eg 18433423 CDH13

[0855] cg06130787 KLK10 cg07696914 TG

[0856] cg08512745 PC cg12981901 KCNIP4 cg03448202 MYBPH cg00753658 BIRC7 cg03070999 LPAR5 cg00045164 AGR2cg01870058 SYTL1 cg15311897 IPCEF1 cg26648185 TIAM1 cg09404790 ARNTL cg12125606 CXCL13 cg16510709 LPAR5 cg15161667 LRP4 cg00825734 SCNN1D eg 13982340 SYTL1 cg07435386 KLHDC8A cg06501964 ARNTL cg16983159 TMEM173 cg03733371 LIPH cg05848175 S100A9 cg22681654 ABR cg11456828 KCNA2 cg08657449 TM7SF4 cg01406381 SLC1A5 eg 19009962 FAM178B cg27397879 FAM20A cg05301052 ACTBL2 cg00714130 SOX5 eg 17888803 FAM178B cg04326298 PDE9A cg23501051 LOC100126784 cg07326586 UBD cg11052439 GJA4 cg24428058 CLU cg23620049 LIPH cg07232945 ADORA1 cg03861143 SYT8 cg05301052 ACTBL2 cg22717825 PLA2G2E cg27086014 SIGLEC12 cg25203007 GALE eg 12304663 KIR3DL2 cg20127859 C8orf73 cg08512745 PC

[0857] cg02192520 RDH5 eg 17302948 LDLRAD3 eg 15073665 PC cg27321439 SLC12A5 cg00045164 AGR2 cg05638787 C15orf62 eg 18338296 THRSP eg 18252433 CHST8 cg02480666 GJA4 cg09013287 CTSB cg09830083 P4HA2 cg21301496 RNASE11 cg00988350 GALE cg07337602 ST6GAL1 cg18122696 SYT8 cg11190327 SLC14A2 cg03074223 LPAR5 cg06196379 TREM1 cg23180713 FAM178B eg 13982340 SYTL1 cg27111077 ENTPD1 cg20127859 C8orf73 cg04709551 ETV5 cg15236354 GREM2 eg 16389610 SERPINA1 cg07554745 LPA cg05368724 FAM178B cg09925075 RTN1 cg26941211 PRSS23 cg09135399 APOBEC3H cg21308767 EPS8 eg 16329896 TNS3 cg13591783 ANXA1 cg00541854 GPR142 cg09968361 SERPINA1 cg02192520 RDH5 cg02392359 RAB27A cg05013507 DUSP14 cg02126235 SERPINA1 eg 18660898 MLLT11 cg15761609 FAM20A cg13286116 ARNTL cg25883083 LPAR5 cg27235662 CLDN16 cg01002857 CTSB eg 19996603 CREB5 cg09275205 GPM6A eg 10447638 S100A4 cg00988350 GALEcg01724702 SIRPG cg07689834 S100A1 cg23180713 FAM178B cg21188037 MYO1G cg12919039 DPP10 cg22681654 ABR cg05859264 MAPK13 cg10918527 NOD1 cg09089187 KCNIP1 cg04071779 NOD1 cg15908125 SYTL1 eg 17267976 SRD5A2 cg07774533 LDLRAD3 cg20749059 SPATA12 cg00797102 KCNJ2 cg07196577 MYBPHL cg14515826 TNS3 cg11832805 FAM19A5 cg18305394 PCK1 cg13778865 IPCEF1 cg25788793 SLC2A9 cg08043750 RETN cg03190661 TMC6 cg15161667 LRP4 cg25638829 FAM19A5 cg13171103 LPAR5 cg25286393 NAPSA cg03733371 LIPH cg05966699 C3orf52 cg03290977 C1 QTNF7 cg01870058 SYTL1 cg04766084 SPATA12 cg03214064 C12orf42 eg 19953634 KCNJ13 cg09087222 CTSB cg06501964 ARNTL cg06588028 LPAR5 cg26336350 C1orf116 cg06935531 TNS3 cg23080769 FMO1 cg01853665 RNASE11 cg02088943 CTSB cg06048910 SLC22A18AS eg 16389610 SERPINA1 cg09968361 SERPINA1 cg24562785 SIRPG eg 15073665 PC

[0858] cg15625172 LRFN5cg24505973 KRT85

[0859] cg13591783 ANXA1 cg11702403 TNFRSF10B cg07381930 GOLGA7B cg03070999 LPAR5 cg09975110 IL16 cg23698062 RPSAP52 cg23698062 HMGA2 cg04877695 FCER1A cg02126235 SERPINA1 cg23699588 CYP19A1 eg 18338296 THRSP cg25700533 C15orf62 cg03547181 ADORA1 cg23620049 LIPH cg19104050 CMKLR1 cg16505204 THRSP cg09830083 P4HA2 cg06244990 TNFRSF10B cg20943461 TMC6 cg23739048 RXRG eg 19737327 CREB5 cg27111077 ENTPD1 cg01418424 MACC1 cg20826957 ETV5 cg03448202 MYBPH eg 17888803 FAM178B eg 18458509 SLC22A18AS cg02392359 RAB27A cg20616883 TMEM117 cg04709551 ETV5 cg03074223 LPAR5 cg16779321 CREB5 cg26941211 PRSS23 cg23190089 SLC22A18AS cg21308767 EPS8

[0860] Interestingly, we found a negative correlation between DNA methylation and gene expression affecting 1080 CpGs - mainly located in the 5’ region - most of which showed CpG hypomethylation and gene upregulation (Figure 5). We also found positive correlations between DNA methylation and gene expression, mainly involving CpGs located in gene bodies (Table 13 and Table 14).Table 13. Correlation between the differentially methylated CpGs in our discovery series associated with differentially expressed genes in The Cancer Genome Atlas (TCGA).

[0861]

[0862] negative correlation positive correlation

[0863]

[0864] hypo - up hyper - down hypo - down hyper - up 5' region 262 11 18 13

[0865] gene body 251 12 44 40

[0866]

[0867] TOTAL 513 23 62 53

[0868]

[0869] negative correlation positive correlation hypo - up hyper - down hypo - down hyper - up 5' region 148 116 56 81

[0870] gene body 176 104 125 138

[0871]

[0872] TOTAL 324 220 181 219

[0873] Hypo: hypomethylated. Hyper: hypermethylated. Up: upregulated. Down: downregulated.

[0874] Table 14. Gene ontology (GO) enrichment analysis of the differentially expressed genes in The Cancer Genome Atlas (TCGA) associated with differentially methylated CpGs at their 5' region in our discovery series.

[0875] Low-risk PTC vs NT

[0876] genesEP Fold

[0877] IC - Cou Exp Enric P- REFLIST nt ecte hmen valu Reduce GO term (17968) (165) d t e FDR d Term secretion 4.83 2.51 secretio (G0:0046903) 498 17 4.57 3.72 E-06 E-02 n cell Biological process cell adhesion 1.51 2.35 adhesio (G0:0007155) 921 25 8.46 2.96 E-06 E-02 n develop tissue development 15.1 2.58 2.01 mental (G0:0009888) 1644 35 0 2.32 E-06 E-02 process basal plasma

[0878] membrane 1.10 2.47 membra (G0:0009925) 243 10 2.23 4.48 E-04 E-02 ne basal basal part of cell 1.93 3.54 part of (G0:0045178) 261 10 2.40 4.17 E-04 E-02 cell cell junction 18.5 5.67 2.29 cell (G0:0030054) 2022 37 7 1.99 E-05 E-02 junction extracell extracellular vesicle 18.2 8.42 2.84 ular Cellular component (GO:1903561) 1983 36 1 1.98 E-05 E-02 region extra cell extracellular organelle 18.2 8.45 2.44 (G0:0043230) 1984 36 2 1.98 E-05 E-02 extracellular

[0879] membrane-bounded

[0880] organelle 18.2 8.45 2.14 (G0:0065010) 1984 36 2 1.98 E-05 E-02 extracellular exosome 18.0 1.36 2.75

[0881]

[0882] (G0:0070062) 1962 35 2 1.94 E-04 E-02extracell extracellular region 36.1 4.07 4.12 ular (G0:0005576) 3941 68 9 1.88 E-08 E-05 region extra cell extracellular space 28.0 1.51 7.62 ular (G0:0005615) 3050 51 1 1.82 E-05 E-03 region cell cell periphery 52.0 1.96 3.96 peripher (G0:0071944) 5671 90 8 1.73 E-09 E-06 plasma membrane 47.8 1.63 1.10 membra

[0883]

[0884] (G0:0005886) 5208 80 3 1.67 E-07 E-04

[0885]

[0886] High-risk PTC vs NT

[0887] genesEP Fold

[0888] IC - Cou Exp Enric P- Redu REFLIST nt ecte hmen valu ced GO term (17968) (259) d t e FDR Term homophilic cell

[0889] adhesion via plasma

[0890] membrane adhesion cell molecules 2.64 4.11 adhes Biological process (G0:0007156) 164 14 2.36 5.92 E-07 E-03 ion cell-cell adhesion via

[0891] plasma-membrane cell adhesion molecules 4.20 3.27 adhes (G0:0098742) 257 17 3.70 4.59 E-07 E-03 ion calciu calcium-dependent m ion protein binding 3.31 4.06 (G0:0048306) 79 8 1.14 7.03 E-05 E-02

[0892] calciu m ion calcium ion binding 10.0 1.17 5.75 bindin (G0:0005509) 697 34 5 3.38 E-09 E-06

[0893] signal Molecular function ing recep cation transmembrane tor transporter activity 1.56 3.82 activit (G0:0008324) 616 24 8.88 2.70 E-05 E-02 y signal ing recep ion transmembrane tor transporter activity 11.5 2.65 4.33 activit (G0:0015075) 803 28 7 2.42 E-05 E-02 y cell calyx of Held 2.38 2.09 juncti (G0:0044305) 18 4.26 15.42 E-04 E-02 on potassium channel

[0894] complex 9.08 1.83 mem (G0:0034705) 87 9 1.25 7.18 E-06 E-03 brane integral component of

[0895] presynaptic

[0896] membrane 1.53 1.72 mem (G0:0099056) 74 7 1.07 6.56 E-04 E-02 brane voltage-gated

[0897] potassium channel

[0898] complex 1.92 1.77 mem Cellular component (G0:0008076) 77 7 1.11 6.31 E-04 E-02 brane intrinsic component of

[0899] presynaptic

[0900] membrane 2.57 2.17 mem (G0:0098889) 81 7 1.17 6.00 E-04 E-02 brane cell axon terminus 2.84 2.30 juncti (G0:0043679) 110 8 1.59 5.05 E-04 E-02 on neuron projection cell terminus 1.21 1.63 projec (G0:0044306) 124 9 1.79 5.04 E-04 E-02 tion cation channel

[0901] complex 4.12 1.19 mem

[0902]

[0903] (G0:0034703) 210 14 3.03 4.62 E-06 E-03 branepresynaptic

[0904] membrane 3.46 2.50 mem (G0:0042734) 144 9 2.08 4.34 E-04 E-02 brane integral component of

[0905] synaptic membrane 4.59 3.09 mem (G0:0099699) 150 9 2.16 4.16 E-04 E-02 brane ion channel complex 5.77 1.46 mem (G0:0034702) 281 16 4.05 3.95 E-06 E-03 brane intrinsic component of

[0906] synaptic membrane 7.75 4.61 mem (G0:0099240) 162 9 2.34 3.85 E-04 E-02 brane cell distal axon 1.75 1.77 juncti (GO:0150034) 264 13 3.81 3.42 E-04 E-02 on neuronal cell body 9.38 3.79 cell (G0:0043025) 481 23 6.93 3.32 E-07 E-04 body synaptic membrane 2.99 5.03 mem (G0:0097060) 360 17 5.19 3.28 E-05 E-03 brane postsynaptic

[0907] membrane 4.58 3.19 mem (G0:0045211) 255 12 3.68 3.26 E-04 E-02 brane transmembrane

[0908] transporter complex 7.88 1.14 mem (GO:1902495) 352 16 5.07 3.15 E-05 E-02 brane cell body 2.29 7.70 cell (G0:0044297) 547 24 7.88 3.04 E-06 E-04 body transporter complex 1.63 1.74 mem (GO:1990351) 376 16 5.42 2.95 E-04 E-02 brane cellul ar somatodendritic anato compartment 11.7 1.47 2.70 mica I (G0:0036477) 815 29 5 2.47 E-05 E-03 entity cell dendrite 3.30 2.57 projec (G0:0030425) 595 21 8.58 2.45 E-04 E-02 tion cell dendritic tree 3.36 2.51 projec (G0:0097447) 597 21 8.61 2.44 E-04 E-02 tion cell postsynapse 6.45 4.07 juncti (G0:0098794) 594 20 8.56 2.34 E-04 E-02 on plasma membrane 17.2 7.42 1.67 mem region (G0:0098590) 1193 38 0 2.21 E-06 E-03 brane integral component of

[0909] plasma membrane 23.1 1.46 9.85 mem (G0:0005887) 1603 51 1 2.21 E-07 E-05 brane intrinsic component of

[0910] plasma membrane 24.2 1.73 8.75 mem (G0:0031226) 1680 52 2 2.15 E-07 E-05 brane cell neuron projection 18.9 6.04 9.39 juncti (G0:0043005) 1317 38 8 2.00 E-05 E-03 on cell synapse 18.4 6.36 4.15 juncti (G0:0045202) 1282 34 8 1.84 E-04 E-02 on cell cell junction 29.1 7.14 4.37 juncti (G0:0030054) 2022 48 5 1.65 E-04 E-02 on plasma membrane 75.0 5.78 1.17 mem (G0:0005886) 5208 117 7 1.56 E-08 E-04 brane cell cell periphery 81.7 5.88 5.94 perip (G0:0071944) 5671 124 4 1.52 E-08 E-05 hery intrinsic component of

[0911] membrane 77.3 1.21 1.53 mem (G0:0031224) 5365 107 3 1.38 E-04 E-02 brane integral component of

[0912] membrane 75.2 1.45 1.73 mem

[0913]

[0914] (G0:0016021) 5218 104 1 1.38 E-04 E-02GO analyses of the differentially expressed genes associated with differentially methylated CpGs at their 5’ region showed an enrichment in BP cell adhesion, MF calcium-dependent processes, and CC membrane, both in non-mPTC and mPTC, as well as a specific enrichment in the CC extracellular region in non-mPTC, confirming our results in the discovery series.

[0915] Taken together, these findings provide evidence for the regulatory role of DNA methylation in the metastatic progression of thyroid cancer.

[0916] Identification of a unique DNA methylation signature associated with mDTC Notwithstanding the different DNA methylomes associated with PTC and FTC, we identified a signature of 156 CpGs that were differentially methylated between non-mDTC and mDTC, 69.2% of which were hypomethylated in mDTC. An unsupervised hierarchical cluster analysis of this 156-CpG signature clearly separated samples into two clusters independently of histology and BRAF and RAS mutations: one cluster contained all the normal tissue and most of the non-mDTC (28 / 30), while the other contained most of the mDTC (31 / 34) and all the LNM and DM derived from mDTC, independently of TERTp mutation status (Figure 6A).

[0917] Importantly, the metastases did not form a distinct subcluster and most of them (10 / 15) clustered next to their paired primary tumor (Figure 17), indicating that the DNA methylation alterations were generated in primary tumors and maintained in metastatic lesions. Intriguingly, the primary tumors of the three patients with metachronous DM grouped with the primary tumors of patients with synchronous DM, suggesting that the 156-CpG signature could potentially predict the risk of developing DM. The average DNA methylation of the 156 CpGs was lower in mDTC and DM than in non-mDTC and normal tissue, and was able to discriminate between non-mDTC and mDTC with an area under the curve (AUC) of 0.835, sensitivity of 86.7% and specificity of 68.6% (Figure 6B-C).The area under a receiver operating characteristic (ROC) curve (AUC) and their sensitivity and specificity was calculated for each of the individual 156 CpG biomarkers of the gene signature (Table 15).

[0918] Table 15. ROC Curve values of the 156 differentially methylated CpGs between non-mDTC and mDTC. Criterion: Youden

[0919] EPIC array ID Location (Feb. 2009 AUC

[0920] (GRCh37 / hg19)) (confidence interval) Sensitivity Specificity cg26785499 chr1:2820433-2820433 0.822 (0.719, 0.924) 0.867 0706 cg01431482 chr1:2989085-2989085 0.726 (0.601, 0.852) 0.647 0.8 cg26139866 chr1:4239037-4239037 0.805 (0.694, 0.915) 0.767 0.824 cg14035553 chr1:21029443-21029443 0.783 (0.67, 0.896) 0.967 0.529 cg09551984 chr1:30229976-30229976 0.793 (0.683, 0.904) 0.883 0.647 cg24364593 chr1:34533065-34533065 0.782 (0.666, 0.899) 0.967 0.559 cg06437931 chr1:38180356-38180356 0.749 (0.627, 0.871) 0.833 0.618 cg10776919 chr1:50889442-50889442 0.747 (0.624, 0.87) 0.588 0.9 cg18279094 chr1:63790044-63790044 0.841 (074, 0.942) 0.824 0.8 cg04360793 chr1:79472361-79472361 0742 (0.619, 0.866) 0.618 0.833 cg15084543 chr1:79472408-79472408 0768 (0.651, 0.884) 0.676 0.833 cg21618157 chr1:92202671-92202671 0.737 (0.615, 0.859) 0.765 0.633 cg07530798 chr1:119535589-119535589 0.79 (0.679, 0.901) 0.588 0.933 cg24398933 chr1:151492569-151492569 0.776 (0.657, 0.895) 0.647 0.9 cg19387132 chr1:158254116-158254116 0.795 (0.681, 0.909) 0.933 0.676 cg27534912 chr1:158258802-158258802 0.781 (0.667, 0.895) 0.867 0.647 cg14275423 chr1: 158656484-158656484 0766 (0.648, 0.884) 0.8 0.676 cg05304729 chr1:158800024-158800024 0.839 (0.743, 0.935) 1 0.559 cg25743622 chr1:223403446-223403446 0.821 (0.715, 0.927) 0.8 0.824 cg26404722 chr2:9353201-9353201 0.74 (0.617, 0.863) 0.5 0.933 cg11448675 chr2:9375570-9375570 0727 (0.604, 0.851) 0.529 0.833 cg15366555 chr2:74601926-74601926 079 (0.671, 0.91) 0765 0.867 cg11412468 chr2:84433405-84433405 0.801 (0.691, 0.911) 0.8 0.735 cg14231073 chr2:121347537-121347537 0.749 (0.624, 0.874) 0.618 0.833 cg22540067 chr2:125305001-125305001 0791 (0.678, 0.904) 0.9 0.618 cg21445398 chr2:125425654-125425654 0.818 (0.709, 0.926) 0.767 0.824 cg05053811 chr2:126401437-126401437 0.775 (0.658, 0.893) 1 0.588 cg12906989 chr2: 137925983-137925983 0.806 (0.699, 0.913) 0.867 0.676 cg12341429 chr2:161887154-161887154 0.75 (0.626, 0.874) 0.853 0.633 cg18621766 chr2: 194689205-194689205 0.801 (0.689, 0.913) 0733 0.853 cg04314280 chr2:218687427-218687427 0.767 (0.648, 0.885) 0.765 0.733 cg15446845 chr3:5470909-5470909 0754 (0.631, 0.877) 0.967 0.559 cg22438338 chr3:14634313-14634313 0.778 (0.665, 0.892) 0.867 0.647 cg18780425 chr3:69250024-69250024 0763 (0.642, 0.884) 0.676 0.9 cg11963062 chr3: 100633394-100633394 0781 (0.67, 0.892) 0.529 0.933 cg16768018 chr3:147108843-147108843 0.8 (0.685, 0.915) 0.647 0.933 cg04871469 chr4:2257119-2257119 0.702 (0.572, 0.832) 0.529 0.933 cg25175009 chr4:63265254-63265254 0.827 (0.729, 0.926) 0.967 0.559 cg18549952 ch r4: 129762662-129762662 0764 (0.648, 0.88) 0.676 0767 cg04783185 chr4:190521538-190521538 0.792 (0.678, 0.907) 0.8 0.765 cg01653189 chr5:2569982-2569982 0.779 (0.663, 0.896) 0.8 0.706 cg07826516 chr5:3106230-3106230 0.821 (0.717, 0.924) 0.767 0.833 cg11925561 chr5:3304588-3304588 0.896 (0.819, 0.973) 0.833 0.853 cg11328079 chr5:3504656-3504656 0.861 (0.77, 0.951) 0.833 0765 cg26976756 chr5:3510004-3510004 0.81 (0.706, 0.914) 0.9 0.618 cg18916149 chr5:3532269-3532269 0.855 (0.758, 0.952) 0.867 0.794 cg20009810 chr5:3537215-3537215 0.852 (0.754, 0.95) 0.933 0.735 cg09696091 chr5:3564465-3564465 0.766 (0.648, 0.884) 0.8 0.676

[0921]

[0922] cg13453847 chr5:5034274-5034274 0.834 (0731, 0.937) 0.9 0706cg03290223 chr5:5063319-5063319 0.781 (0.668, 0.895) 0.933 0.559 cg06587623 chr5:5079159-5079159 0.8 (0.69, 0.91) 0.733 0.735 cg14672100 chr5:6344259-6344259 0.804 (0.696, 0.912) 0.967 0.618 cg20997704 chr5:13143814-13143814 0.807 (0.693, 0.92) 0.967 0.618 cg00086247 chr5:29568333-29568333 0.795 (0.685, 0.905) 0.8 0.735 cg17738861 chr5:32233621-32233621 0.758 (0.636, 0.88) 0.647 0.833 cg04206637 chr5:63202277-63202277 0.829 (0.728, 0.93) 0.9 0.735 cg07497602 chr5:105881003-105881003 0.789 (0.678, 0.9) 0.767 0.735 cg01664151 chr5:132945292-132945292 0.781 (0.671, 0.892) 0.8 0.647 cg20876760 chr6:9399379-9399379 0.741 (0.618, 0.865) 0.676 0.8 cg00237606 chr6:29455256-29455256 0.852 (0.762, 0.942) 0.833 0.735 cg22424903 chr6:35454434-35454434 0.789 (0.679, 0.899) 0.618 0.933 cg24749265 chr6:41216090-41216090 0.799 (0.69, 0.908) 0.933 0.618 cg01378044 chr6:85462137-85462137 0.897 (0.819, 0.975) 0.867 0.882 cg25901805 chr6:144286162-144286162 0.85 (0.757, 0.943) 0.706 0.9 cg25075684 chr7:2933297-2933297 0.798 (0.684, 0.912) 0.867 0.647 cg08241694 chr7:50633896-50633896 0.798 (0.689, 0.907) 0.9 0.618 cg16455187 chr7:50635841-50635841 0.766 (0.646, 0.885) 1 0.529 cg06915574 chr7:76639592-76639592 0.763 (0.643, 0.882) 0.794 0.667 cg01611504 chr7:108525503-108525503 0.825 (0.721, 0.93) 0.867 0.706 cg10589286 chr7: 154916724-154916724 0.85 (0.756, 0.944) 0.933 0.706 cg15983026 chr7: 155325552-155325552 0.779 (0.663, 0.896) 0.882 0.633 cg03974619 chr7: 156851553-156851553 0.822 (0.72, 0.924) 0.933 0.618 cg06437651 chr7:157466201-157466201 0.791 (0.677, 0.905) 0.833 0.706 cg15819924 chr7: 157886456-157886456 0.787 (0.675, 0.899) 0.867 0.618 cg06593940 chr7: 157950633-157950633 0.749 (0.629, 0.869) 0.833 0.588 cg05647867 chr7: 158244818-158244818 0.778 (0.665, 0.892) 1 0.529 cg07178994 chr8:816998-816998 0.745 (0.623, 0.867) 0.867 0.618 cg07340894 chr8:8238774-8238774 0.738 (0.614, 0.863) 0.912 0.533 cg14793063 chr8:30557532-30557532 0.759 (0.639, 0.879) 0.529 1 cg15325961 chr8:49183143-49183143 0.806 (0.7, 0.912) 0.967 0.559 cg10780897 chr8:55528099-55528099 0.794 (0.677, 0.911) 0.8 0.765 cg17406248 chr9:97807693-97807693 0.765 (0.644, 0.886) 0.794 0.733 cg13811240 chr9:97851374-97851374 0.762 (0.643, 0.88) 0.599 0.933 cg04027576 chr9:117996948-117996948 0.795 (0.678, 0.913) 0.733 0.824 cg19537600 chr9:119836588-119836588 0.802 (0.692, 0.912) 0.567 0.971 cg16681908 chr9:140782525-140782525 0.788 (0.675, 0.902) 0.9 0.647 cg03478630 chr10:10436568-10436568 0.873 (0.785, 0.96) 0.733 0.882 cg11601663 chr10:15666820-15666820 0.814 (0.703, 0.925) 0.933 0.676 cg09918581 chr10:68812190-68812190 0.853 (0.759, 0.947) 0.767 0.853 cg15438734 chr10:77811027-77811027 0.731 (0.604, 0.859) 0.824 0.6 cg20069378 chr10:98117247-98117247 0.773 (0.657, 0.889) 0.933 0.588 cg27567836 chr10:124091022-124091022 0.822 (0.716, 0.927) 0.667 0.912 cg11174855 chr10:134598352-134598352 0.769 (0.65, 0.888) 0.647 0.833 cg08264839 chr11:3191337-3191337 0.765 (0.648, 0.882) 0.867 0.647 cg05133853 chr11:5265008-5265008 0.786 (0.666, 0.907) 0.9 0.676 cg10092198 chr11:6129548-6129548 0.748 (0.626, 0.87) 0.7 0.765 cg16573636 chr11:19917499-19917499 0.714 (0.58, 0.847) 0.824 0.633 cg18157275 chr11:56345140-56345140 0.848 (0.755, 0.941) 0.8 0.794 cg04562441 chr11:56380921-56380921 0.858 (0.764, 0.952) 0.8 0.853 cg22151881 chr11:68082621-68082621 0.749 (0.627, 0.871) 0.676 0.8 cg05355757 chr12:16762831-16762831 0.758 (0.638, 0.878) 0.667 0.794 cg14096615 chr12:16762843-16762843 0.77 (0.65, 0.889) 0.667 0.853 cg19104050 chr12:108701064-108701064 0.797 (0.685, 0.909) 0.767 0.765 cg26998150 chr12:130734158-130734158 0.815 (0.702, 0.927) 0.967 0.647 cg15858166 chr12:130953551-130953551 0.795 (0.684, 0.906) 0.767 0.765 cg06626892 chr13:78995230-78995230 0.841 (0.743, 0.939) 0.867 0.706 cg24690071 chr13:100635352-100635352 0.751 (0.627, 0.875) 0.559 0.967 cg23117999 chr14:96508302-96508302 0.815 (0.707, 0.923) 0.867 0.706 cg23450358 chr14:97581404-97581404 0.841 (0.74, 0.942) 0.833 0.765 cg04690464 chr14:104711952-104711952 0.789 (0.671, 0.907) 0.733 0.794 cg20925031 chr15:29988135-29988135 0.767 (0.651, 0.883) 0.529 0.967 cg09262300 chr15:86519955-86519955 0.835 (0.736, 0.935) 0.933 0.676

[0923]

[0924] cg09783969 chr15:100886002-100886002 0.806 (0.7, 0.912) 0.7 0.794cg01252526 chr16:711033-711033 0.763 (0.645, 0.88) 0.826 0.6 cg03287339 chr16:711135-711135 0.774 (0.656, 0.891) 0.588 0.9 cg27039467 chr16:1730561-1730561 0.809 (0.704, 0.913) 0.794 0.733 cg13985307 chr16:25435877-25435877 0.828 (0.726, 0.93) 1 0.618 cg19530424 chr16:51795292-51795292 0.733 (0.61, 0.857) 0.676 0.733 cg04770195 chr16:53114571-53114571 0.766 (0.647, 0.884) 0.735 0.733 cg04980805 chr16:82377691-82377691 0.786 (0.675, 0.897) 0.833 0.676 cg16122778 chr16:86671056-86671056 0.785 (0.671, 0.899) 0.833 0.647 cg16596367 chr16:89070757-89070757 0.765 (0.647, 0.882) 0.529 0.9 cg21851534 chr17:3907994-3907994 0.74 (0.619, 0.862) 0.588 0.833 cg14601733 chr17:6280370-6280370 0.774 (0.66, 0.887) 0.8 0.676 cg03916694 chr17:18150707-18150707 0.784 (0.668, 0.9) 0.882 0.633 cg06603761 chr17:32691699-32691699 0.771 (0.652, 0.89) 0.833 0.706 cg23368579 chr18:519265-519265 0.767 (0.644, 0.89) 0.967 0.618 cg19291618 chr18:4599319-4599319 0.816 (0.709, 0.923) 0.833 0.706 cg15832062 chr19:3282509-3282509 0.79 (0.68, 0.9) 0.767 0.706 cg19058865 chr19:4327350-4327350 0.741 (0.618, 0.864) 0.529 0.933 cg19839825 chr19:7616011-7616011 0.76 (0.641, 0.879) 0.647 0.867 cg20185461 chr19:18899082-18899082 0.798 (0.691, 0.905) 0.559 0.933 cg19386774 chr19:22467847-22467847 0.764 (0.648, 0.88) 1 0.471 cg08584061 chr19:22718926-22718926 0.753 (0.629, 0.877) 1 0.529 cg17699869 chr19:29421611-29421611 0.838 (0.739, 0.937) 0.767 0.824 cg20804418 chr19:32614622-32614622 0.795 (0.681, 0.909) 0.967 0.588 cg22836826 chr19:42194016-42194016 0.79 (0.678, 0.902) 0.8 0.676 cg11083325 chr19:49244668-49244668 0.793 (0.676, 0.91) 0.882 0.667 cg16509509 chr19:51151383-51151383 0.794 (0.684, 0.904) 0.933 0.588 cg26324258 chr20:40612816-40612816 0.787 (0.669, 0.905) 0.733 0.824 cg21939911 chr20:44644971-44644971 0.79 (0.675, 0.906) 0.8 0.765 cg25535351 chr20:50956552-50956552 0.763 (0.645, 0.881) 0.9 0.588 cg02762107 chr20:51335826-51335826 0.796 (0.68, 0.912) 0.9 0.676 cg18305394 chr20:56134788-56134788 0.819 (0.715, 0.923) 1 0.618 cg16493240 chr20:58156194-58156194 0.808 (0.698, 0.918) 0.867 0.735 cg16326851 chr20:60094855-60094855 0.855 (0.764, 0.946) 0.6 0.941 cg05198244 chr21:26732918-26732918 0.781 (0.668, 0.895) 0.7 0.794 cg02442623 chr21:30673051-30673051 0.78 (0.663, 0.898) 0.9 0.647 cg14861920 chr22:35019460-35019460 0.76 (0.638, 0.882) 0.967 0.588 cg01586506 chr22:38379506-38379506 0.729 (0.602, 0.857) 0.647 0.833 cg13210663 chr22:47790645-47790645 0.795 (0.683, 0.907) 0.867 0.676 cg04768709 chr22:48080590-48080590 0.815 (0.703, 0.927) 0.9 0.706 cg16689883 chr22:48441270-48441270 0.765 (0.647, 0.882) 1 0.529 cg27334056 chr22:48479650-48479650 0.778 (0.657, 0.9) 0.867 0.735 cg17659636 chr22:48581422-48581422 0.811 (0.707, 0.915) 0.933 0.588

[0925]

[0926] cg19491717 chr22:48843289-48843289 0.8 (0.69, 0.91) 0.9 0.706

[0927] We validated these results using an independent series of 13 mDTC patients, nine of whom had metachronous DM (validation series I), providing further evidence for the predictive value of the 156-CpG signature. All but one mDTC from this validation series clustered with the mDTC of the discovery series. Additionally, we analyzed three thyroid cancer invasive cell lines (TPC1, BCPAP, and HTh83) and found that they formed a subcluster within the mDTC cluster (Figure 18).

[0928] To explore the potential functional relevance of these DNA methylation alterations in the 156 CpGs, we first determined their location relative to different genomic features (Figure6D). We found that most hypomethylated CpGs were outside CGIs in intergenic regions, while hypermethylated CpGs were primarily located in CGIs at gene body regions. Interestingly, around 20% of the CpGs were located in 5’ regions and 28% in potential enhancers. GO analysis did not show significantly enriched GO terms, probably due to the low number of CpGs and associated genes (Table 16); however, some of these genes have been reported to play a role in tumorigenesis, including CCL1, DEC1, ELTD1, F0XD3, and MMP9.

[0929] Table 16. Genes associated with the CpGs in our 156-CpG signature.

[0930]

[0931] non-mPTC vs mPTC MAS1L

[0932] WDR90

[0933] TBX18

[0934] PRDM16

[0935] SOX10

[0936] C7orf66

[0937] FSTL4

[0938] BACH1

[0939] LINC01020

[0940] FLIIDEC1

[0941] TNS1

[0942] ELTD1

[0943] OR5M1

[0944] LOC284930

[0945] CHD9

[0946] MXD4

[0947] HBBP1

[0948] MNDAPTPRN2

[0949] CCL1

[0950] RNF219-AS1

[0951] DTX2P1-UPK3BP1-PMS2P11

[0952] PRAGMINDDCFLJ42289

[0953] LRRTM3

[0954] OR56B4

[0955] DMRTA2

[0956] RP1

[0957] IZUMO1

[0958] NKX6-2LINC01017 ASAP2 ITGA8 ABI3BP THSD7B LOC340094 KIF17 SPTA1

[0959] GSR DCTN1 C10orf11 CELF5 RIMBP2 CNPY1 CDH4 PHACTR3 C19orf81 NAV2

[0960] LOC100133077 ZIC4

[0961] C9orf3 MTMR12 OR5M10 FOXD3 PCK1 JADE1 FRMD4B LINC01019 STAP2 CMKLR1 ZNF729 LINC01571 ASTN2 PNPLA6 OPALIN COMP CNTNAP5 TGFBR3 ZZEF1 MMP9

[0962] LRP5 TEAD3 C14orf132 CSMD2 CGN ZIC2 TREML2P

[0963] SUSD4PLAGL1

[0964] HN1L

[0965] CD1C

[0966] BTBD16

[0967] We observed some negative correlations between gene expression and DNA methylation: for example, hypermethylation of two CpGs (cg04360793, also identified as chr1:79472361-79472361 and cg15084543, also identified as chr1:79472408-79472408) in the 5’ region of ELTD1, an orphan G-protein-coupled receptor (GPCR) of the adhesion family. ELDT1 has been related to angiogenesis in different cancer types and its upregulation has been correlated with favorable prognosis (Masiero M, et al. A Core Human Primary Tumor Angiogenesis Signature Identifies the Endothelial Orphan Receptor ELTD1 as a Key Regulator of Angiogenesis. Cancer Cell. 2013;24(2):229-241). Interestingly, these two CpGs are located inside a CGI, and analysis of the CGI showed that the hypermethylation was regional and affected the entire CGI, especially in FTC (Figure 19A-B). TCGA data confirmed the hypermethylation of the two CpGs and downregulation of ELTD1 in high-risk PTC compared to low-risk PTC (Figure 19C-D), suggesting that ELTD1 expression may be regulated through DNA methylation.

[0968] Collectively, these results indicated the potential prognostic value of the 156-CpG signature as well as the role of some CpGs in regulating cancer- related genes and the metastatic progression of thyroid cancer.

[0969] Clinical application of a CpG-based predictive model to identify DTC patients at high risk of developing DM

[0970] Our results showed the potential predictive power of the 156-CpG signature to identify DTC with a high risk of developing DM. In order to optimize the applicability of the signature to routine clinical practice, we reduced the 156-CpG signature to four CpGs that were highly methylated in normal tissue and non-mDTC but hypomethylated inmDTC, LNM and DM: cg01378044 or chr6:85,462,137-85,462, 137; cg18305394 or chr20:56, 134,788-56, 134,788; cg11925561 or chr5:3, 304, 588-3, 304, 588; and cg14672100 or chr5:6, 344, 259-6, 344, 259 (Figure 7A). cg01378044 is located in the gene body of TBX18, cg18305394 is located in the promoter region of PCK1, chr20:56, 134,788-56, 134,788; and cg11925561 or chr5:3, 304, 588-3, 304, 588 and cg14672100 or chr5:6, 344, 259-6, 344, 259 are located in separate intergenic regions of chromosome 5. Interestingly, cg11925561 or chr5:3, 304, 588-3, 304, 588 is located in a CGI shore and overlaps with a potential enhancer. This group of only four CpGs was still able to discriminate between non-mDTC and mDTC (Figure 7B).

[0971] Using these four CpGs, we generated a predictive model to determine the risk of a primary DTC developing DM. When we applied this predictive model to the 65 DTC samples included in the discovery series, using the beta values from the EPIC arrays, it correctly classified 61 of them with an accuracy of 89% (Figure 7C). In addition, in validation series I, it correctly classified 11 of 13 mDTC. We also included normal tissue and metastases in the validation analysis, and as expected, the 4-CpG predictive model indicated a low risk of DM for the former and a high risk for the latter.

[0972] We then developed simple and quantitative assays based on DNA bisulfite pyrosequencing, which is easily implementable in the clinical setting, to measure the DNA methylation of the four CpGs. The comparison of the DNA methylation levels measured by EPIC arrays and bisulfite pyrosequencing in 62 samples from the discovery series showed a strong correlation for all four CpGs (Pearson r = 0.76-0.88; P < 0.001) (Figure 20). Nevertheless, we retrained the predictive model using bisulfite pyrosequencing values to improve accuracy and obtained an adjusted predictive model that was able to correctly classify 48 of 50 DTC from the discovery series with an accuracy of 96% (Figure 7D). Finally, we validated the adjusted predictive model by analyzing the four CpGs by bisulfite pyrosequencing in the validation series II and ROCcurve analysis showed a high discriminatory performance with an accuracy of 90%, sensitivity of 81% and specificity of 91% (Figure 7E).

[0973] ITEMS OF THE INVENTION

[0974] 1. A method of determining the methylation level of a set of biomarkers in a subject suspected of having thyroid cancer, the method comprising:

[0975] а. from an extracted genomic DNA obtained from a biological sample from the subject suspected of having thyroid cancer, determining the methylation level of one or more biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), cg11925561 (located in an intergenic region of chromosome 5, chr5:3, 304, 588-3, 304, 588), and cg14672100 (located in an intergenic region of chromosome 5, chr5:6,344,259- 6,344,259) from the extracted genomic DNA.

[0976] 2. The method according to item 1, wherein step a) comprises determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PC K1 gene, chr20:56, 134,788-56, 134,788), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259.

[0977] 3. The method of any one of items 1 or 2, wherein the methylation level is determined by a quantitative assay based on bisulfite pyrosequencing to measure the DNA methylation.

[0978] 4. An in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:a. from an extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the methylation level of one or more biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), cg11925561 (located in an intergenic region of chromosome 5, chr5:3,304,588-3,304,588), and cg14672100 (located in an intergenic region of chromosome 5, chr5:6,344,259-6,344,259) from the extracted genomic DNA; and

[0979] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by a predictive model which correlates the methylation level of one or more of the biomarkers identified in step (i) with representative methylation level profiles from samples obtained from subjects having thyroid cancer previously identified as suffering from metastasis, preferably of distant metastasis (DM), said predictive model having been generated by training a computer with a plurality of methylation level profiles from previously identified subjects having thyroid cancer and suffering from metastasis, preferably of distant metastasis (DM), by machine learning on said plurality of methylation level profiles so as to obtain representative methylation level profiles associated with metastasis, preferably of distant metastasis (DM).

[0980] 5. An in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:

[0981] a. from an extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), intergenic region cg11925561 inchromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259; and

[0982] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by a predictive model which correlates the methylation level of at least the biomarkers identified in step (i) with representative methylation level profiles from samples obtained from subjects having thyroid cancer previously identified as suffering from metastasis, preferably of distant metastasis (DM), said predictive model having been generated by training a computer with a plurality of methylation level profiles from previously identified subjects having thyroid cancer and suffering from metastasis, preferably of distant metastasis (DM), by machine learning on said plurality of methylation level profiles so as to obtain representative methylation level profiles associated with metastasis, preferably of distant metastasis (DM).

[0983] 6. The method of any one of items 4 or 5, wherein the methylation level is determined by a quantitative assay based on bisulfite pyrosequencing to measure the DNA methylation.

[0984] 7. An in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:

[0985] a. from an extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the methylation level of one or more biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), cg11925561 (located in an intergenic region of chromosome 5, chr5:3,304,588-3,304,588), and cg14672100 (located in an intergenic region of chromosome 5, chr5:6,344,259-6,344,259) from the extracted genomic DNA; andb. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by (ii) comparing said methylation level / s with a reference value, wherein a deviation in the methylation level of said at least one biomarker selected from the list shown in a) with respect to said reference value, is indicative that the subject is at risk of metastasis, preferably of distant metastasis (DM).

[0986] 8. An in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:

[0987] a. from an extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259; and

[0988] b. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by (ii) comparing said methylation levels with one or more reference values, wherein a deviation in the methylation level of said biomarkers selected from the list shown in a) with respect to said one or more reference values, is indicative that the subject is at risk of metastasis, preferably of distant metastasis (DM).

[0989] 9. The method of any one of items 7 or 8, wherein the methylation level is determined by a quantitative assay based on bisulfite pyrosequencing to measure the DNA methylation.

[0990] 10. The method according to any one of the preceding items, wherein the biological sample is a sample containing DNA derived from thyroid tumor.11. The method according to item 10 wherein the sample containing DNA derived from thyroid tumor is a biopsy.

[0991] 12. In vitro method for selecting a subject having thyroid cancer at risk of metastasis, preferably of distant metastasis (DM), as a candidate to receive an adequate therapy to treat thyroid cancer, the method comprising:

[0992] a. determining whether the patient is at risk of metastasis, preferably of distant metastasis (DM), following the method of any one of items 4 to 11; and b. selecting said patient as a candidate to receive an adequate therapy to treat thyroid cancer if the patient is at risk of metastasis, preferably of distant metastasis (DM).

[0993] 13. The method according to item 12 wherein the therapy adequate for the treatment of thyroid cancer is selected from the group consisting of surgery, radiation therapy, chemotherapy, hormonal therapy or targeted therapy.

[0994] 14. The method according to items 12 or 13, wherein the targeted therapy includes one or more multikinase inhibitors selected from the group consisting of lenvatinib, sorafenib, cabozantinib, selpercatinib, pralsetinib, larotrectinib and entrectinib.

[0995] 15. A method for the treatment of thyroid cancer in a subject in need thereof comprising the administration of a therapy adequate for the treatment of thyroid cancer, wherein the patient to be treated has been identified by a method according to any of items 4 to 11.

[0996] 16. A kit, package and / or device comprising reagents adequate for implementing the methods according to any one of items 1 to 11.17. The method of any one of items 1 to 15, wherein the subject has differentiated thyroid cancer (DTC).

[0997] 18. The method of any one of items 1 to 15, wherein the subject has papillary thyroid cancer (PTC) or follicular thyroid cancer (FTC).

[0998] 19. A computer implemented method for classifying a subject having thyroid cancer as a subject at risk of developing DM, wherein the method classifies said subject on the basis of the quantitative methylation levels of one or more biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), intergenic region cg11925561 in chromosome 5 or chr5:3, 304, 588-3, 304, 588, and intergenic region cg14672100 in chromosome 5 or chr5:6, 344, 259-6, 344, 259.

[0999] 20. A computer implemented method for classifying a subject having thyroid cancer as a subject at risk of developing DM, wherein the method classifies said subject on the basis of the quantitative methylation levels of the combination of biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462,137-85,462,137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56,134,788), intergenic region cg11925561 in chromosome 5 or chr5:3,304,588-3,304,588, and intergenic region cg14672100 in chromosome 5 or chr5:6,344,259-6,344,259.

Claims

1. Claims1. An in vitro method of determining the DNA methylation level of a set of biomarkers in a subject having thyroid cancer, the method comprising:a. from an isolated or extracted genomic DNA obtained from a biological sample from the subject, determining the methylation level of at least all of the biomarkers selected from the group consisting of cg01378044 (located in TBX18 gene, chr6:85,462, 137-85,462, 137), cg18305394 (located in the 5’ region of PCK1 gene, chr20:56, 134,788-56, 134,788), cg11925561 (located in an intergenic region of chromosome 5, chr5:3, 304, 588-3, 304, 588), and cg14672100 (located in an intergenic region of chromosome 5, chr5:6, 344, 259-6, 344, 259).

2. The method according to claim 1, wherein the methylation level is determined by a quantitative assay selected from the group consisting of bisulfite-based methods, methylation-sensitive restriction enzyme-based methods, microarrays- based methods, next-generation sequencing (NGS) methods, methylation-specific melting curve analysis, mass spectrometry, immunoprecipitation methods, pyrosequencing, digital polymerase chain reaction (digital PCR), quantitative polymerase chain reaction (qPCR) or any combination thereof to measure the DNA methylation level, preferably, bisulfite pyrosequencing and digital PCR.

3. An in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:a. from an isolated or extracted genomic DNA obtained from a biological sample from the subject, determining the DNA methylation level of the set of biomarkers as defined in step a) of any one of claims 1 or 2; andb. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by a classification method that identifies the subject as a subject at risk of developing metastasis based on the biomarkers identified in step (a).

4. An in vitro method for determining the risk of metastasis, preferably of distant metastasis (DM), in a subject having thyroid cancer, the method comprising:a. from an isolated or extracted genomic DNA obtained from a biological sample from the subject having thyroid cancer, determining the DNA methylation level of the set of biomarkers as defined in step a) of claim 1 or 2; andb. identifying the subject as a subject at risk of developing metastasis, preferably distant metastasis (DM), by comparing said methylation level / s with a reference value, wherein a deviation in the methylation level of said at least one biomarker selected from the list shown in a) with respect to said reference value, is indicative that the subject is at risk of metastasis, preferably of distant metastasis (DM).

5. A computer implemented method for determining the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer, wherein the method comprises the steps:a. receiving the methylation levels of at least all of the biomarkers identified in step a) of claim 1;b. integrating the methylation levels of at least all of the biomarkers identified in step a) of claim 1 into a predictive model configured to determine the risk of metastasis, preferably of distant metastasis (DM) in a subject having thyroid cancer; andc. generating with the predictive model of step (b) an output indicative of whether the subject is at risk of having metastasis, preferably of distant metastasis (DM).

6. A computer program comprising processor readable instructions which, when the program is executed by a computer, cause the computer to carry out steps of the method of claim 5.

7. The method according to any one of the preceding claims, wherein the biological sample is a sample containing DNA derived from thyroid tumor.

8. An in vitro method for classifying a subject as a subject having thyroid cancer at risk of metastasis, preferably of distant metastasis (DM), the method comprising:a. determining whether the patient is at risk of metastasis, preferably of distant metastasis (DM), following the method of any one of claims 3 to 6; andb. classifying said patient as a subject having thyroid cancer at risk of metastasis, preferably of distant metastasis (DM), if the patient is identified as a subject at risk of metastasis, preferably of distant metastasis (DM), according to step a).

9. A composition suitable for the treatment of a subject having thyroid cancer at risk of metastasis, for use in the treatment of thyroid cancer in a subject classified as having thyroid cancer at risk of metastasis according to step b) of claim 8.

10. The composition for use according to claim 9, wherein the composition comprises one or more therapies selected from the group consisting of radiation therapies, hormonal therapies, targeted therapies and chemotherapies, or any combination thereof.

Images