Comprehensive transcriptomic analysis of Papillary Thyroid Cancer: potential biomarkers associated with tumor progression

Abstract

Background Identification of stage-specific prognostic/predictive biomarkers could lead to the more efficient clinical management in papillary thyroid carcinoma (PTC). The main objective of this study was to characterize the stage-specific deregulations in gene and miRNA expression in PTC and also to identify potential prognostic biomarkers.

Methods 495 RNASeq and 499 miRNASeq PTC samples (stage I-IV) as well as, respectively, 56 and 57 normal samples were retrieved from The Cancer Genome Atlas (TCGA). DESeq2 was used to identify deregulation of genes and miRNAs between sequential stages of PTC. To identify the minority of patients who progress from stage I to higher stages, we additionally performed a clustering analysis focused on the stage I RNASeq data. Moreover a validation study was done on an independent PTC RNASeq study.

Results Large amount of heterogeneity in gene expression patterns was observed in stage I PTC patients. LTF and PLA2R1 were identified as two promising biomarkers that exhibited down-regulation in a small subgroup of stage I (both in TCGA and in the validation data set) and in the majority of stage IV patients (in TCGA data set). hsa-miR-205, hsa-miR-509-2, hsa-miR-514-1 and hsa-miR-514-2 were also identified as up-regulated miRNAs in both PTC patients with stage I and stage III.

Conclusion Common gene expression alterations in early and advanced stages of PTC, could be used for individualized risk stratification and personalized treatment approach.

Introduction

The incidence of thyroid cancer has been increased worldwide during the past decades. According to Globocan 2018, thyroid cancer is in ninth place for incidence with 567,000 cases worldwide. It is 3 times more prevalent in women with the incidence rate of 10.2 per 100,000 [1]. Four major types of thyroid cancer, including Papillary (PTC), Follicular (FTC), Medullary (MTC) and Anaplastic Thyroid Carcinoma (ATC) have already been characterized, among which PTC is the most frequent histotype, comprising 85-90 % of all thyroid cancer cases. Despite good prognosis, lymph node or distant metastasis is observed in approximately 10% of PTC patients. Although it is still controversial, it has been suggested that there are specific PTC subtypes related to the risk of loco regional recurrence or distant metastasis. [2].

Thyroid cancer staging provides prognostic information related to disease surveillance, risk stratification and therapeutic strategies. Accurate initial staging requires detailed data obtained from preoperative work up, intra operative findings and postoperative follow-up. Moreover risk stratification based on clinical and pathology risk factors or molecular profiling can be efficiently used to guide follow-up management decisions. In addition identification of at-risk patients is fundamental to suggest additional therapies or select a more conservative management approach.

According to the American Joint Committee on Cancer (AJCC) TNM staging system, based on the size of tumor (T), spread to nearby lymph nodes (N), and metastasis to distant sites (M), four stages for thyroid cancer (I-IV) can be specified. Moreover the staging system uses a combination of risk factors such as age at diagnosis, size of the primary tumor, specific tumor histology, direct extension of the tumor to the outside of thyroid gland, loco-regional metastases, and/or distant metastases to stratify patients into various risk categories with differing risks of death from thyroid cancer. PTC has already been subjected to the large-scale genetic characterization studies through which the genetic alterations associated with its tumorigenesis were partly uncovered [3, 31,32]. It has been indicated that activation of mitogen-activated protein kinase (MAPK) and phosphoinositide 3-kinase AKT (PI3K-AKT) pathways through gene fusions or somatic mutations could be considered as primary alterations in PTC [3]. Recently Tengs and colleagues, in a comprehensive transcriptomics study, characterized a signature of deregulated genes and lncRNAs associated with extra-thyroidal extension. Moreover they suggested several biomarkers with the potential use in the early detection of thyroid cancer [4].

