Introduction

Papillary thyroid carcinoma (PTC) derives from follicular cells and is the most common thyroid malignancy, comprising more than 80% of all malignant thyroid tumors [1,2,3]. PTC is generally considered a tumor with indolent behavior, slow progression, and low mortality rates [4]. However, a subset of cases may present an aggressive clinical course. Clinical and morphological features and gene alterations allow us to predict the aggressive behavior of this tumor, but only up to a certain extent [5]. Compared to classical PTCs, the hobnail, tall cell, columnar, and solid subtypes predominate in aggressive cases [6, 7]. Since several years, it has been ascertained that prognosis is worsened by higher disease stages, tumor size above 3–4 cm, presence of distant metastases, extrathyroidal extension, neoplastic infiltration of the resection margins, and/or older age at diagnosis [8, 9].

In terms of the impact of the molecular signature on prognosis, various studies have attempted to better characterize the genetic make-up of these neoplasms. The BRAF p.V600E mutation is the most common genomic alteration in PTC (found in approximately 70% of cases). It is currently considered that the BRAF mutation by itself is not a high-risk factor, but it has a strong adverse prognostic impact in association with TERT promoter mutations [10, 11]. Mutations of NRAS, KRAS, and HRAS are also common [12], especially in the follicular variant [13], but lack a significant and independent prognostic role. Other molecular alterations accumulate in PTC during cancer progression, such as those in TP53, PIK3CA, and AKT1 [14, 15]. In this context, TERT promoter mutations (C228T and C250T) are present in more advanced PTCs and are strongly associated with a higher risk of distant metastases, recurrence, and cancer-related mortality [10, 16].

Apart from genomic alterations, some studies investigated gene or microRNA expression profiling in PTC. However, only a few studies were designed with the specific purpose of characterizing PTCs with aggressive clinical features. The most recent data are obtained by analysis of publicly available datasets (i.e., TCGA) but are affected by the lack of a detailed characterization of cases and by an unbalanced selection of cases with or without aggressive features [17,18,19,20,21,22].

We therefore aimed to test the profiles of gene expression through the NanoString nCounter® technology in a series of well-characterized clinically aggressive PTC, matched with non-aggressive PTC cases, to identify molecular biomarkers potentially predictive of an aggressive clinical course that might be integrated with the current pathological characterization. A summary of the study design and main results is illustrated in Fig. 1.

Fig. 1
figure 1

Summary of design and main results of the study

Full size image

Material and Methods

Case Selection

A group of 43 cases with aggressive clinical behavior, defined based on the presence of metastasis at diagnosis, development of distant metastases during follow-up, and/or biochemical recurrence (> 2 ng/ml levels of thyroglobulin after TSH stimulation), was selected from a large series of PTC undergoing thyroidectomy from 2002 to 2017 at “Città della Salute e della Scienza” Hospital (Turin, Italy), San Luigi Gonzaga Hospital (Orbassano, Turin, Italy), and Mauriziano Hospital (Turin, Italy). Six cases had a biochemical recurrence, only. Thirty-seven cases had metachronous metastases. These latter cases included one case with synchronous lung metastases, only, four cases with metachronous multiorgan dissemination (including lymph nodes lung, brain, bone), and 32 cases with metachronous lymph node involvement, with or without infiltration of soft tissues in the cervical neck. These cases were matched for age, sex, and T and N parameters with 43 controls, undergoing surgery in the same period, which resulted in free of disease during the clinical follow-up (median follow-up 12.1 years, range 4.4–19 years). Characteristics of lymph node metastases in both aggressive and non-aggressive cases are reported in Supplementary Table 1.

