Evolution of intra-tumoral heterogeneity across different pathological stages in papillary thyroid carcinoma

Intra-tumor heterogeneity (ITH) results from the continuous accumulation of mutations during disease progression, thus impacting patients’ clinical outcome. How the ITH evolves across papillary thyroid carcinoma (PTC) different tumor stages is lacking. We used the whole-exome sequencing data from The Cancer Genome Atlas Thyroid Cancer (TCGA-THCA) cohort to track the ITH and assessed its relationship with clinical features through different stages of the PTC progression. We further assayed the expression levels of the specific genes in papillary thyroid cancer cell lines compared to an immortalized normal thyroid epithelial cell line by qRT-PCR. We revealed the timing of mutational processes and the dynamics of the temporal acquisition of somatic events during the lifetime of the PTC. ITH significantly influences the PTC patient’s survival rate and, as genetic heterogeneity increases, the prognosis gets worse in advanced tumor stages. ITH also affects the mutational architecture of each clinical stage which is subject to periodic fluctuations. Different mutational processes may cooperate to shape a stage-specific mutational spectrum during the progression from early to advanced tumor stages. Moreover, different evolutionary paths characterize PTC progression across pathological stages due to both mutations recurrently occurring in all stages in hotspot positions and distinct codon changes dominating in different stages. A different expression level of specific genes also exists in different thyroid cancer cell lines. Our findings suggest ITH as a potential unfavorable prognostic factor in PTC and highlight the dynamic changes in different clinical stages of PTC, providing some clues for the precision medicine and suggesting different diagnostic decisions depending on the clinical stages of patients. Finally, complete clear guidelines to define risk stratification of PTC patients are lacking; thus, this work could contribute to defining patients who need more aggressive treatments and, in turn, could reduce the social burden of this cancer.


Background
Tumorigenesis is often initiated by a single mutated neoplastic cell evolving through a series of sequential clonal or subclonal mutational events, thus influencing the course of disease [1,2]. In this way, cell clones diversify, resulting in intratumor heterogeneity (ITH) [3,4]. As