In the present study we aimed to investigate the transcriptomic alterations between different PTC stages at both RNA and miRNA levels.  to identify gene specific signatures of different stage transitions. Moreover, since most PTC patients do not progress to the higher stages, we were interested to identify the subgroup of PTC patients in stage I with the potential to progress into the higher stages. Having this knowledge about different PTC high-risk subgroups could help to stratify patients based on the molecular signatures of tumor. This may lead to the adoption of a personalized approach for treating PTC patients and to reduce the financial burden of unnecessary treatment modalities.

Methods

RNASeq and miRNASeq open access data on Thyroid Cancer Adenocarcinoma (PTC type) was retrieved from TCGA (http://portal.gdc.cancer.gov/). Tumor samples were categorized into four stages (“stage I-IV”) according to the clinical data provided by TCGA (Table 1 and 2).

Differential Expression Analysis on RNASeq and miRNASeq data

495 Illumina Genome Analyzer RNASeq row counts (htseq-count) files of PTC in which the expression value of 60,483 transcripts had been measured, and 499 Illumina HighSeq raw counts covering 1,046 miRNAs, beside respectively 56 and 57 RNASeq and miRNASeq htseq-count files of normal thyroid tissues were included in this study. . Figure 1 shows the overall work flow of this study. Normalization steps and differential expression analysis between different stages, including “normal to stage I”, “stage I to II”, “stage II to III” and “stage III to IV” was performed separately using DeSeq2 package [25] in R 3.5.0. Transcripts and miRNAs with row counts less than 10 were removed from all samples before normalization by DeSeq2. The list of the resulted differential expressed genes (DEGs) and miRNAs was then filtered for ≤ 0.05 adjusted p-value and |log₂ FC | ≥ 1.

In order to validate the final results of our RNASeq analysis, we also analyzed another independent PTC RNASeq data set (bio project PRJEB11591) which was the only available data set including 77 PTC (49 stage I, 18 stage III and 2 stage IV ) beside 31 adjacent normal tissue samples. Trimming of paired end Fastq files was performed using Trimmomatic algorithm. Trimmed Fastq files were then mapped on human genome (homosapiens) (b37): hg19 using BWA algorithm. The resulted bam files were then introduced in to htseq-count algorithm. DeSeq2 program was then used to get the list of DEGs between normal vs. stage I and also between stage I vs. the integration of stage III and stage IV htseq-count files, due to the insufficient number of stage IV samples. Similarly, as a result of the small number of stage IV samples, differential expression analysis could not be carried out between stage III and stage IV.

Contrasting and comparing deregulated genes and miRNAs across individual tumor stages

In order to identify differences and similarities in the expression patterns of deregulated genes and miRNAs between different stages (including comparison of changes in direction, i.e., up to down-regulated, or down to up-regulated), the lists of the resulting differentially expressed genes and miRNAs were contrasted and compared between the results of the sequential stage comparisons.

To identify interactions between deregulated miRNAs and their corresponding targets, the lists of deregulated miRNAs in four group studies were separately introduced into miRDB.V6 (http://mirdb.org) and their corresponding target genes were identified. These target genes were then searched in the list of deregulated genes with the opposite pattern of expression (for example target genes of down-regulated miRNAs in each group were searched in the list of up-regulated genes in that group).

“Stage I” sample clustering of RNASeq and miRNASeq data

We clustered all genes in “stage I” samples (278 RNASeq and 282 miRNASeq data sets) normalized htseq-count files,  using hierarchical clustering by using “Euclidean distance” measure in GENE-E program([email protected]).Visually separated clusters of patients with distances ≤ 0.4 were then identified. The expressional changes of each cluster of stage I tumor samples were compared with normal samples and the lists of deregulated genes and miRNAs in each ”stage I” clusters were identified using DESeq2. Finally the lists of down and up-regulated genes and miRNAs in each cluster were compared with the lists of deregulated genes in different stage transitions (from normal to stage IV).

Similarly the clustering of RNASeq stage I samples was carried out for 49 stage I samples of the validation data set.

Gene enrichment analysis of identified deregulated genes

To find deregulated pathways in each stage comparisons we performed pathway enrichment using Enrichr [26] and pathways with Benjamini-Hochberg adjusted p-value ≤ 0.05 was reported.

Results

mRNAs and miRNAs deregulations

The list of deregulated genes and miRNAs in different sequential stage progressions are respectively represented in Table 3-4 and Table 5. Deregulated genes and miRNAs were ranked based on their log₂ FC of expression values. The results of the RNASeq validation study (normal to stage I and stage I to the integration of stage III-IV) showed the statistical significant results (adjusted p-value ≤ 0.05) for only normal to stage I comparison.

mRNA

Figure 2 and 3 show KEGG signaling pathways identified for up and down-regulated genes, respectively. Simultaneous down-regulation of RAC2, VAV1 and DOCK2 genes was observed in stage I to II progression. These proteins are involved in leukocyte transendothelial migration pathway [33]. ITGAL and ITGB2 were also two down-regulated integrins [33, 34] (in stage I to II) that are present at the leukocyte membrane and involved in leukocyte binding to the endothelial membrane which in turn leads to the direct activation of VAV1 and indirect activation of Rac2. Concerted up-regulation of several extra cellular matrix (ECM) proteins, including THBS1, THBS2, COL6A3, SPP1, TNC, FN1, LAMB3 and COL1A2, along with ITGA11 (an ECM associated receptor) was observed in stage II to III transition. The up-regulation of these proteins results in the considerable activation of PI3K-AKT pathway and therefore leads to the continuous cell proliferation [33]. Moreover the activation of TGF-beta signaling pathway which contributes to the epithelial-mesenchimal transition (EMT), has previously been indicated following the up-regulation of THBS1( up-regulated in stage II to III progression) [35].

The up-regulation of three MHC class II members, including HLA-DOA, HLA-DQA1 and HLADRB1, was detected in stage II to III progression. The higher frequency of HLA-DRB1 in PTC patients compared to controls has been already demonstrated [36]. Similar results were demonstrated by JO and colleagues. They also indicated that the expression of MHC class II antigens (HLADQ and HLA-DR) is inversely correlated with the recurrence of PTC [37].

Comparing the list of up and down-regulated genes between pairs of continuous stages resulted in the identification of 3 down-regulated genes that were also detected as down-regulated in higher stages. These include FBLN1, identified as down-regulated in the “normal to stage I”, and also in “stage I to II” progression. Similarly, LTF and PLA2R1 were identified as down-regulated in both “normal to stage I” and “stage III to IV” progressions. No common up-regulated gene was identified in more than one stage comparison.

In order to identify any gene expression inversion between two continuous stages (e.g. down/up-regulated in “normal to stage I”, but up/down-regulated in stage I to II), the list of DEGs between each pair of continuous stages was compared with its immediate next stages pair. This led to the identification of a few genes with inverse pattern of expression (Table 6). Notably, 28 down-regulated genes in “stage I to II” progression were identified as up-regulated in “stage II to III”. The corresponding KEGG signaling pathways for these 28 genes has been shown in Figure 4. Several numbers of these genes have a direct or indirect role in the activation of apoptosis. These include HLA-DOA, HLA-DQA, BIRC3, COL1A1 and IL2RG. Down-regulation of these genes in stage II is expected in the primary stages of tumor development. However their over-expression, which leads to the increase of apoptotic activity, has also been previously confirmed in some malignancies including breast, endometrial and thyroid cancers [38]. This increase in the apoptotic activity has previously been shown to be due to the presence of hypoxic conditions in the core regions of solid tumors and indicating the large area occupied by tumor cells [24]. Similarly, changing the pattern of expression from down to up-regulated was observed for RAC2 which along with DOCK2 and VAV1 (identified as down-regulated in stage I to II progression), are associated with the regulation of actin cytoskeleton, production of reactive oxygen species (ROS) and metastasis [33].

miRNA

There is large discrepancy among the results of studies on miRNA expression in thyroid cancer [20]. However the deregulation of several identified miRNAs in this study has previously been well confirmed in PTC (bold fonts in Table 5). Among previously reported deregulated miRNAs, the up-regulation of hsa-miR-146b and hsa-miR-222 and the down-regulation of hsa-miR-144 have been also reported in the blood of PTC patients [21, 22]. A similar comparison to that of mRNAs was also carried out for differential expressed miRNAs to identify further deregulation of their expression in higher stages. hsa-miR-205,  hsa-miR-509-2, hsa-miR-514-1, and hsa-miR-514-2 were four miRNAs whose up-regulation was detected in both ”normal to stage I” and in “stage II to III” progressions. Similarly, hsa-miR-9-1 and hsa-miR-9-2 were identified as down-regulated in “normal to stage I” and “stage I to II” progressions. Moreover, gene expression pattern inversion was observed for hsa-miR-1247, which was detected as down-regulated in “stage I to II”, and as up-regulated in “stage II to III” progression.

Results of stage I tumor sample clustering

Clustering of stage I subjects revealed 6 clusters for RNASeq, and 5 clusters for miRNASeq (Table 7) (Figures 5 and 6). Stage I RNAseq sample clustering revealed much more heterogeneity than that of miRNASeq. Therefore 6 identified clusters (totally comprising 265 out of 278 tumor samples, after removing outliers) of stage I RNASeq samples were subjected to further analyses. Comparing the lists of down and up-regulated genes in each “stage I” cluster (Supplementary table 1) with the lists of deregulated genes in higher stages led to the identification of few deregulated genes (Table 8).

LTF which its down-regulation was also identified in “stage III to IV” progression, belongs to the cluster 3 of “stage I”, which includes 102 samples. Since we did not expect tumor progression to the highest stage (stage IV) in a large number of PTC patients, we further sub-clustered these 102 samples in order to find a specific sub-cluster with the most contribution in the observed down-regulation of LTF. A total of 7 sub-clusters were identified in “cluster 3” of “stage I”. Running DeSeq2 showed that in only one sub-cluster (containing 21 samples), LTF was significantly down-regulated (log₂ FC =-3). This level of down-regulation (8 fold reduction in the expression) appeared to influence LTF average expression in the whole cluster 3 (Table 8).

No deregulated genes were detected in “cluster 4” of stage I compared with the higher stages. Therefore it seems that patients in clusters “1, 2, 5 and 6” of stage I, as well as 21 patients in “cluster 3” (totally comprising 71 subjects out of 278 stage I tumor samples) which have shown further deregulation in higher stages, would probably be in the higher risk for tumor progression. No statistically significant difference between the potential high and low risk (71 vs. 207) groups of stage I patients was detected based on their age, age of onset, race and sex.

The changes in the expressional pattern of two genes i.e. LTF and PLA2R1 was also checked in different stage I clusters’ of the validation data set. In total 3 clusters of samples (including C1=12, C2=10 and C3=6 samples) were identified in stage I. While the expression of LTF and PLA2R1 showed mild down-regulation in C1 and C2, their significant down-regulation was observed in C3 which includes only 6 samples (log₂ FC = -3.78 and -3.57 for LTF and PLA2R1 respectively, equals to approximately 13.5 fold expressional reduction). This indicates almost similar percentage of stage I patients with significant down-regulation of LTF and PLA2R1 (for LTF, 21 out of 262 samples (8%) in TCGA data set and 6 out of 49 (12%) in the validation data set and for PLA2R1, 42 out of 262 samples (16%) in TCGA data set and 6 out of 49 (12%) in the validation data set).

Results of comparing down and up-regulated miRNAs in each miRNAs clusters with the list of deregulated miRNAs obtained from stage comparisons was resulted in to the identification of hsa-miR-1247 whose down-regulation was observed in cluster 4 of “stage I”, and also in “stage I to II” progression. However it was also detected as up-regulated in “stage II to III” progression.

Discussion

Despite good prognosis of PTC, it has been indicated that a sub-group of PTC patients show an aggressive phenotype characterized by local recurrence and/or distant metastasis [13]. Characterizing DNA level molecular alterations of PTC including mutations and structural variations has already been the subject of several studies [2, 30]. Recently, a few transcriptomic investigations have also been conducted to generate a comprehensive view of PTC molecular alterations [4, 39]. Results of these studies have suggested several biomarkers that may have application in the early detection of PTC or as novel therapeutic targets. However no systematic study has been conducted to unravel the stage specific transcription alterations occurring in PTC. The identified stage-specific transcriptional alterations in our results can significantly improve our understanding about the gene expression signature of different PTC stages and suggests potential biomarkers which could potentially be implemented to classify subgroups of patients who might progress to the more advanced stages.

Biomarkers with prognostic potential in PTC

The deregulation of a total of 12 genes in different sub-clusters of stage I were further identified in higher stages (Table 7). As we expected, for the majority of stage I PTC samples (207), no further deregulation was observed in the higher stages.

LTF and PLA2R1 were two down-regulated tumor suppressors in “stage I” that, interestingly, exhibited further down-regulation in “stage III to IV” progression. The association of down-regulation of LTF with the poor prognosis of PTC has recently been indicated [5].The role of PLA2R1 in the activation of signaling pathways related to apoptosis and senescence and its down-regulation in some cancers including thyroid, renal, mammary and leukemia has been also previously established [6]. As stated earlier, the significant down-regulation of LTF and PLA2R1 was also confirmed in respectively 12% and 16% of stage I samples of the validation data set. However as a result of the small number of stage IV samples (only 2 samples), we could not check their expressional changes in stage III to IV progression. Nevertheless it seems that these two genes can be considered as the potential biomarkers associated with the poor prognosis in PTC.

TRIB3 and KIF5C are the other two down-regulated genes in stage I which were also detected as down-regulated in stage II to III progression. TRIB3 is a tumor suppressor which its loss has been shown to lead to the enhanced phosphorylation of AKT and in turn to the hyperphosphorylation of FOXO3 which results to its inactivation [7]. There are contradictory reports about  KIF5C pattern of deregulation in thyroid cancer (both the down and up-regulation of KIF5C in thyroid tumors has already been reported)[8, 9].

There is still no report about the deregulation of SNX20 in thyroid cancer. However in a very recent case report study, somatic mutations of this gene have been indicated in both tissue DNA and in circulating cell free DNA (cfDNA) of a follicular thyroid carcinoma patient with lung and bone metastasis [10].

The further up-regulation of all seven up-regulated genes in different sub-clusters of stage I were identified in stage II to III progression (Table 8). Among these, CST7 and FN1 were two genes that have previously been suggested as biomarkers associated with poor prognosis of PTC [4,11]. Similarly the over-expression of TNFSF13B in tumor cells has been shown to be related to the poor prognosis in non small cell lung cancer [12].The up-regulation of SCEL in PTC has previously been demonstrated by Huang and colleagues [14]. The genetic polymorphism of IL1RN has been previously associated with the poor prognosis in cancer [15]. Moreover IL1RN has been already suggested as a biomarker involved in tumor differentiation of thyroid carcinoma [13].

It has been confirmed by several studies that the hypermethylation of Runx3 is in association with poor prognosis in thyroid cancer [17, 18, 26]. However we identified RUNX3 as a up-regulated gene in both a very small sub-cluster of “stage I” (including 3 samples) and also in stage “I to II” progression. Interestingly the results of a recent comprehensive meta-analysis conducted on three different cancers (including Breast, Colon and Lung), provide further interesting evidence about the direct relationships between the level of hypermethylation of genes and the levels of genes expression. They indicated that except for the genes that are lowly expressed in normal tissues, the hyper or hypomethylation of genes in cancer tissue does not necessarily result in respectively their down or up-regulation [19].

No significant difference was identified between the age, age of onset, race and sex of potential 71 low risk and 207 high risk stage I patients. This indicates the importance of using transcriptional biomarkers for classifying PTC patients and select efficient personalized treatment.

Prognostic miRNAs for PTC tumor progression

As stated earlier, the results of miRNA expression studies vary significantly due to the several factors including the heterogeneity of tumor samples, the variations in sample size and in the quantitative method implemented for expression measurements. The problem becomes more complex considering the fact that there is a direct correlation between the deregulation of some miRNAs with some oncogenic mutations. For instance, the up-regulation of hsa-miR-146b, hsa-miR-221, hsa-miR-222, hsa-miR-21-5p, hsa-miR-31-5p, hsa-miR-34a-5p, hsa-miR-181b-5p and hsa-miR-551b-3p have been confirmed in BRAF V600E samples [23].

In spite of the large number of analyzed tumor samples in the present study, due to the imbalance between the number of tumor vs. normal samples, the results of this study are not comparable to other studies performed on miRNA expression alterations in PTC. However as it is shown in Table 5, the deregulation of a considerable number of identified miRNAs have already been confirmed by several studies. As stated in the results section, the down-regulation of hsa-miR-205, hsa-miR-509-2, hsa-miR-514-1, and hsa-miR-514-2 were detected in both stage I and stage III and therefore after further confirmation they can be implemented as new predictors of tumor progression. There are contradictory reports about the expression pattern of hsa-miR-1247 (down-regulated in stage I and up-regulated in stage III). Down-regulation of hsa-miR-1247 has previously been reported in thyroid tumors with follicular pattern of growth [29]. Similarly the down-regulation of hsa-miR-1247 has been confirmed in the lymph node metastasis in breast cancer [27]. However the association between the up-regulation of this miRNA with lung metastasis has been shown in liver cancer. It has been indicated that tumor derived exosomal hsa-miR-1247 affects cancer associated fibroblasts and activates lung metastasis [28].

Expression measurements of other identified deregulated miRNAs with no previous reports on their association with PTC, in a large number of PTC patients, will reveal if they can be truly considered as the novel prognostic biomarkers in PTC. The result of searching for the targets of deregulated miRNAs among the list of deregulated mRNAs led to the identification of only one interaction between HSPA1B (down-regulated in normal to stage I) and hsa-miR-34a (up-regulated in normal to stage I). This is however in accordance with studies demonstrating that the most alterations in mRNA expression are short term and independent of miRNAs [23].

We believe that the strength of this study lies in its stage-wise approach, which to the best of our knowledge, has not already been performed in PTC. Moreover the large sample size used in this study increases the confidence of our results. The main limitations of this study include the small number of normal samples compared to the large number of tumor samples and also the lack of complete demographic data for the PTC patients included in this study.

Declarations

Acknowledgement:

The authors gratefully acknowledge the supports provided by Iran University of Medical Sciences, (Grant number 33444).

Disclosure statement:

Authors declare no competing interests.

Author Contributions:

M.K designed and supervised the project. N.H performed the data analyses and wrote the manuscript. K.B, M.K and M.H commented on the manuscript. T.M reviewed the paper.

References

[1] Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians, 2018;68(6):394-424.

[2] Katoh H, Yamashita K, Enomoto T, Watanabe M. Classification and general considerations of thyroid cancer., Ann Clin Pathol. 2015;3(1):1045.

[3] Liu R, Zhang T, Zhu G, Xing M. Regulation of mutant TERT by BRAF V600E/MAP kinase pathway through FOS/GABP in human cancer., Nature communications. 2018;9(1):579.

[4] Teng H, Mao F, Liang J, Xue M, Wei W, Li X, et al. Transcriptomic signature associated with carcinogenesis and aggressiveness of papillary thyroid carcinoma., Theranostics. 2018;8(16):4345.

[5] Brennan K, Holsinger C, Dosiou C, Sunwoo JB, Akatsu H, Haile R, et al. Development of prognostic signatures for intermediate-risk papillary thyroid cancer. BMC cancer. 2016;16(1):736.

[6] Menschikowski M, Hagelgans A, Nacke B, Jandeck C, Sukocheva O, Siegert G. Epigenetic control of phospholipase A 2 receptor expression in mammary cancer cells. BMC cancer. 2015;15(1):971.

[7] Salazar M, Lorente M, García-Taboada E, Gomez EP, Dávila D, Zúñiga-García P, et al. Loss of Tribbles pseudokinase-3 promotes Akt-driven tumorigenesis via FOXO inactivation. Cell death and differentiation. 2015;22(1):131.

[8] Lu C, Mishra A, Zhu YJ, Meltzer P, Cheng S-y. Genomic profiling of genes contributing to metastasis in a mouse model of thyroid follicular carcinoma. American journal of cancer research. 2011;1(1):1.

[9] Rodriguez MdCP, Freire MdCC, Cameselle-Teijeiro JM, Puente FD, Figueroa AV. Methods for diagnosing follicular thyroid cancer. Google Patents; 2012.

[10] Song J, Yang Z. Case report: Whole exome sequencing of circulating cell-free tumor DNA in a follicular thyroid carcinoma patient with lung and bone metastases. Journal of circulating biomarkers. 2018;7:1849454418763725.

 [11] Cerutti JM, Oler G, Michaluart P, Delcelo R, Beaty RM, Shoemaker J, et al. Molecular profiling of matched samples identifies biomarkers of papillary thyroid carcinoma lymph node metastasis. Cancer research. 2007;67(16):7885-92.

[12] Qian Z, Qingshan C, Chun J, Huijun Z, Feng L, Qiang W, et al. High Expression of TNFSF13 in Tumor Cells and Fibroblasts Is Associated With Poor Prognosis in Non–Small Cell Lung Cancer. American journal of clinical pathology. 2014;141(2):226-33.

[13] Blanco CG, Matute EM, de Leiva Hidalgo A. Molecular biomarkers involved in the tumor dedifferentiation process of thyroid carcinoma of epithelial origin: perspectives. Endocrinología y Nutrición (English Edition). 2012;59(7):452-8.

[14] Huang Y, Prasad M, Lemon WJ, Hampel H, Wright FA, Kornacker K, et al. Gene expression in papillary thyroid carcinoma reveals highly consistent profiles. Proceedings of the National Academy of Sciences. 2001;98(26):15044-9.

[15] Niedźwiecki S, Stępień T, Kuzdak K, Stępień H, Krupiński R, Seehofer D, et al. Serum levels of interleukin-1 receptor antagonist (IL-1ra) in thyroid cancer patients. Langenbeck’s archives of surgery. 2008;393(3):275-80.

[16] Chen F, Liu X, Bai J, Pei D, Zheng J. The emerging role of RUNX3 in cancer metastasis. Oncology reports. 2016;35(3):1227-36.

[17] Joung KH, Shong M.Epigenetic regulation of RUNX3 in thyroid carcinoma. The Korean journal of internal medicine. 2012;27(4):391.

[18] dos Reis M, Beltrami C, Kowaslki L, Rogatto S. Epigenetic alterations in well-differentiated thyroid cancer. J Clin Epigenet. 2016;1:1.

[19] Moarii M, Boeva V, Vert J-P, Reyal F. Changes in correlation between promoter methylation and gene expression in cancer. BMC genomics. 2015;16(1):873.

[20] Celano M, Rosignolo F, Maggisano V, Pecce V, Iannone M, Russo D, et al. MicroRNAs as biomarkers in thyroid carcinoma. International journal of genomics. 2017;2017.

[21] Poulios E, Zografos G, Tseleni-Balafouta S, Linos D. Identification of a Novel, Diagnostic microRNA Signature in Papillary Thyroid Cancer. Journal of Medical & Surgical Pathology. 2018, 3(3): 162.

[22] Yu S, Liu Y, Wang J, Guo Z, Zhang Q, Yu F, et al. Circulating microRNA profiles as potential biomarkers for diagnosis of papillary thyroid carcinoma. The Journal of Clinical Endocrinology & Metabolism. 2012;97(6):2084-92.

 [23] Saiselet M, Pita JM, Augenlicht A, Dom G, Tarabichi M, Fimereli D, et al. miRNA expression and function in thyroid carcinomas: a comparative and critical analysis and a model for other cancers. Oncotarget. 2016;7(32):52475.

[24] Gottfried, Y., et al., Expression of the pro‐apoptotic protein ARTS in astrocytic tumors: Correlation with malignancy grade and survival rate. Cancer: Interdisciplinary International Journal of the American Cancer Society, 2004. 101(11): p. 2614-2621.

[25] Love, M.I., W. Huber, and S. Anders, Moderated estimation of fold change and dispersion for RNAseq data with DESeq2. Genome biology, 2014. 15(12): p. 550.

[26] Chen, E.Y., et al., Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC bioinformatics, 2013. 14(1): p. 128.

[27] Zhang, P., et al., Association of miR‐1247‐5p expression with clinicopathological parameters and prognosis in breast cancer. International journal of experimental pathology, 2018. 99(4): p. 199-205.

[28] Fang, T., et al., Tumor-derived exosomal miR-1247-3p induces cancer-associated fibroblast activation to foster lung metastasis of liver cancer. Nature communications, 2018. 9(1): p. 191.

[29] Mancikova, V., et al., MicroRNA deep-sequencing reveals master regulators of follicular and papillary thyroid tumors. Modern Pathology, 2015. 28(6): p. 748.

[30] Acquaviva, G., et al., Molecular pathology of thyroid tumours of follicular cells: a review of genetic alterations and their clinicopathological relevance. Histopathology,  2018. 72(1): p. 6-31.

[31] Agrawal, N., et al., Integrated genomic characterization of papillary thyroid carcinoma. Cell, 2014.159(3): p. 676-690.

[32] Haraldsdottir, S. and M.H. Shah, New era for treatment in differentiated thyroid cancer. The Lancet, 2014. 384(9940): p. 286-288.

[33] Kanehisa, M., et al., KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic acids research, 2016. 45(D1): p. D353-D361.

[34] UniProt: the universal protein knowledgebase. Nucleic acids research, 2016. 45(D1): p. D158-D169.

[35] Seliger, C., et al., Lactate-modulated induction of THBS-1 activates transforming growth factor (TGF)-beta2 and migration of glioma cells in vitro. PloS one, 2013. 8(11): p. e78935.

[36] Amoli, M.M., et al., HLA-DR association in papillary thyroid carcinoma. Disease markers, 2010. 28(1): p. 49-53.

 [37] Jo, Y.S., et al., Significance of the expression of major histocompatibility complex class II antigen, HLA‐DR and‐DQ, with recurrence of papillary thyroid cancer. International journal of cancer, 2008. 122(4): p. 785-790.

[38] Soini, Y., P. Pääkkö, and V. Lehto, Histopathological evaluation of apoptosis in cancer. The American journal of pathology, 1998. 153(4): p. 1041-1053.

[39] Liyanarachchi, S., et al., Genome-wide expression screening discloses long noncoding RNAs involved in thyroid carcinogenesis. The Journal of Clinical Endocrinology & Metabolism, 2016. 101(11): p. 4005-4013.

Tables

Table 1: Number of samples in stage I-IV as well as normal thyroid tissue samples.

Table 2: Clinical information of samples (Stage I-V).

Table 3:  Up-regulated genes involved in sequential stage progression.

* log₂ FC of expression

Table 4:  Down-regulated genes in PTC continuous stage progressions.

* log₂ FC of expression

Table 5:  Deregulated miRNAs in the different stage progressions.

Known deregulated miRNAs are shown in bolds

* log₂ FC of expression

Table 6:  Lists of DEGs with the opposite pattern of expression between different stages.

Table 7: Different clusters identified by hierarchical clustering in stage-I samples.

Table 8: Deregulated genes in stage I clusters with the same pattern of deregulation in higher stages.

* LTF significant down-regulation (log₂ FC = −3) was observed in 21 out of 102 PTC patients in cluster 3.