Representative hematoxylin and eosin (H&E) stained slides were re-evaluated in all cases enrolled in the study by three of us (FM, MV, MP). The following parameters were re-assessed and recorded: histological subtype, TNM stage according to AJCC 8th edition, tumor diameter, uni- or multifocal presentation (including bilaterality), presence of tumor capsule, presence and extent of microscopic extrathyroidal invasion, presence of angioinvasion, status of surgical margins, presence of necrosis, presence of tumor-infiltrating lymphocytes (TILs), and presence of fibro-sclerotic changes. TILs were coded as absent, present “non-brisk,” and present “brisk” according to standard evaluation methods applied in breast cancer [23]. For statistical purposes, cases were grouped as TILs absent or present, in this latter case including both non-brisk and brisk patterns. Fibro-sclerotic changes were coded as mild (in the presence of thin and rare bands of collagen deposition) or extensive (in the presence of more abundant and thick fibrous bands). Representative images of the different patterns of distribution of TILs and fibro-sclerotic changes are illustrated in Fig. 2. Mitotic count was determined as the number of mitotic figures in 2 mm2. Five cases (four in the aggressive PTC group and one in the non-aggressive PTC group) met the criteria for high-grade differentiated thyroid carcinomas because of the presence of coagulative necrosis. One case had also a mitotic index of 5 in 2 mm2. However, due to their limited number, these cases were not further considered as a separate group but necrosis and mitotic index were considered individually for the sake of statistical analyses.

Fig. 2
figure 2

Representative images of the distribution of tumor-infiltrating lymphocytes and fibro-sclerotic changes in papillary thyroid carcinoma cases. Top left: tumor-infiltrating lymphocytes and fibro-sclerotic changes absent; top right: tumor-infiltrating lymphocytes present, brisk, and fibro-sclerotic changes absent; bottom left: tumor-infiltrating lymphocytes present, brisk, and fibro-sclerotic changes present mild; bottom right: tumor-infiltrating lymphocytes present, non-brisk, and fibro-sclerotic changes present extensively

Full size image

A pool of thyroid tissues devoid of macroscopically recognizable alterations from patients operated on due to follicular nodular disease was used as the baseline for gene expression testing and consisted of 3 female (aged 33, 50, and 71) and 3 male patients (aged 40, 50, and 69).

Before the study started, all cases were de-identified and coded by a pathology staff member not involved in the study, and all data were accessed anonymously. The study was approved by the local Ethical Committee (#610, date December 20th, 2017) and conducted in accordance with the principles set out in the Declaration of Helsinki. Considering the retrospective nature of this research protocol and that it had no impact on patients’ care, no specific written informed consent was required.

Molecular Analyses

For molecular testing, 24 pairs of aggressive and their non-aggressive match were analyzed by means of gene expression profiling (for a total of 48 cases tested) along with the pool of 6 samples of non-neoplastic thyroid tissue. In this sub-series, the aggressive cases were also matched with non-aggressive cases for the histological subtype.

For RNA extraction, two 10-μm-thick formalin-fixed and paraffin-embedded (FFPE) tissue sections were obtained from each tissue block and collected in sterile Eppendorf tubes by means of manual microdissection. The same procedure was adopted for normal thyroid tissues, making sure that only healthy/unaltered parenchymal areas were dissected. RNA isolation was performed using the FFPE RNA Isolation Kit (Roche Diagnostics GmbH, Mannheim, Germany), according to the manufacturer’s protocols. Total RNA concentration was assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Inc., Wilmington, DE, USA).

NanoString nCounter® technology was used to detail the transcriptomic cancer-associated landscape in PTC. Three hundred nanograms of total RNA from each sample were hybridized to the nCounter PanCancer Pathways Panel, according to the manufacturer’s instructions (NanoString Technologies, Seattle, WA, USA). The nCounter PanCancer Pathways panel includes 730 genes from 13 canonical pathways and 40 selected reference genes. The analyses were set up according to the protocol provided by the manufacturer. Expression data were normalized and analyzed with the nSolver Analysis Software (version 4.0.62). For background correction, the mean count of negative controls plus two times the standard deviation was subtracted from the counts of each gene. The means of the supplied positive controls and the geometric mean of the housekeeping genes were used to normalize the measured expression values. Both positive and negative controls were included in the panel, according to the manufacturer’s instructions. Additionally, the advanced analysis module (version 2.0.134) was used to perform differential expression analyses. Briefly, by the differential expression analysis tool (nSolver Advanced Analysis module), for each gene, a single linear regression was fit using all selected covariates to predict expression. Volcano plots were generated to display each gene’s -log10 (p-value) and log2 fold change with the selected covariate. Highly statistically significant genes fell at the top of the plot above the horizontal lines, and highly differentially expressed genes fell to either side. Horizontal lines indicated various p-value thresholds.