Open Access
Cancer Cell International † Giuliana Salvatore and Monica Franzese contributed equally to this work *Correspondence: ornella.affinito@synlab.it a result, cancer evolves as a mosaic entity composed of a mixture of cells with distinct genetic, phenotypic or behavioral features within the same tumor [5]. Changes in selection pressure within tumor microenvironment may also enhance clonal diversity or allow the positive selection of more aggressive clones with subclonal mutations that increase the tumor's fitness to its environment and determine a more aggressive evolution of diseases [6]. ITH may provide detailed records of tumor clonal evolution history and evolutionary dynamics of mutations, chronologically accumulated and selected during the lifetime of a tumor [7][8][9]. Deciphering and tracking the clonal evolution pattern could provide valuable information regarding crucial genetic events in tumorigenesis and progression, thus greatly benefiting therapy selection and prognosis management. Indeed, ITH represents a significant challenge in the implementation of precision medicine because it is responsible for the progression from an early tumor stage to a more aggressive cancer. Moreover, ITH is one of the major causes of poor prognosis, treatment failure and drug resistance [8,[10][11][12]. Some evidence suggests that ITH influences the clinical outcome in several cancers, such as chronic lymphocytic leukemia [6], head and neck cancer [10,11], colorectal [12] and lung adenocarcinoma [13]. Moreover, ITH has the potential to be a valuable predictor for clinical outcomes [14] and may serve as a clinically useful biomarker in the development of personalized therapies and clinical outcomes [15,16]. Thus, exploring intratumor heterogeneity and tumor evolution is of great clinical importance because, depending on the tumor evolution, different clinical approaches may exist in terms of diagnosis, prognosis, and treatment of the individuals [17].
Despite its potential clinical relevance, ITH was poorly investigated in papillary thyroid carcinoma (PTC). PTC is the most frequent endocrine tumor accounting for almost 80% of thyroid cancer cases. Somatic mutations of genes involved in the mitogen-activated protein kinase (MAPK) signaling pathway, including point mutations in BRAF, RAS genes [18][19][20][21] and RET/PTC rearrangement [22]. Although PTC has in general a good prognosis with a 5-year survival of over 95% [23], in about 85% of case subjects, a small fraction (10-20%) shows higher aggressiveness with either local or distant relapse [24,25]. A large Cancer Genome Atlas (TCGA) Research Network study [26], based on genetic, epigenetic and transcriptomic analysis of almost 500 PTCs, demonstrated substantial inter-tumor heterogeneity in PTCs concerning their molecular alterations. Moreover, they found that driver mutations in genes such as BRAF, NRAS, HRAS, KRAS, and EIF1AX were present in most PTC cases and were largely clonal. On the other side, as regards the intra-tumoral heterogeneity, most studies focused on the most prevalent genetic alterations in PTC, the BRAF V600E mutation and RET rearrangements [27][28][29][30]. However, while each cancer cell may contain a driver mutation, the PTC cells are still proliferating and can diversify through the gain of subclonal genetic alterations that can be positively selected and might impact prognosis and clinical outcome [31]. Moreover, these non-driver alterations may provide insight into a tumor's history and help to identify mutational processes occurring during tumorigenesis [32]. Some evidence supports the occurrence of ITH in PTC [33], but the relevance of subclonality in PTC is still debated [34]. However, a comprehensive analysis taking into account the changes and the impact of ITH through different pathological stages of the PTC, along with the timing of mutational processes and the dynamics of the temporal acquisition of somatic events, is lacking.
In the current study, we analyzed the whole-exome sequencing data from The Cancer Genome Atlas Thyroid Cancer (TCGA-THCA) cohort to efficiently and comprehensively evaluate somatic variants across clinical stages of PTC. We tracked the intra-tumoral genetic heterogeneity by mutant-allele tumor heterogeneity (MATH) algorithm [35] and assessed its relationship with clinical features through different stages associated with the PTC progression. We also revealed how mutational processes vary over time during tumor evolution and identified several evolutionary patterns depending on the behavior of the mutated genes. We further investigated the relationship between the mutational status and RNA expression level. These findings highlighted the dynamic changes of oncogenesis in different clinical stages of PTC, providing some clues for the development of precision medicine and the improvement of diagnostic strategies in PTC patients.
There were 487 patients with mutation data and 507 with available clinical and survival data. Data were filtered to exclude patients without mutation or clinical information. We also excluded patients without available information on tumor stage and retained patients with the following histological types: classical, follicular and tall cell variant papillary thyroid carcinoma. The final dataset was composed of 474 patients consisting of 265 stage 1, 50 stage 2, 108 stage 3 and 51 stage 4.

Somatic mutation analysis and measurement of heterogeneity
We used the R/Bioconductor package Maftools [38] to efficiently and comprehensively analyze somatic variants in PTC. Genetic intra-tumor heterogeneity (ITH) was assessed by mutant-allele tumor heterogeneity (MATH) score [35]. To estimate MATH scores, the InferHeterogeneity function of the Maftools package was used, which infer clonality by clustering variants with similar allele frequencies and removing patients with less than 3 variants. The MATH score for each PTC patient was calculated according to the method described by Mroz et al. [10]. The MATH score is defined as the percentage ratio of the median absolute deviation (MAD) and the median of the distribution of MAFs among the tumor's mutated genomic loci: MATH = 100*MAD/median. For each tumor stage, patients were categorized according to tumor heterogeneity into high-and low-MATH groups by the median of the MATH scores. Highand low-MATH groups were then analyzed with respect to clinicopathologic features and Progression-free survival (PFS) data.

