Introduction
Papillary thyroid cancer (PTC) is the most common type of differentiated thyroid cancer deriving from the follicular cell of the thyroid gland, accounting for about 80% of all thyroid malignancies [1]. Over a decade, high-throughput methods, including next-generation sequencing and microarrays, have been used for the evaluation of PTC genotype and gene expression profile to find characteristic features related to neoplastic transformation and turn it into diagnostic and clinically valuable information [2-5]. The analysis of 500 PTC tumors, based on next-generation sequencing, under the genotyping project Cancer Genome Research Network, confirmed four major molecular events present in the majority of PTC cases: BRAFV600E, RAS, and TERT point mutations as well as RET/PTC rearrangements [4]. Although the PTC transcriptome is well characterized [3, 6-8] and correlates with some prognostic and genotypic factors [3, 4], biological mechanisms of neoplastic transformation in PTC are not unequivocally elucidated. Nevertheless, these data do not explain the whole biology of PTC tumorigenesis and are only to some extent helpful in diagnostics (BRAF, TERT, RET/PTC rearrangements) [9-14] and prognosis (BRAF, TERT) [15-18].
It is well known that the tumor microenvironment (TME) contributes to cancer development and promotes angiogenesis, invasion, metastasis, and chronic inflammation [19]. The possible role of TME in response to chemotherapy [20] or antiangiogenic therapy [21] has been addressed recently. In healthy tissues the wound healing process is promoted by the activation, migration, and proliferation of fibroblasts, whereas cancer-associated fibroblasts play a role in tumor growth and therefore are crucial for disease progression [19, 22].
Already published data which addressed the interaction between PTC and stromal tissue, lymphocytic infiltrate, or normal thyrocytes and their influence on gene expression profile are scarce. The results of our previous analysis demonstrated a strong impact of stromal cells on PTC gene expression profile. We observed that immunity-related genes provided the most intensive confounding signal, probably related to lymphocyte tumor infiltration [3]. The proportion of particular cells may vary between tumors, which may also change the pattern of a transcriptome analyzed in whole tissue slide.
In the present study, we aimed to assess the impact of TME on gene expression profile in PTC. To achieve this goal, we decided to evaluate the gene expression pattern of PTC and normal thyroid cells, isolated by laser capture microdissection and in whole tissue slides. We believe that such an approach allows for wider analyses and, in consequence, for better understanding of the biology of PTC.
Materials and Methods
Analysis of PTC Transcriptome
Tumor tissues samples were removed during surgical resection, snap-frozen on dry ice, and stored at –80°C until use.
Fifteen microdissected PTC samples, 11 samples containing normal thyrocytes taken from contralateral thyroid lobe and 15 pairs of whole slide PTC/macroscopically unchanged tissue, were subjected to transcriptome profiling with oligonucleotide microarrays HG-U133 PLUS 2.0 (Affymetrix, currently Thermo Fisher Scientific USA). Total RNA was extracted from 1 × 106 μm3 of microdissected tissue (PALM system) and extracted by RNeasy Micro Kit (Qiagen GmbH, Hilden, Germany). Ten to thirty nanograms of total RNA were used for double round amplification (Two Cycle Target Labeling Kit, Affymetrix/Thermo Fisher, USA). From the whole slides 100 ng of total RNA were used for one round amplification (GeneChip 3’ IVT Express Kit, Affymetrix/Thermo Fisher, USA). All reactions, washing, staining, and scanning were performed according to the manufacturer’s manual.
Validation of Obtained Results by Quantitative Real-Time Polymerase Chain Reaction
RNA subjected to quantitative real-time polymerase chain reaction (qRT-PCR) validation was extracted from 31 pairs of whole slide PTC/normal thyroid, 12 microdissected PTC samples, 6 samples of microdissected normal thyrocytes, 12 PTC slides, and 13 whole normal thyroid slides.
Total RNA was extracted using RNeasy Mini Kits (Qiagen GmbH) according to the manufacturer’s protocol. Reverse transcription reaction was performed using Qiagen SensiScript RT Kits and 5 ng of input RNA according to the manufacturer’s protocol. All genes were amplified by the 7900HT Fast Real-Time PCR system (Life Technologies/Thermo Fisher, USA). For amplicon design, Roche Universal Probe Library was used. For normalization of the qRT-PCR data, EIF3S10, HADHA, and UBE2D2 genes were used. The normalization factor was obtained using the GeNorm applet for Microsoft Excel.
Statistical Analysis
Microarray Data Analysis. All bioinformatic analyses of gene expression profile were performed in an R/Bioconductor environment. The GCRMA algorithm was used for normalization. Unsupervised analysis was perform using principal component analysis. Supervised analysis was based on gene expression comparison using Student t test. Differentiated transcripts were selected in three separate comparisons: (1) transcripts differentially expressed in microdissected and whole slides – nominal p < 0.001 (FDR < 2%) and fold change >4 (both directions); (2) transcripts differently expressed in microdissected samples (criteria as above) and not changed in whole slides – nominal p > 0.01 (FDR > 5%); (3) transcripts differentially expressed in whole slides and not changed in microdissected samples (criteria are given in point 2). Hypergeometric test was used for association of Kyoto Encyclopedia of Genes and Genomes (KEGG) biological pathways.
qRT-PCR Analysis. Normalized qRT-PCR data were assessed by the nonparametric Mann-Whitney U test adjusted by Bonferroni correction.
Results
In the first step, in order to assess the main source of variation, principal component analysis was performed separately in two groups: in a set containing thyrocytes obtained by microdissection and in a set containing whole slides of PTC tumor and normal thyroid (Fig. 1a, b). In both groups the main source of variation was the difference in gene expression profile between PTC samples and normal thyroid samples. Principal component analysis carried out in a combined set containing both microdissected cells and whole slides also revealed the difference in the transcriptome between PTC and normal thyroid. However, the differences between microdissected cells and whole slides were also visible. These differences were shown in the first principal component (Fig. 1c).
Fig. 1.
Principal component analysis showed that the main variability in gene expression profile was related to sample used (microdissected cells vs. whole slide) (a, b) and to sample origin (PTC vs. normal thyroid regardless of material used – microdissected cells or whole slide) (c). PC, principal component; PTC, papillary thyroid cancer.
In the second step, supervised analysis was carried out in the following comparisons to select: (1) transcripts differentiating PTC and normal thyroid in microdissected cells and whole tissue slides; (2) transcripts differentiating PTC and normal thyroid in microdissected cells and not changed in whole tissue slides; (3) transcripts differentiating PTC and normal thyroid in the whole slides and not changed in microdissected cells.
The first part of the supervised analysis identified 558 transcripts showing significant differences in expression between PTC and normal thyroid both in microdissected samples and whole slides (Fig. 2a; p < 0.001, FDR < 2%, and fold change > 4 in both directions).
Fig. 2.
Hierarchical clustering performed on the transcripts differentiating PTC and normal thyroid in microdissected cells and whole slides (a), in microdissected cells only (b), and in whole slides only (c). PTC, papillary thyroid cancer.
The goal of the second part of the supervised analysis was to select genes differentially expressed in the whole sections (nominal p < 0.01, FDR < 2%, and fold change > 4 in both directions) and not altered in microdissected cells (p > 0.01). A total of 159 significantly deregulated transcripts were identified (Fig. 2b).
The third step of supervised analysis focused on the selection of transcripts with dysregulated expression only in microdissected cells (nominal p < 0.01, FDR > 2%, and fold change > 4 in both directions) but not significant in whole tissue slides (p > 0.01). This analysis identified 106 differentially expressed transcripts (Fig. 2c).
Hierarchical clusterization of transcripts selected in a group combined of microdissected cells and whole tissue slides correctly separated PTC and normal thyroid samples regardless of the sample analyzed. However, clusterization of differentiating transcripts obtained in the analyses of whole slides or microdissected cells alone separated samples in an incomplete manner (Fig. 2).
The last part of the study aimed to validate obtained microarray data. It was performed with the use of qRT-PCR on an independent set of samples. Eleven transcripts were selected for qRT-PCR validation, among them seven transcripts differentially expressed between PTC and normal thyroid in both groups (microdissected cells and the whole slides) and four transcripts significantly differentiating PTC and normal thyroid only in microdissected cells (Table 1).
All seven selected genes demonstrating significant differences in expression between PTC and normal thyroid were confirmed by qRT-PCR (Fig. 3). Only one among four transcripts differentially expressed between PTC and normal thyroid in microdissected cells was not positively validated by qRT-PCR (FOS). Three other genes (PTCSC1, CTGF, EGFR1) were confirmed on an independent set of samples. In addition, the EGFR1 gene showed a significant difference in expression between PTC and normal thyroid in whole slides (Fig. 4).
Fig. 3.
Validation of seven genes differentially expressed between PTC and normal thyroid in both laser microdissected cells and whole slides in an independent set of samples. T, N, normal thyroid; PTC, papillary thyroid cancer; T, tumor.
Fig. 4.
Expression of four genes differentiating between PTC and normal thyroid only in microdissected cells, validated by qRT-PCR. N, normal thyroid; ns, not significant; PTC, papillary thyroid cancer; qRT-PCR, quantitative real-time polymerase chain reaction; T, tumor.
All three transcript sets were also analyzed by the KEGG concerning biological pathways to determine possible implications of transcripts associated with microdissected cells and whole tissue slides. The pathways characteristic for each comparison are presented in Table 2.
Table 2.
KEGG biological pathways significant in PTC: microdissected cells and whole slides, only in microdissected cells, and only in whole slides
Discussion
The evaluation of the transcriptome of selected transformed neoplastic thyrocytes and whole PTC slide compared to follicles of normal thyroid gland and macroscopically unchanged thyroid tissue demonstrated that the main source of variation was the difference between the whole section and the cell subjected to microdissection, while the differences in gene expression profile between neoplastic samples and normal thyroid remained significant. It means that the gene expression we observed during the analyses of the whole section was related to other cells than neoplastic ones. It has also been demonstrated that stromal cells interact with neoplastic cells. It is evidenced by the fact that genes selected in a set of microdissected cells and in a set of the whole slides allowed for a correct classification of PTC cells and normal thyroid cells, regardless of the analyzed material – microdissected cells or whole slide. Finally, eleven genes were selected to be verified in an independent set of samples; among them only four genes differentiated microdissected neoplastic cells and normal cells. Two genes, whose expression was significantly lower in PTC cells (PTCSC1 and CTGF), were confirmed in an independent validation set. Regarding two other genes, FOS (FBJ murine osteosarcoma viral oncogene homolog) was not confirmed by the validation, whereas for early growth response 1 (EGR1) a significant result was also obtained in the analysis of the whole slides.
The remaining genes (TFF3, FN1, MPPED2, MET, KCNJ2, TACSTD2, and GALE) showed varied expression in microdissected thyrocytes and in the whole tumor slides. It is worth noticing that these genes are frequently involved in classifiers characterizing PTC profile [3, 4, 8].
The papillary thyroid carcinoma susceptibility candidate 1 gene (PTCSC1) is a noncoding RNA belonging to a group of long noncoding RNA, >200 nucleotides in length [23]. Recent studies demonstrated that long noncoding RNA could function both as oncogenes or suppressor genes [23].
The PTCSC1 gene is localized in a 3 intron of the Src-like adapter (SLA) gene (sense strand), which is in turn found in an intron of thyroglobulin gene (antisense strand) [24, 25]. In 2009, He et al. [24] proposed the PTCSC1 gene as a potential candidate gene related to a predisposition to PTC. This study also demonstrated lower expression of this gene in the majority of samples obtained from members of the family with a frequent occurrence of PTC. However, the authors did not find any mutation in a gene sequence which could be linked with predisposition to PTC or its lower expression in PTC.
The connective tissue growth factor gene (CTGF/CCN2) belongs to the CNN gene family. It encodes a mitogen (mitosis-promoting factor) involved in cell division, differentiation, and cell adhesion [26]. The gene plays an essential role in the biology of connective tissue. It is expressed in most of the stromal cells (fibroblasts, endothelial cells, smooth muscle cells), stimulates the production of extracellular matrix (ECM) [27], and participates in wound healing. Its high expression is observed in stromal cells adjacent to different epithelial neoplasms, including breast or prostate cancer [26]. It is believed that the CTGF gene is a marker of wound healing and also a factor related to the progression of neoplastic disease. The study, carried out on a mouse prostate cancer xenograft model, showed that expression of the CTGF gene was related to the TME, that it played a crucial role in tumor-reactive stroma, and that it was a stroma-related factor enhancing carcinogenesis [26].
The EGFR1 gene encodes a nuclear transcription regulatory protein. Trinh et al. [28] demonstrated that a decreased expression of this gene in tissues under hypoxic conditions promoted angiogenesis. EGFR1 may be regulated both by the MAPK/ERK pathway and hypoxia-inducible factor-1α (HIF-1), which activated EGFR1 promoter. The activation of MAPK pathway is characteristic for PTC [29], whereas the activation of the HIF-1 gene under hypoxic conditions is frequently noticed in solid tumors, among them thyroid cancer [30]. Angiogenesis is essential for tumor growth and progression, which is probably reflected by MAPK activation observed in PTC. Similarly, hypoxia leads to activation of proangiogenic factors, including the EGFR1 gene.
The remaining genes (TFF3, FN1, MPPED2, MET, TACSTD2, KCNJ2, and GALE), whose expression characterized the whole PTC slide and microdissected neoplastic transformed thyrocytes, reflect the reaction of the TME. The majority of them (MET, FN1, KCNJ2, TACSTD2, and GALE) showed significantly increased expression in normal thyroid tissue. The MET gene is a protooncogene encoding receptor tyrosine kinase which plays a role in signal transduction from the ECM to the cytoplasm [31]. Activation of the MET receptor is a consequence of a binding to hepatocyte growth factor (HGF) ligand. HGF ligand is a stromal-derived factor secreting by fibroblasts of the TME. Matsumoto et al. [31] demonstrated in 1989 that the TME was necessary to activate the MET receptor. The Japanese group also proved that fibroblasts were crucial for the invasion of oral squamous cell cancer cells and that HGF ligand was a factor responsible for this process [31, 32]. Cell activation by stromal fibroblasts via stimulation of the MET receptor was also noticed in other malignancies, including breast, prostate, esophageal, and gallbladder carcinomas [33-36]. Moreover, the MET receptor showed a higher activity in hypoxic areas of the tumor. It could be related to a higher invasiveness of tumors showing a high level of hypoxia [37]. This dependence was also noticed in PTC. A high activity of MET protein was more pronounced at the tumor periphery, in regions of cell invasion promoted by HIF-1 [38].
Similar cancer-stromal interactions concern the expression of the fibronectin 1 gene (FN1). This gene encodes fibronectin, an ECM glycoprotein secreted by fibroblasts and present on the cell surface of fibroblasts or isolated from human plasma (fibronectin circulating in human plasma derives from hepatocytes). By binding to collagen, receptors on the cell surface (integrins), plasma fibrin, and glycosaminoglycans (heparin), FN1 is involved in cell adhesion and activates cell migration [39, 40]. It could be expressed by epithelial cells [41]. High expression of the FN1 gene is frequently observed in PTC. Following the results of different studies concerning PTC transcriptome, the FN1 gene has been proposed as a molecular PTC marker in numerous papers [3, 42-45]. However, attempts to use FN1 protein in immunostaining were not so successful. Expression of FN1 protein is observed in up to 40% of benign follicular lesions. Better results are obtained if the expression of several proteins is evaluated at the same time [46]. Tumor-stroma crosstalk, including FN1 protein expression, was also analyzed in melanoma and breast carcinoma [47, 48]. These studies demonstrated an increased synthesis of FN1 protein and collagen and their release to the ECM by stromal cells, such as myofibroblasts [47, 48]. Moreover, it is believed that stromal fibroblasts, stimulated by the tumor cells, secrete different factors promoting tumor progression and growth [49]. The presence of FN1 protein in epithelial cells is also linked to the so-called epithelial-mesenchymal transition process, based on changes in biochemical and biological phenotype of epithelial cells, playing an essential role in physiology and pathology, including cancer progression. The transition of epithelial cells to mesenchymal ones leads to loss of cell-cell adhesion, cell polarity, and acquisition of migratory and invasive properties [50, 51].
Another gene characterized by a significantly higher expression in the whole slide and selected PTC cells was tumor-associated calcium signal transducer 2 (TACSTD2). First data concerning its function pointed to its potential role as an oncogene [52]. Its expression was observed in different neoplasms, including PTC. However, according to recent data, the expression of the TACSTD2 gene was also noticed in normal, healthy cells [53], whereas in neoplastic cells, the TACSTD2 gene did not demonstrate oncogenic activity, but was related to the pro-cancer effect of the TME and epithelial-mesenchymal transition [54].
The potassium inwardly-rectifying channel, subfamily J member gene (KCNJ2) encodes Kir2.1 protein, a potassium channel that allows K+ to move into the cell. In physiology, Kir2.1 protein is present in the majority of healthy cells, including neurons, muscle cells, and the immunological system. It is also expressed by neoplastic cells [55]. The UDP-galactose-4-epimerase gene (GALE) encodes UDP-galactose-4-epimerase that converts UDP-galactose to UDP-glucose. Its overexpression in PTC was confirmed by our previous study [43].
Significantly decreased expression of the trefoil factor 3 gene (TFF3) was also reported in other studies. Takano and Yamada [56], who analyzed the gene expression profile in thyroid tumors, proposed TFF3 as a marker of malignancy. Expression of this gene was observed in normal thyroid tissue, benign follicular neoplasms, and minimally invasive follicular carcinoma, whereas its decreased expression was characteristic for PTC. The authors tried to evaluate the expression of the protein encoded by TFF3 in immunostaining, but precise assessment of the intensity of such a reaction was not possible.
Decreased expression of the metallophosphoesterase domain containing 2 gene (MPPED2) in the whole PTC slide was demonstrated for the first time in our previous study [43]. MPPED2 was reported as a top gene which allowed for a proper classification of benign and malignant thyroid tumors in 96.7% of samples. This gene encodes a metallophosphoesterase and demonstrated high expression in the brain during mammal fetal development. New data regarding its potential role as a suppressor gene substantially involved in the neoplastic transformation of follicular thyroid cells have been published recently [57].
To date, only one study using whole-transcriptome profiling in PTC cells selected by laser microdissection has been reported [58]. However, one should stress that the comparison of gene expression between PTC cells and thyrocytes obtained from macroscopically unchanged thyroid tissue was not based on any statistical tests considering significance level or multiple-comparison correction. The authors selected up- or downregulated genes based on an expression ratio >5.0 or <0.2 in at least 50% of informative cases. There are no data regarding the number of analyzed microdissected cells. The FN1 gene was listed as the first among all genes commonly upregulated in thyroid cancer. The authors reported overexpression of KCNJ2, whereas the FOS gene was downregulated.
One should notice that most analyses carried out with the use of microdissected cells regarding the thyroid focused on the evaluation of single molecular markers [59-63], whereas our approach allows for a more comprehensive look at the PTC transcriptome and surrounding microenvironment composed of different cells and ECM.
It is currently believed that cancer progression is promoted by the orchestrated interaction between malignant neoplastic cells and their environment. This interaction is, in some part, related to molecular alterations of mutated genes, which specifically drive the composition of the TME. It is widely accepted that dynamic crosstalk between the TME and tumor cells is involved in cancer progression, metastatic spread, recurrence, and response to drugs applied [64]. However, most studies which analyze the contribution of TME to cancer progression concern its single component and mainly focus on immune response [65, 66].
The crosstalk between the tumor and its microenvironment might also be visible in pathways activated by transcripts significantly differing in PTC microdissected cells and whole slides. The top pathway was “metabolism of xenobiotics by cytochrome P450,” which mainly related to drug metabolism in cells affected by the TME, while the top KEGG pathway characterized for genes deregulated only in microdissected cells was “vascular muscle contraction.” It refers to vascular smooth muscle cells that in pathological conditions may transfer to a proliferative state [67]. Regarding genes differentiating PTC and normal tissue observed only in whole tissue slide, the top KEGG pathway was “complement and coagulation cascades,” crucial for the inflammatory network [68]. It could be related to immune response as an initial phase of a neoplastic transformation in the interaction between the tumor and its microenvironment.
For a better understanding of gene interactions between tumor cells and microenvironment, in the present study, we compared the gene expression profiles of PTC and normal thyroid in microdissected cells and whole tissue slides. To our best knowledge, this is the first study aimed at identifying and describing transcripts interacting between PTC and other cells composing the TME. We hope that our results on one hand reflect the complexity of tumor-microenvironment crosstalk and on the other hand contribute to a better understanding of PTC biology.
To conclude, most of the identified genes differentiating PTC and normal thyroid are related to the tumor-microenvironment interaction. One should answer the question of whether analysis of selected cancer cells is sufficient to know tumor biology or whether analysis of the whole tumor may better characterize tumor molecular profile.
Statement of Ethics
All subjects participating in this study gave their written informed consent and the study protocol was approved by the institute’s committee on human research.
Disclosure Statement
The authors have no conflicts of interest to declare.
Funding Sources
The study was supported by the Polish National Center of Research and Development MILESTONE project – Molecular diagnostics and imaging in individualized therapy for breast, thyroid and prostate cancer, grant no. STRATEGMED 2/267398/4/NCBR/2015.
Author Contributions
M. Oczko-Wojciechowska: study concept, study supervision, evaluation of gene expression profiles, data analysis, writing and correction of the manuscript. A. Pfeifer: statistical analyses. M. Jarzab: study concept, data analysis. M. Swierniak: statistical analyses. D. Rusinek: microdissection. T. Tyszkiewicz: qPCR validation. M. Kowalska: evaluation of gene expression profiles. E. Chmielik and E. Zembala-Nozynska: histopathological analysis. A. Czarniecka: surgical procedures, tumor sampling. B. Jarzab: study concept, study supervision, data analysis. J. Krajewska: study concept, study supervision, writing and correction of the manuscript. All authors discussed the results and contributed to the final version of the manuscript.