Statistical Analyses

All analyses were performed using Stata/MP 15.0 Statistical Software (STATA, College Station, TX). The differences in the distribution of clinical-pathological variables were analyzed using parametric and non-parametric tests (student’s t test, Pearson’s chi-square test, Bonferroni’s correction, Wilcoxon’s rank test). Overall survival (OS) was calculated from the surgical excision date of the primary tumor to the date of death or last check-up. The disease-free survival (DFS) was calculated from the date of surgical excision of the primary tumor to the date of the first relapse or last check-up. Survival curves were plotted using the Kaplan-Meier method, and the statistical comparisons were performed using the log-rank test. The Cox regression analyses were carried out on DFS and OS to calculate HRs and 95% CIs for the different study groups. Cox proportional hazards regression model was used in multivariate survival analysis to assess the independent role of parameters significant at univariate analysis. All statistical tests were two sided. A p-value < 0.05 was considered significant.

Results

Clinical and Pathological Characteristics

As shown in Table 1, aggressive PTC was more frequently multifocal (26/43, 60.5%) as compared to non-aggressive cases (16/43, 37.2%) (p = 0.031). The tumor capsule was absent in more than half of aggressive cases (29/43, 67.4%) and in only 17/43 (39.5%) of non-aggressive ones (p = 0.015). Furthermore, aggressive PTCs had more frequent TILs (26/43, 60.5%) as compared to the non-aggressive group (11/43, 25.6%) (p = 0.001). Only 10/43 (23.3%) non-aggressive PTCs had extensive fibro-sclerotic changes, as compared to 21/43 (48.8%) aggressive cases (p = 0.047). Although not statistically significant, angioinvasion was observed more commonly in aggressive (30/43, 69.8%) as compared to non-aggressive cases (21/43, 48.8%) (p = 0.077). Of note, at follow-up, 11/43 (25.6%) patients with aggressive PTC died (8 of their PTC and 3 due to other causes), as compared to 4/43 (9.3%) patients with non-aggressive PTC who died of other causes, all not related to their PTC (p = 0.047).

Table 1 Clinico-pathological features of the papillary carcinoma case series investigated
Full size table

By means of Cox univariate analyses (Table 2), in terms of DFS, the presence of angioinvasion (p = 0.019), TILs (p = 0.001), necrosis (p = 0.050), and higher mitotic counts (p = 0.001) were all correlated with an adverse prognosis. Extensive fibro-sclerotic changes were also associated with shorter DFS (p = 0.038 and p = 0.021 in the log rank test, Fig. 3). On the other hand, the presence of a complete tumor capsule was a favorable morphologic feature (p = 0.029). In terms of OS, tumor diameter larger than 4 cm (p = 0.004), presence of necrosis (p < 0.001), and higher mitotic index (p = 0.001) were parameters correlated with adverse prognosis in terms of OS. However, none of the parameters significant in univariate analyses showed independent statistical significance in multivariate analyses.

Table 2 Univariate analysis of disease-free and overall survival in the complete series of non-aggressive and aggressive PTC
Full size table
Fig. 3
figure 3

Kaplan-Meier estimates of disease-free survival (DFS) according to the fibro-sclerotic intratumoral changes (coded as absent, mild or extensive)

Full size image

OS in patients with aggressive PTC was significantly shorter than in patients with non-aggressive PTC (p = 0.038) (Fig. 4).

Fig. 4
figure 4