Mutation signature analysis
Somatic mutations of each sample were classified according to base substitution (C > A, C > G, C > T, T > A, T > C, T > G) and immediately 5′/3′ base information. We identified mutational context based on the human reference genome hg38 through BSgenome.Hsapiens.UCSC.hg38 R package. The resulting triplet SNV spectra were analyzed for contributions of known mutational signatures in Catalogue of Somatic Mutations in Cancer (COSMIC; https:// cancer. sanger. ac. uk/ cosmic/ signa tures).
Mutational signatures were predicted using the deconstruct-Sigs R package [39]. This tool evaluates the contribution of the signatures reported in COSMIC (used as reference mutational signatures) to the mutational profile of the somatic SNVs in each tumor stage. Mutations not included in previously identified signatures were classified as "unknown".

Expression level analysis using web resources
The expression level of the common mutated genes (JMJD1C, MALAT1, MUC16, PDZD2, PKHD1, RYR1, SLA and TTN) was analyzed across thyroid tumor tissues and normal thyroid samples and across thyroid cancer stages deposited in TCGA project using UALCAN data portal (http:// ualcan. path. uab. edu) [40].

RNA extraction and qRT-PCR
Total RNA was extracted from an immortalized normal thyroid epithelial cell line (Nthy-ori 3-1) and from three papillary thyroid cancer cell lines (K1, TPC-1 and BCPAP) using TRIzol reagent (Catalog number 15596026) purchased from Thermo Fisher Scientific (Waltham, USA).
After extraction, RNA was quantified using the Nan-oDrop spectrophotometer (Thermo Fisher Scientific). Next, cDNA was synthesized using 1 µg of total RNA from each sample using QuantiTect Reverse Transcription (Catalog number 205311) purchased from Qiagen (Hilden, Germania) according to the manufacturer's instructions.
The relative expression level (expressed as fold change) was determined by applying the formula 2 −ΔΔCt [41] where the Ct-value of each gene was technically normalized using β-ACTIN. The experiment was carried out three independent times in triplicate.

DNA sequencing
Genomic DNA was purified from cultured thyroid cells (Nthy-ori 3-1, K1, TPC-1 and BCPAP) using QIAamp DNA mini kit (Catalog number, 51304) purchased from Qiagen (Hilden, Germany) and quantified with the Nanodrop spectrophotometer (Thermo Fisher Scientific, Massachusetts, USA). Next, we designed specific primers able to amplify the exonic region of genes containing the high impact mutations. PCR was carried out on Applied Biosystems 2720 Thermal Cycler (Thermo Fisher Scientific) using Taq DNA polymerase (Catalog Number EP0401, Thermo Fisher Scientific) with the following thermal cycles: denaturation at 95 °C (for 3 min), 40 cycles at 95 °C (for 30 s), 60 °C (for 30 s), 72 °C (for 1 min), followed by 72 °C (for 5 min).
The genes, the exons harboring the high impact mutations found in PTC patients, and the primer sequences used are listed in Additional file 1: Table S1.
The amplicons sequencing was performed as follows: 8 μl of primer forward (concentrated 2 μM) and 20 μl of DNA (concentrated 16 ng/μl). Sanger electropherograms were analyzed to evaluate the presence/absence of the indicated mutations.

Statistical analysis
The normal distribution assumption for continuous variables was assessed by Shapiro-Wilk test. For comparisons between two groups, a two-tailed t-test for independent samples (for normally distributed data) or a Wilcoxon rank sum test (for not normally distributed data) was used. Comparisons among more than two groups were made using the Kruskal-Wallis test. p-values were also corrected for multiple testing by the Bonferroni method. Progression-free survival (PFS) data were used as endpoints for survival analysis. Survival rates were analyzed with the Kaplan-Meier method and the statistical relevance of the differences between survival curves was assessed by log-rank test. Univariate analysis Cox proportional hazards regression models were applied to evaluate the prognostic value of the MATH score and to analyze the effects of clinical confounding factors on survival. All confidence intervals (CIs) were stated at the 95% confidence level. Statistical significance was set at p-value ≤ 0.05. All analyses were performed using R 3.6.0 (https:// www.r-proje ct. org/).

Tumor stage-related survival analysis
We analyzed exome sequencing data from 474 patients with PTC obtained from the TCGA data portal. Clinicopathological information of the patients included in the current analysis is summarized in Table 1 and Additional file 2: Table S2.
A Kaplan-Meier survival analysis was conducted using PFS data to analyze the survival probability in different tumor stages (Fig. 1). The results showed a statistically significant difference in survival probability depending on the tumor stage (Log-rank test; p-value < 0.0001), with a worse survival odd for patients in stage 4. In addition, we performed an univariate Cox proportional hazards regression to evaluate the association between different tumor stages and prognosis ( Table 2). The analysis showed that stages 3 and 4 were significantly associated with a poor prognosis (p-value = 0.0073 and p-value = 1.9e-05, respectively). Instead, patients in stage 2 also had a poor survival but not significant (p-value = 0.61).

Potential role of the MATH score in PTC for tumor stage-related survival
We evaluated whether intratumor heterogeneity (ITH) could explain the different prognosis of tumor stages. ITH can be assessed using the mutant allele tumor heterogeneity (MATH) score. As reported in Fig. 2A, the MATH scores distribution showed a broad spectrum of values (from 0.49 to 89.01 with a mean ± standard deviation (SD) of 28.64 ± 16.33), suggesting that PTC patients exhibit remarkable intra-tumoral heterogeneity. MATH score distribution according to tumor stage highlighted a statistically significant difference (Kruskal-Wallis rank sum test; p-value = 0.033) (Fig. 2B). According to the  (Fig. 2C, D). Next, we investigated the prognostic significance of MATH for each MATH group through a Kaplan-Meier survival analysis (Fig. 2E, F) and an univariate Cox proportional hazards regression (Table 3). Kaplan-Meier survival analysis showed that the MATH values were significantly associated with tumor stage (Fig. 2E, F). Specifically, we found that patients in advanced stages (stage 3 and stage 4) exhibited shorter PFS rates than those in early stages (stage 1 and stage 2) in both MATH groups (Log-rank test; p-value = 0.015 for the low-MATH group; p-value = 0.0017 for the high-MATH group) (Fig. 2E, F). The univariate Cox proportional hazards regression analysis underlined the role of MATH score as potential risk factor associated with PTC aggressiveness. Specifically, low MATH score showed a statistical significance in stage 3 (p-value = 0.021) and stage 4 (p-value = 0.007), whereas high MATH score showed a statistical significance in stage 4 (p-value = 0.0005) ( Table 3).

Relationship between MATH and clinical features
We also investigated whether the tumor heterogeneity was associated with clinical features of PTC patients (Table 4, Additional file 3: Table S3 and Additional file 4: Table S4). We performed a Kaplan-Meier analysis followed by a univariate Cox regression model to establish the impact of different clinical features on survival and prognosis in low and high MATH groups. In the low-MATH group for stage 4, M stage was significantly associated with survival (Log-rank test; p-value = 0.02) with a poorer prognosis (p-value = 0.04) for M1 stage (p-value = 0.04) ( Table 4). Moreover, the primary neoplasm focus type and neoplasm length clinical features were significantly associated with PFS (Log-rank test; p-value = 0.007 and p-value = 0.01, respectively), but their prognostic value was not significant (p-value > 0.05) ( Table 4). All the other clinical features were not significantly associated with survival and prognosis in the low-MATH group (Additional file 3: Table S3). None of the clinical features was significantly associated with survival and prognosis in the high-MATH group (Additional file 4: Table S4).

Exploring the PTC intra-tumoral heterogeneity
To better understand the intratumor heterogeneity of PTC, we explored the mutational landscape in different tumor stages. We calculated the mutation rates of the top 20 mutated genes for each stage and identified numerous somatic mutations. We identified several somatic mutations that frequently occur in PTC, along with the most common and well-described [18][19][20][21] with particularly high frequencies ( Fig. 3A-D; Additional file 5: Table S5). Indeed, the serine/threonine kinase BRAF gene presented relatively high mutation rates in all stages (57% in stage 1, 38% in stage 2, 69% in stage 3, 75% in stage 4). The most frequent mutation was the well-known hotspot BRAF c.1799 T > A (p.V600E). In addition, in stage 1 one patient had the simultaneous occurrence of two BRAF mutations at positions c.1799 T > A (p.V600E) and c.1800G > A (p.V600V). Also, one patient had the BRAF c.1801A > G (p.K601E) mutation, while another one had the BRAF c.1467_1481delACC TAC ACC TCA GCA (p.P490_Q494del) deletion. In stage 2, one patient had the BRAF c.1801A > G (p.K601E) mutation.
We also reported altered pathways in different PTC tumor stages. The RTK/RAS/MAP kinase signaling pathway (hereafter RTK/RAS) was the most frequently affected by somatic mutations in all stages. The overall results were reported in supplementary Additional file 6: Figure S1.

Mutational processes vary dynamically during tumor evolution
Somatic alterations in cancer genome may provide insights into the mutational processes occurring during tumorigenesis. These processes leave a peculiar pattern of mutations, called as mutational signatures [32]. We analyzed these mutational signatures across tumor stages to determine the mutational processes contributing to intratumor heterogeneity and shaping PTC tumor evolution. We identified two mutational signatures shared among all tumor stages (Fig. 4): Signature 1, related to endogenous spontaneous deamination of 5-methylcytosines, and Signature 25, related to chemotherapy treatment. While the Signature 25 was prominent in stage 1 and decreased in stage 4, the Signature 1 had an opposite trend. Furthermore, the Signature 5 was enriched in stage 2, although the etiologies of this signature have not been elucidated. Although its etiology is unknown, Signature 5 shows a clock-like behavior in many cancer types in that the number of mutations increases with the age [42]. Moreover, Signature 5 exhibits transcriptional strand bias for T > C substitutions, potentially indicating that some of these mutations arise from adducts subject to transcription-coupled repair [43]. In addition, we found that distinct defective DNA repair mechanisms might play a role in the tumor progression towards more aggressive stages. Indeed, we found a contribution of DNA mismatch repair (MMR) deficiency signatures (6,15) in all tumor stages. Signatures 6 and 15 are two of the seven mutational signatures associated with defective MMR and microsatellite instability (MSI). As seen in more aggressive forms of thyroid cancer [44], tumors harboring MMR deficiency signatures completely lack loss-of-function mutations in MMR genes (MLH1, MSH2 and MSH6). In addition, mutations in stages 2 and 3 were associated with Signature 30, related to defective DNA base excision repair due to inactivating mutations in NTHL1, while mutations in stages 1 and 4 were associated with Signature 7, related to Ultraviolet light exposure. Finally, advanced tumor stages (stage 3 and stage 4) were characterized by mutations assigned to Signature 2, associated with the activity of APOBEC family of cytidine deaminase.

Common mutated genes in 4 tumor stages
To characterize changes in the genetic architecture of PTC, we extracted the 12 common mutated genes across clinical stages, regardless of mutation type. We identified three patterns depending on the behaviors of these genes' mutation changes: (i) recurrently mutated genes in all the stages; (ii) some genes mostly mutated in early stages while disappeared in later stages; (iii) other mutated genes emerged in a dominant way in advanced stages.
In the first scenario, we found the BRAF gene, with the lowest mutation frequency in stage 2 (38%) and the highest one in stage 4 (75%). However, compared to other mutated genes, it has the highest mutation rate in all clinical stages (Fig. 5).
In the second scenario, we found genes mostly mutated in the early stages (stage 1 and stage 2). Specifically, in stage 1, the dominant mutated gene was NRAS (9%); its frequency decreased in stage 3 (5%) and stage 4 (6%). JMJD1C and PDZD2 reached a plateau from stage 2 with a mutation frequency of 2% and remained constant in advanced stages (Fig. 5). In stage 2, the dominant mutated genes were NRAS (12%) and TG (12%) genes. Both of them decreased in advanced stages: NRAS had a mutation frequency of 5% in stage 3 and 6% in stage 4, while TG had a mutation frequency of 3% in stage 3 and 4% in stage 4 (Fig. 5).
In the last scenario, we found genes mostly mutated in advanced stages (stage 3 and stage 4). In stage 3, the dominant mutated genes were NRAS (5%) and TTN (5%). Both of them had a tendency to slightly increase in stage 4 (6% for NRAS and 6% for TTN) (Fig. 5). In stage 4, the dominant mutated genes were MUC16 (8%), followed by HRAS and TTN genes showing a similar mutation frequency (about 6%). Taken together, our results suggest that: (i) BRAF driver gene shows a dominant or "wave" role across all stages, (ii) one to three driver genes exhibit a dominant role in certain stages; (iii) common genes, differently mutated across clinical stages, drive the cancer cells developing from early to advanced stages and highlight the dynamic changes of PTC oncogenesis (Fig. 5).
It is noted that mutation frequency changes were also sometimes associated with a different codon change (Table 5). Indeed, the codon changes also showed certain dynamics across stages. Specifically, some codon changes highly occurred in all stages in specific hotspots positions, such as the c.1799 T > A (p.V600E) mutation in BRAF gene, the c.182A > G (p.Q61R) mutation in HRAS gene, and the c.182A > G (p.Q61R) and c.181C > A (p.Q61K) mutations in NRAS gene. For the other genes, different codon changes dominated in different stages.

Expression level analysis using web resource
We aimed to evaluate the common mutated genes (JMJD1C, MALAT1, MUC16, PDZD2, PKHD1, RYR1, SLA and TTN) between thyroid tumor and normal samples. To this aim, we performed an in-silico analysis with UALCAN interactive computational tool [40]. The UALCAN data portal generated boxplots for each gene across 505 primary tumors and 59 normal thyroid tissues. The analysis showed that the gene expression level of JMJD1C, MUC16 and SLA were statistically downregulated while the gene expression level of PDZD2 and   RYR1 were up-regulated in thyroid cancer with respect to normal thyroid tissues (Fig. 6).
We also compared the expression of these genes across different cancer stages. Based on AJCC (American Joint Committee on Cancer) 284 patients were in stage 1, 52 in stage 2, 112 in stage 3 and 55 in stage 4. Box-whisker plots in Fig. 7 showed that JMJD1C, PDZD2, PKHD1, RYR1 and SLA presented a different expression among thyroid cancer stages, suggesting that the expression of these genes could reflect the intra-tumoral genetic heterogeneity across specific stages. Interestingly, since RYR1 was differentially expressed across all the stages, it is possible that this gene could play a key role in the pathogenesis of thyroid cancer.

Expression level of the common mutated genes in human thyroid cancer cell lines
Finally, we investigated the expression levels of the common mutated genes (JMJD1C, MUC16, PDZD2, RYR1, and SLA) in papillary thyroid cancer cell lines (K1, TPC-1 and BCPAP) compared to an immortalized normal thyroid epithelial cell line (Nthy-ori 3-1) by qRT-PCR. We found a different expression levels of these genes among the different cell lines (Fig. 8). In line with our previous findings based on TCGA, RYR1 was highly expressed in an aggressive cell line, K1, derived from metastasis of a well-differentiated PTC compared to Nthy-ori 3-1 (Fig. 8).
Finally, to verify whether the different expression level of the investigated genes in the thyroid cancer cell lines was due to their different mutational status, we purified the genomic DNA from cultured cells and sequenced the exonic region of the genes carrying high impact mutations (Table 5). We found that the high impact mutations in JMJD1C, SLA, PDZD2, identified from in-silico analysis in PTC samples, are missing in the thyroid cancer cell lines analyzed, thus suggesting  that the altered expression of these genes is due to other regulative biological mechanisms.

Discussion
Tumor heterogeneity results from the continuous accumulation of mutations during disease progression, leading to the occurrence of genetically different tumor subpopulations and influencing the clinical outcome of patients [1,4,5,45]. Due to tumor heterogeneity, the mutational landscape may differ considerably among different clinical stages of the same tumor. Based on literature results, our study is the first to directly address the evolution of intra-tumoral genetic heterogeneity during PTC progression across tumor stages. Using the whole-exome sequencing data from the TCGA-THCA cohort, we tracked the intra-tumoral heterogeneity and assessed its clinical relevance through different stages of the PTC progression to better understand the genetics of PTC carcinogenesis. We also revealed the timing of mutational processes and the dynamics of the temporal acquisition of somatic events during the lifetime of the PTC. Finally, we assayed the relationship between the mutational status and RNA expression level in a panel of thyroid cancer cell lines. Overall, the findings support the occurrence of the genetic instability in tumor progression so that the mutational landscape of each tumor stage is not stable; rather, it is subject to periodic fluctuations and multiple mutational processes shape the mutational spectra of the lesion. This appears to occur in a temporal-specific manner, thus painting a stage-specific mutational picture in a continuous dynamic manner and determining ITH to some extent. The expression levels of selected genes also showed a different expression level among cell lines characterized by different origins (primary PTC or metastasis from PTC) and by different mutation status [46]. Pan-cancer studies have shown that ITH impacts the clinical outcome and contributes to the risk of poor survival in tumor patients [14,15]. However, the association between ITH, measured by MATH score, and different tumor stages of PTC have not been thoroughly explored. Our study represents the first attempt to describe the prognostic value of MATH score. We found that ITH significantly influences the PTC patient's survival rate and as genetic heterogeneity increases, the prognosis gets worse in advanced tumor stages, mainly in stage 4. A high MATH score indicates a high percentage of heterogeneous subclones [47], which also likely makes the tumor more aggressive. Indeed, high MATH scores were associated with advanced-stage tumors and shorter overall survival in head and neck squamous cell carcinoma [10,11], with tumor stage and triple-negative or basal-like subtypes in breast cancer [48,49] and with higher risks of metastasis in stage 2 and 3 colon cancer [50]. Moreover, a high MATH score was a potential unfavorable prognostic factor in PTC, not affected by any other clinic-pathologic features. Conversely, the MATH value was not the only factor affecting the prognosis in patients with low tumor heterogeneity. Indeed, the low MATH score showed a potential M stage-dependent prognostic power, thus associated with a worse prognosis when the tumor in an advanced stage spreads to distant organs and tissues. Although the mechanisms underlying the links between high genetic heterogeneity and short overall survival cannot be deduced from these results, the strong association between ITH and disease stage supports a possible role for ITH as a prognostic biomarker in PTC. Further investigation will be required to use ITH as a novel potential biomarker for survival prediction and therapy selection.
Both somatic mutations and mutational processes can account for this heterogeneity, providing mutational fuel upon which selection can act. Some of these mutational processes are active throughout the lifetime of the cancer cell while other ones are active in a temporal-specific manner [42,51]. In accordance with this evidence, we found that the PTC genome harbored a mixture of signatures from different mutational processes. Multiple endogenous or exogenous mutational forces can operate simultaneously or successively in tumor stage diversification during PTC progression. Worthy of note is that several mutational signatures related to DNA repair deficiency (APOBEC-mediated mutagenesis, mismatch repair and base excision repair) increased later in tumor evolution or were specific for advanced stages [52], thus supporting the notion that biological processes driving the development of PTC may differ from those resulting in their progression.
The well-known driver genes (BRAF, RAS) belonging to MAPK or PI3K pathways [18][19][20][21] harbored the naturally occurring mutations and were commonly mutated across tumor stages. Specifically, mutations in BRAF and RAS were established at the early stages and confirmed in advanced stages; thus they can be considered early clonal events. On the other hand, some genes with moderate frequency emerged with a stage-by-stage expansion. Therefore, a rapidly growing clone may arise from several early key mutations, with subsequent mutations improving the fitness of tumor cell subpopulations in each stage. These subclonal alterations may likely predispose to the advancement towards more aggressive stages and potentially poorer prognosis tumors.
To better characterize the changing genetic architecture of PTC and highlight the dynamic changes of PTC oncogenesis, we followed the mutational patterns of 12 common mutated genes across clinical stages, regardless of mutation type. We identified three evolutionary paths of gene mutations that could drive PTC progression across pathological stages: (i) recurrently mutated genes in all the stages, with BRAF driver gene showing a dominant role across all stages; (ii) some genes (NRAS and TG) mostly mutated in the early stages while disappearing in the advanced stages; (iii) other mutated genes (TTN, MUC16 and HRAS) emerged dominantly in advanced stages, thus suggesting their impact in the aggressive phenotype of this tumor. It also appeared that distinct somatic events affecting the same gene occurred in distinct tumor subclones. Indeed, the mutation frequency changes, sometimes associated with a different codon change, also showed certain dynamics across stages. Specifically, some codon changes highly occurred in all stages in specific hotspots positions (BRAF c.1799 T > A, HRAS c.182A > G, NRAS c.182A > G and c.181C > A) while for the other genes, different codon changes dominated in various stages. The theory of clonal evolution of tumor [5] argues that the accumulation of mutations drives early slow-growing subclones into fast-growing subclones, thus accelerating the tumor progression. In accordance with this theory, mutations occurring at hotspot positions may drive the PTC progression towards more aggressive stages. Meanwhile, each stage acquires a specific mutational landscape that does not directly drive cancer progression but may have a strong cumulative effect, thus contributing to the aggressiveness of cancer.
Except for BRAF, HRAS, NRAS and TG, whose involvement in PTC is deeply revealed [18][19][20][21], for the other commonly mutated genes in all tumor stages, we further assayed the relationship between TCGA data, the RNA expression level of thyroid cancer-derived cell lines and their mutational status. We found different expression levels of the common mutated genes among cell lines. These differences were not due to high impact mutations harboring by these genes but probably to additional altered mechanisms such as transcriptional regulatory networks (e.g., expression of miRNAs and transcriptional factors) and epigenetic networks (e.g., DNA methylation or chromatin remodeling complex). Additionally, the BCPAP cell line is derived from a poorly differentiated PTC; the TPC-1 cell line is derived from differentiated PTC, and K1 cell line derives from metastasis of a well-differentiated PTC. Also, the mutation status of these cell lines is different: BCPAP harbored BRAF and TP53 mutations, TPC-1 cells harbored the RET/PTC1 gene rearrangement, and K1 harbored mutations in PIK3CA and BRAF genes [46]. All these data confirmed the intra-tumoral genetic heterogeneity across different papillary thyroid carcinoma models. In accordance with our results, Shen et al. [53] found an up-regulation of the RYR1 gene in 100 thyroid cancer samples compared to 64 normal thyroid tissue samples. RYR1 protein belongs to ryanodine receptors family and acts as a calcium (Ca 2+ ) release channel. Since the altered homeostasis of intracellular Ca 2+ is correlated to several hallmarks of cancer Box-Whisker plots show the expression level of the indicated genes in 59 normal thyroid tissues (in blue), 284 PTC in stage 1 (in orange), 52 PTC in stage 2 (in brown), 112 PTC in stage 3 (in green) and 55 PTC in stage 4 (in red). Statistical differences were calculated using Student's t-test. *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001 cells, the study of RYR1 could be interesting to better understand the pathogenesis of thyroid cancer. Although it is not clear how these mutations regulate RNA expression levels, the here screened genes were meaningful and worthy of further studies to better investigate the functional mechanism of how mutations work.
One of the most notable strengths of this study is the stratification of patients by tumor stage. This allowed us to obtain homogeneous patient subpopulations to be profiled for mutational landscape. However, our study has several limitations. First, we used samples taken from one small area of the tumor and at one point in and SLA was analyzed in immortalized normal thyroid epithelial (Nthy-ori 3-1) and papillary thyroid cancer cell lines (K1, TPC-1 and BCPAP) by qRT-PCR. Fold-change was determined by 2 −ΔΔCt formula and β-actin was used as internal control. Results are the means of three independent experiments ± standard error of the mean (SEM). **p ≤ 0.01 the disease course. The analysis of a single biopsy could not be representative of the whole tumor. Thus, we likely underestimated the true extent of heterogeneity within each tumor stage. Comprehensive longitudinal studies, coupled with deep multi-region sequencing, are required to better understand the mutational evolution of PTC over time. Second, our study is a retrospective analysis of a single publicly available data. Thus, a multi-data analysis and a prospective study will be needed to validate the prognostic value of the MATH score and the mutational processes driving the PTC progression towards more advanced and aggressive tumor stages. Thirdly, further in-depth studies in vitro and/or in vivo are required to gain a deeper understanding of the functional impact of JMJD1C, MALAT1, MUC16, PDZD2, PKHD1, RYR1, SLA and TTN on changes of RNA expression levels. As such, our results should be considered as hypothesis-generating findings.
Despite these limitations, our results improve the knowledge of intra-tumoral heterogeneity in different tumor stages of PTC, revealing for the first time that ITH measured by MATH could be a potential unfavorable prognostic factor in PTC and a potential risk factor for shorter survival. Moreover, this study shed light that different biological processes contributed to tumor heterogeneity of PTC, by both adding to the mutational burden and promoting molecular diversification of PTC in different tumor stages. In conclusion, our study may contribute to the development of precision medicine and the improvement of diagnostic strategies in PTC patients.

Conclusions
Emerging data, mainly due to new technologies, showed that the intra-tumoral heterogeneity represents an important feature of human cancer. Our research unveiled that the intra-tumoral heterogeneity characterized the PTC aggressiveness among different stages opening new scenarios about its impact on diagnosis, prognosis and treatment of PTC patients. Indeed, the high ITH profile associated with advanced clinical stages could be responsible for the drug resistance, thus posing clinical challenges.
The identification of specific genetic alterations in PTC patients could represent a novel tool capable to increase the diagnostic sensitivity, providing novel implications for therapeutic decisions. Thus, the implication of this work could be important in clinical practice and highlight that pathologists should consider the ITH for a precise diagnosis. The findings of this work could also give insight into the molecular mechanisms associated with the progression from the less to the more aggressive form of PTCs.