Kaplan-Maier estimates of overall survival (OS) of 43 aggressive versus 43 non-aggressive papillary thyroid carcinoma cases

Full size image

Gene Expression Profile of Aggressive and Non-Aggressive PTC Cases

The clinical and pathological characteristics of the 48 patients analyzed for gene expression profiling are summarized in Supplementary Table 2.

In aggressive PTC, 38 genes were significantly upregulated and 3 were significantly downregulated as compared to the non-neoplastic tissue (Supplementary Table 3; volcano plot image in Fig. 5a). In non-aggressive PTC, 80 genes were significantly upregulated and 9 were significantly downregulated as compared to non-neoplastic thyroid tissue (Supplementary Table 4; volcano plot image in Fig. 5b). Three genes, only, were significantly de-regulated in aggressive but not in non-aggressive PTC, namely two upregulated (PKMYT1 and WNT10A) and one downregulated (GLI3) (Table 3). By contrast, 45 genes were significantly upregulated and seven were downregulated in non-aggressive, but not in aggressive PTC. Thirty-six upregulated and two downregulated genes were in common between aggressive and non-aggressive PTC cases. However, even in this subset, four genes (DUSP4, INHBB, JAG2, and MLF1) were significantly upregulated in the group of non-aggressive cases as compared to the aggressive ones (Fig. 6).

Fig. 5
figure 5

Volcano plot representing differential gene expression profile in a aggressive papillary thyroid carcinomas versus normal thyroid tissue and b non-aggressive papillary thyroid carcinoma versus normal thyroid tissue

Full size image
Table 3 Distribution of deregulated genes, as compared to normal thyroid tissue, in aggressive and non-aggressive PTC cases
Full size table
Fig. 6
figure 6

Differential expression of DUSP4, INHBB, JAG2, MLF1, in aggressive versus non-aggressive papillary thyroid carcinomas

Full size image

The differential gene expression between aggressive and non-aggressive cases according to different pathways is presented in Supplementary Fig. 1. None of the genes related to DNA damage repair (POLD4, H2AFX, POLR2J, ATM, MAD2L2, PRKDC) that were upregulated in non-aggressive PTC cases were found differentially expressed in aggressive PTC (Supplementary Fig. 1a1/2).

Two genes related to the NOTCH pathway, DTX2 and JAG2, were upregulated (Supplementary Fig. 1b1/2) in both groups. All genes de-regulated in the aggressive PTC group and belonging to the transcriptional misregulation pathway were also de-regulated in the non-aggressive PTC group (Supplementary Fig. 1c1/2).

Deregulation of the hedgehog pathway involved different genes in aggressive and non-aggressive PTC. In fact, WNT10A and GLI3 were found up- and downregulated in aggressive PTCs, whereas GSK3B was found upregulated in non-aggressive PTCs (Supplementary Fig. 1d1/2).

The analyses of the WNT pathway revealed overlap in only one gene between the two groups (PPP3CB), whereas several other genes were found to be de-regulated in non-aggressive cases (CCND3, GSK3B, SFRP4, SMAD3, AXIN1, RAC2, JUN) (Supplementary Fig. 1e1/2).

Regarding the JAK-STAT pathway, upregulation of EPOR, BCL2L1, and CBL genes was found upregulated in both groups, while PIK3R2, CCND3, SOS2, OSM, JAK3, PIM1, and IL2RB genes were up- and TPO down-regulated in the non-aggressive group, only (Supplementary Fig. 1f1/2).

Considering the RAS pathway, PDGFC, KRAS, RAC2, SOS2, and PIK3R2 genes were upregulated in non-aggressive cases, only, whereas RASA4, MET, PDGFA, RIN1, TIAM1, LAT, and BCL2L1 were upregulated in both groups (Supplementary Fig. 1g1/2).

Regarding the genes labeled as “driver genes,” there was an overlap for ATRX, TRAF7, STK11, MET, SMAD2, CBL, and RUNX1, with deregulation of several other genes in non-aggressive PTC, only (Supplementary Fig. 1h1/2). A similar picture was observed for the sub-analysis of the MAPK pathway, with upregulation of TGFB1, PPP3CB, CACNB3, DUSP4, HSB1, and PDGFA in both groups, but upregulation of KRAS, MAP3K1, SOS2, DUSP5, RAC2, and CDC25B and downregulation of CACNA2D1 and JUN in the non-aggressive PTC group, only (Supplementary Fig. 1k1/2).

Finally, the PI3K pathway was significantly de-regulated in both groups, with 11 genes upregulated in both groups and several others de-regulated in non-aggressive cases (Supplementary Fig. 1l1/2).

Discussion

In the present study, we reported that clinically aggressive PTC is associated with peculiar pathological features and gene expression profiles.

PTCs in general show indolent behavior with localized disease and typically do not recur or metastasize beyond locoregional lymph nodes. A recent meta-analysis, using univariate comparisons, has identified male sex, advanced age, tumor size, multifocality, vascular invasion, extra-thyroidal extension, and lymph node metastasis as risk factors for distant metastasis in well-differentiated thyroid cancers [24]. In our series, we could not confirm sex, age, and TNM stage as parameters of aggressiveness in the aggressive vs non-aggressive PTC group comparison, since these parameters were balanced in the case match procedure. Moreover, the relatively small sample size for each group—designed to be balanced for molecular analyses—and the not complete homogeneity of compositions of the two groups (i.e., RAI treatment) limit the statistical impact of clinical pathological correlates, as demonstrated by the low level of statistical significance in DFS and OS of non-aggressive and aggressive PTC cases. We could confirm multifocality and angioinvasion as parameters significantly more prevalent in aggressive PTC, the latter also showing a significant association with shorter disease-free survival.

Concerning additional pathological parameters, we found that the presence of TILs was more frequent in aggressive PTCs as compared to non-aggressive cases. This finding is in contrast with previous data on the association of immune cell infiltrate (including B and T cells) and a more favorable outcome [25,26,27]. As a limitation of our study, we did not perform subtyping of the tumor-associated immune compartment, but it is conceivable that our case series selection approach influenced our results as compared to literature studies. Moreover, our data are strengthened by the significant association of the presence of TILs with shorter disease-free survival (with a trend to significance also in terms of shorter overall survival).

Another feature associated with an aggressive clinical course in PTC was the presence of fibro-sclerotic tissue within the tumor. Our findings are supported by a few previous studies [28, 29] but are in contrast with a recent report on a large unselected series [30]. However, the value of our results is supported by a disease-free survival analysis that showed a significant adverse impact of the presence of extensive fibrotic changes. These fibrous tissue modifications were appearing as trabecular, net-like structures permeating tumor architecture. In a setting of parathyroid tumors, the presence of thick fibrous bands is commonly referred to as a sign of malignancy [31, 32].

Finally, increased mitotic activity and the presence of necrosis were identified as strong adverse prognostic indicators, both in disease-free and overall survival analyses, although their distribution was not statistically different comparing aggressive and non-aggressive PTC cases. This latter finding is strongly supportive of the concept of the high-grade differentiated thyroid cancer group, as a specific prognostic category within differentiated thyroid cancers and is in line with the novel category proposed by the fifth edition of the WHO book on endocrine tumors, which embeds aggressive differentiated carcinomas and poorly differentiated carcinomas in a single group named high-grade thyroid carcinomas [33, 34].

A second part of our study was focused on the definition of gene expression profiles in a subset of aggressive and non-aggressive PTC cases, selected from the pathological series and matched also for the histological subtype. Despite other studies employed NanoString nCounter technology to explore the molecular make-up of thyroid malignances [35,36,37], to the best of our knowledge, this is the first study using the nCounter® PanCancer Pathways Panel to explore different canonical pathways in correlation with PTC aggressiveness. Our study falls in the scarce literature aimed at testing gene expression profiles as determinants of aggressiveness in PTC. The weakness of our approach is mostly based on the targeted approach (although covering over 700 genes), the lack of genomic data of the cases investigated, and the lack of RNA validation using an alternative approach or protein analysis by means of immunohistochemistry. By contrast, the strengths of our study are the robustness of the nanostring technology and the study design in terms of case selection. In fact, as compared to similar studies, we designed a comparative study on aggressive and non-aggressive PTC cases based on a strongly balanced and matched selection. A recent study based on TCGA data identified a specific signature of immune-related genes (namely HSPA1A, NOX5, and FGF23) associated with prognosis in PTC, but using data from a series of nearly 500 PTC cases with 3.2% of patients, only, died of disease and less than 2% of cases with distant metastases [38]. Another very recent study, also based on TCGA data, compared gene expression profiles in 455 non recurrent vs 46 recurrent PTC tumors specimens, identifying differential expression profiles in 40 genes with a significant influence of the tumor genotype [19]. In this study, recurrent PTC cases with BRAF-like signature significantly overexpressed—among other genes—FN1, ITGA3, and MET and downregulated TPO and TG.

Our study clearly shows that clinically aggressive and non-aggressive PTC cases in our series consistently share a gene expression signature as compared to normal thyroid tissue. In terms of pathway analysis, up-regulation of the RAS, MAPK, and PI3K pathways was consistent in both groups, although with an enrichment of de-regulated genes in the non-aggressive PTC group. However, a larger number of genes was deregulated in non-aggressive PTC as compared to normal tissues than those deregulated in aggressive PTC cases. The larger set of de-regulated genes depicted in non-aggressive PTC cases shows a more significant upregulation of genes belonging to the DNA damage repair and JAK/STAT pathways and a downregulation of genes belonging to the WNT pathway. Interestingly, the hedgehog pathway was differentially de-regulated in aggressive PTC as compared to non-aggressive PTC cases. In fact, whereas WNT10A and GLI3 were significantly up- and down-regulated in aggressive PTC, GSK3B was upregulated in non-aggressive PTC cases.

These findings suggest that these genes are likely involved in molecular mechanisms of aggressiveness. In fact, WNT10A/beta-catenin pathway activation has been demonstrated to promote cellular characteristics of aggressiveness, such as proliferation and migration in thyroid cancer [39]. Data on GLI3-driven molecular mechanisms in thyroid cancer are not available. GLI3 activation/overexpression has been shown to promote tumor progression in colon cancer [40] but not in other models such as gallbladder or cervical cancer [41, 42] thus supporting tumor type-specific activities. On the other hand, GSK3B expression has been described to inhibit the progression of several cancer models, such as, among others, lung [43], gastric [44], and hepatocellular [45] cancers. Finally, PKMYT1 was significantly upregulated in aggressive PTC cases, only. Although no data are available on the role of this gene in thyroid cancer, PKMYT1 has been shown to exert an important effect on tumor immunity and progression in several cancer models [46].

In conclusion, we demonstrate that clinically aggressive PTC is characterized by peculiar pathological features, including an increased TILs infiltration and a more prevalent occurrence of intratumoral fibrosis, supporting the value of a careful reporting of these parameters in the PTC diagnostic work-up. Moreover, although non-aggressive PTC in our series showed a larger number of deregulated genes as compared to aggressive PTC cases, the existence of a set of differentially regulated genes in the two groups suggests that different molecular mechanisms are active in promoting clinical aggressiveness in PTC with special reference to the impairment of the hedgehog pathway. Although limited by the relative small sample size and biases in case selection, these results have a prime impact highlighting the need to evaluate the stromal component as part of the diagnostic workup for PTC, with special reference to the presence of fibrosclerotic changes, a pathological characteristic that is currently not considered in standardized pathology reporting for thyroid cancer [47]. Moreover, gene expression data pave the way for molecular validation studies targeting regulators of the hedgehog pathway to identify novel biomarkers of aggressiveness.