Gene Expression Profiling of Papillary Thyroid Carcinoma Identifies Transcripts Correlated with BRAF Mutational Status and Lymph Node Metastasis Free

Purpose: To identify papillary thyroid carcinoma (PTC)–associated transcripts, we compared the gene expression profiles of three Serial Analysis of Gene Expression libraries generated from thyroid tumors and a normal thyroid tissue.

Experimental Design: Selected transcripts were validated in a panel of 57 thyroid tumors using quantitative PCR (qPCR). An independent set of 71 paraffin-embedded sections was used for validation using immunohistochemical analysis. To determine if PTC-associated gene expression could predict lymph node involvement, a separate cohort of 130 primary PTC (54 metastatic and 76 nonmetastatic) was investigated. The BRAF^V600E mutational status was compared with qPCR data to identify genes that might be regulated by abnormal BRAF/MEK/extracellular signal-regulated kinase signaling.

Results: We identified and validated new PTC-associated transcripts. Three genes (CST6, CXCL14, and DHRS3) are strongly associated with PTC. Immunohistochemical analysis of CXCL14 confirmed the qPCR data and showed protein expression in PTC epithelial cells. We also observed that CST6, CXCL14, DHRS3, and SPP1 were associated with PTC lymph node metastasis, with CST6, CXCL14, and SPP1 being positively correlated with metastasis and DHRS3 being negatively correlated. Finally, we found a strong correlation between CST6 and CXCL14 expression and BRAF^V600E mutational status, suggesting that these genes may be induced subsequently to BRAF activation and therefore may be downstream in the BRAF/MEK/extracellular signal-regulated kinase signaling pathway.

Conclusion:CST6, CXCL14, DHRS3, and SPP1 may play a role in PTC pathogenesis and progression and are possible molecular targets for PTC therapy.

Thyroid cancer is the most common endocrine malignancy, with papillary thyroid carcinoma (PTC) and follicular thyroid carcinoma (FTC) being the most frequent tumor subtypes. Although PTC and FTC arise from the follicular epithelial cells, these tumors differ in biological and clinical behavior. Defining molecular differences and similarities between these subtypes of thyroid tumors could provide insights into their pathogenesis and may help identify new therapeutic approaches.

PTC represents 80% of thyroid cancers and therefore is the most extensively investigated thyroid tumor. Thus far, the main molecular finding is the activating mutations of RAS/BRAF or rearrangement of RET, both of which activate the mitogen-activated protein kinase pathway. These genomic events appear to be mutually exclusive and are the earliest mutations occurring in PTC (1–3). In vitro analysis showed that RET/PTC/RAS/BRAF induce changes in the expression of an overlapping set of genes (4), whereas other have suggested that the various alterations in this pathway produce distinct patterns of gene expression (5, 6). Despite the importance of the mitogen-activated protein kinase pathway in PTC, it is likely that there are other mutations and subsequent alterations in gene expression that define this cancer.

The premise of this study is that the identification of novel PTC-associated transcripts, combined with previous molecular data, will lead to a better understanding of the events involved in pathogenesis and progression of PTC. We used Serial Analysis of Gene Expression (SAGE) to identify PTC-associated molecular alterations. SAGE was employed because of its ability to accurately produce comprehensive expression profiles and because we have previously conducted independent in-deep SAGE libraries of PTC, FTC, normal thyroid, and follicular thyroid adenoma (FTA; refs. 7, 8).

In addition to observing the known transcripts altered in PTC, we identified new transcripts expressed in PTC. We selected five genes [cystatin E/M (CST6), chemokine (C-X-C motif) ligand 14 (CXCL14), dehydrogenase/reductase (SDR family) member 3 (DHRS3), Brevican (BCAN), and secreted phosphoprotein 1 (osteopontin; SPP1)] to validate by quantitative PCR (qPCR). Immunohistochemistry was done to determine the immunolocalization of CXCL14. As most of the transcripts were associated previously with progression of other human cancers, we assessed their expression in a panel of 130 primary metastatic and nonmetastatic PTC. Additionally, to investigate if these genes expression might be activated by BRAF/MEK/extracellular signal-regulated kinase/mitogen-activated protein kinase signaling, we combine the mutational status of BRAF gene, the most frequent mutated gene in PTC, with the qPCR expression data of 54 metastatic and 76 nonmetastatic PTC to look for associations.

Tissue samples. A total of 57 thyroid tumors and 14 normal adjacent thyroid tissues were used to find PTC-associated transcripts. Samples were obtained from patients undergoing thyroid surgery at Hospital São Paulo, Universidade Federal de São Paulo and Hospital das Clínicas, Universidade Estadual de São Paulo. Tissue specimens were frozen immediately after surgical biopsy and stored at -80°C until use. Final histologic classifications were obtained from paraffin-embedded sections and comprised 28 PTC (15 classic and 13 follicular variant), 17 FTC, and 12 FTA. For protein analysis, an additional set included paraffin-embedded sections selected from the archives of the Department of Pathology, Federal University of São Paulo. H&E-stained sections were reviewed by a pathologist and included 42 PTC, 10 FTC, 19 FTA, and 28 matched-normal thyroid tissues. To evaluate if the candidate genes were associated with tumor progression, a panel of 130 primary PTC (54 metastatic and 76 nonmetastatic) was included in a separated validation group. The study was approved by the ethics and research committees from both universities and was conducted in accordance with the Declaration of Helsinki principles.

Cell lines. Follicular (WRO) and papillary (NPA) thyroid carcinoma cell lines were used in this study (generous gifts of Dr. Alfredo Fusco, University Federico II). WRO and NPA were grown in DMEM (Invitrogen) supplemented with 10% fetal bovine serum (Invitrogen), 100 units/mL penicillin, and 100 μg/mL streptomycin in a humidified incubator containing 5% CO₂ at 37°C. NPA cell line is known to have several specific genetic abnormalities including homozygous BRAF mutation (V600E). WRO has a wild-type BRAF (1).

Comparison of papillary and follicular thyroid SAGE libraries. The papillary and follicular SAGE libraries, constructed for previous studies (7, 8), were compared for this study using SAGE 2000 software. The full set of tag counts for all thyroid libraries is available for downloading or analysis at the Cancer Genome Anatomy Project SAGE Genie Web site (7).⁶

Pair-wise and Monte Carlo simulations were used to identify transcripts that were statistically different at P ≤ 0.001 and minimum tags account of 5. Transcripts that were overexpressed in the PTC library were assessed in the FTA and normal thyroid libraries (7). Genes with large expression differences in PTC library were selected for further validation by qPCR (Table 1).

RNA extraction, cDNA synthesis, and qPCR. RNA isolation and cDNA synthesis were done as reported previously (7). An aliquot of cDNA was used in 20 μL PCR containing SYBR Green PCR Master Mix (Applied Biosystems) and 200 nmol/L of each primer for the target or reference gene. The PCR was done for 40 cycles of a four-step program: 94°C for 15 s, annealing temperature for 20 s, 72°C for 20 s, and a fluorescence-read step for 20 s. PCR conditions and primer sequences are described in Supplementary Table S1. qPCR was done in triplicate and the threshold cycle was obtained using Applied Biosystems software and was averaged (SD ≤ 1). Gene expression was normalized using the average of two-control genes (ribosomal protein S8 and QPC) shown by SAGE to be at equivalent levels in all four SAGE libraries. Fold changes were calculated as described previously (7, 9).

Immunohistochemical analysis. Immunohistochemical staining was done in a total of 71 formalin-fixed paraffin-embedded thyroid tumors and adjacent normal tissues as described (7). Antigen retrieval was done in 1 mmol/L Tris-HCl (pH 8.0) using a steamer for 15 min. Nonspecific protein binding sites were blocked by incubation for 30 min in 1% bovine serum albumin/PBS. Subsequently, slides were incubated overnight with 5 μg/mL primary monoclonal antibody against human BRAK/CXCL14 (clone 131120; R&D Systems). The DAKOCytomation horseradish peroxidase system was used for immunodetection. Sections were incubated with the horseradish peroxidase substrate 3,3′-diaminobenzidine (DAKOCytomation), counterstained with 1% hematoxylin, and permanently mounted. All slides were scored in a blind fashion and the immunopositivity was determined as follows: negative when <10% of the cells were immunoreactive or immunoreactivity was absent and positive when >10% of the cells were immunoreactive.

Statistical analysis. To identify PTC-associated transcripts, qPCR data were log-transformed before application of the ANOVA test with Bonferroni adjustment to keep the family-wise error rate at 0.5. Immunohistochemistry was analyzed using χ². To identify genes that were statistically different between the metastatic and the nonmetastatic groups, qPCR data was log-transformed before application of the Student's t test (unpaired). The test for association between BRAF and gene expression was done using a Student's t test (unpaired). The tests of qPCR comparisons were done using the parametric (Pearson) correlation. P < 0.05 was used as the cutoff for statistical significance.

SAGE analysis and selection of PTC-associated transcripts. We sought to identify PTC-associated genes. Comparing the expression profile between PTC and FTC revealed that 1,670 transcripts were differentially expressed at a Monte Carlo P ≤ 0.05. Using a more stringent criterion (P < 0.001), we generated a list of 548 transcripts. Of these genes, 131 were exclusively expressed in the PTC library by >10-fold and 61 were expressed by >15-fold. Transcripts that were overexpressed in PTC by >15-fold (excluding tags matched to unknown transcripts and internal tags) were assessed in the FTA and normal thyroid libraries. The most highly expressed transcripts in PTC, not expressed in FTC and not expressed or at low levels in FTA and normal thyroid libraries, are described in Table 1.

As expected, our analysis revealed transcripts that were identified previously as overexpressed in PTC but also provided novel information (Table 1). Based on the higher difference on level of expression, five genes (CST6, CXCL14, BCAN, SPP1, and DRHS3) were chosen for transcript level validation by qPCR (Table 1).

Validation of overexpressed genes in PTC by qPCR. A confirmatory qPCR analysis was done for selected PTC-associated genes in an independent set of samples. Relative expression (RE) value for each gene in each profiled sample is presented as mean of log-transformed ratios (Fig. 1). Briefly, CST6 expression levels were increased in 64% of PTC but not expressed in FTA, most normal thyroid tissues, and FTC (Fig. 1A). CXCL14 and DHRS3 expression was elevated in 68% and 79% of PTC, respectively (Fig. 1B and C). CXCL14 and DHRS3 expression was not observed in normal adjacent thyroid tissues and most of FTA and FTC. Although BCAN was expressed in 64% of PTC, it was expressed in several FTA and FTC (Fig. 1D), suggesting that it is not specifically associated with tumorigenesis of PTC. SPP1, a described RET/PTC-induced transcriptional target gene (10), was expressed in ∼30% of PTC (Fig. 1E). Interestingly, SPP1 was found expressed at similar levels in ∼24% of FTC. This is the first analysis of SPP1 expression in FTC.

To test whether the difference of expression was significant, the log-transformed RE data were used for statistical analysis. This analysis indicates that CST6, CXCL14, and DHRS3 were differentially expressed between PTC and FTC at significant levels (P < 0.0001, ANOVA test with Bonferroni correction). The expression was also considered significant among PTC, normal thyroid tissue, and FTA (P ≤ 0.0025). The results are summarized in Fig. 1. BCAN was differentially expressed at significant levels between PTC and normal thyroid tissues (P = 0.0014). Although SPP1 was expressed at higher levels in PTC, it was not considered significant.

Protein expression. The chemokine CXCL14 has been described as ubiquitously expressed in normal tissues, inflammatory cells, and other stromal cells in multiple cancer types (11–14). To confirm the qPCR data and investigate the immunolocalization of CXCL14, immunohistochemical analysis was done in an independent set of 71 thyroid tumors and adjacent normal thyroid. CXCL14 staining was observed in the cytoplasm of tumor epithelial cells in 30 PTC (71%), whereas adjacent nonneoplastic cells show no expression of this chemokine (Fig. 2A-D). In addition, most follicular adenomas and follicular carcinomas were negative for CXCL14 in both tumor cells and normal adjacent tissue (Fig. 2E and F; P < 0.0001).

Evaluation of CST6, CXCL14, DHRS3, and SPP1 expression in metastatic and nonmetastatic PTC. Because CST6, CXCL14, DHRS3, and SPP1 have been associated previously with progression of other tumor subtypes (13, 15, 16), we sought to investigate whether the expression of these genes could predict lymph node involvement. The expression of these transcripts was assessed in a set of 54 metastatic and 76 nonmetastatic PTC. CST6 and CXCL14 tumor expression levels were higher in metastatic than in nonmetastatic PTC (P < 0.0001 and P = 0.0037 respectively; Fig. 3A and B). Although DHRS3 was expressed in both groups, its expression was higher in nonmetastatic PTC (P = 0.0261; Fig. 3C). Although SPP1 had a higher expression in metastatic than nonmetastatic PTC samples, the difference was significant when we considered only PTC with >1 cm (P = 0.0058; Fig. 3D).

Protein-protein association. To determine possible associations among the PTC overexpressed proteins or to find a pathway association, we used a curated signal transduction pathway resource frequently used to analyze the molecular networks.⁷

This human pathway annotation project revealed that these molecules are not involved in any of 20 pathways listed on this resource.

Correlation between BRAF^V600E mutation and gene expression. The BRAF^V600E mutation is a common event found in the metastatic and nonmetastatic PTC samples used in this study.⁸

Because targets of this pathway are still unknown, combining expression analysis with BRAF status may reveal genes regulated by the RAS/RAF/MEK/extracellular signal-regulated kinase/mitogen-activated protein kinase signal transduction pathway. Interestingly, the expression of CST6 and CXCL14 differed significantly between the groups BRAF mutated and BRAF wild-type (P < 0.0001; Fig. 4A and B). Moreover, when we examined a possible correlation among the four genes identified as overexpressed in PTC, CST6 and CXCL14 were correlated (r = 0.745; P < 0.0001). No additional correlations were found. These findings suggest that CST6 and CXCL14 expressions may be induced subsequently to BRAF activation and therefore may be downstream of the BRAF/MEK/extracellular signal-regulated kinase signaling pathway.

Based on our findings, we tested the level of expression of CXCL14 in thyroid carcinoma cell lines with known BRAF mutation status. The qPCR results showed an ∼30-fold increased expression of CXCL14 in NPA (homozygous BRAF mutation) compared with WRO (homozygous wild-type BRAF; Fig. 5).

To identify PTC-associated genes, a PTC SAGE library was compared with three SAGE libraries generated from different thyroid tumor subtypes and normal thyroid tissue. SAGE was used because it is a comprehensive and quantitative method for analyzing transcriptomes within a cell population.

Our profiling analysis led to the identification of tags that were most abundant in PTC and were not associated previously with thyroid tumorigenesis process such as DHRS3, NRGN, NUM, GFAP, and BCAN. Moreover, we identified transcripts that were described previously as overexpressed in PTC, such as CST6, CXCL14, SPP1, and LGALS3 (4, 6, 16–18). CXCL14, DHRS3, BCAN, SPP1, and CST6 were chosen for validation. Our qPCR results showed that CST6, DHRS3, and CXCL14 were overexpressed in PTC at significant levels compared with all tumor subtypes. Although BCAN was overexpressed in several PTC, it was also expressed in several FTC and FTA. We then compared the expression of BCAN in PTC with normal thyroid tissues and found a significant difference.

CST6, also named cystatin M, is a secreted inhibitor of lysosomal cysteine proteases, similar in class to cathepsin B. Interestingly, CST6 was first described as a tumor suppressor gene that is epigenetically silenced as a consequence of promoter hypermethylation in metastatic breast cancer cell lines (19–21). Others have found CST6 markedly elevated in metastatic oropharyngeal squamous carcinoma cell lines and in primary tumor from which these cell lines were originated (15). They suggested that CST6 expression aids metastasis by blocking intrinsic cathepsin B and rescuing tumor cells from tumor necrosis factor–induced apoptosis (15). Although the role of overexpression of CST6 in PTC is still unclear, our findings are supported by the fact that CST6 was listed previously as one of the genes up-regulated in PTC and not expressed in FTC (17, 18).

DHRS3, also named retSDR1, encodes an enzyme that is involved in the generation of storage forms of retinol. The gene is located in 1p36.1, a region commonly deleted or rearranged in human cancers (22). Although some authors suggest that DHRS3 is involved in a growth/tumor-suppressive pathway, DHRS3 was found constitutively expressed in breast cancer cell lines that were capable of storage of retinol through esterification (22, 23). As predicted by SAGE, DHRS3 was expressed in normal cells, but its expression was higher in PTC. A better understanding of the properties and the physiologic role of DHRS3 in PTC need to be clarified.

CXCL14, also named BRAK, is a new CXC chemokine with unknown function and receptor selectivity (24). CXCL14 was described as constitutively expressed in normal tissues, whereas it is not expressed in the majority of head and neck squamous cell carcinoma, some cervical squamous cell carcinoma, and many established tumor cell lines and human cancers (14, 24). It was also shown that CXCL14 is overexpressed in the inflammatory cells adjacent to the tumors (11, 12, 14). Recently, CXCL14 was found as overexpressed in myoepithelial cells of in situ breast carcinoma (13). It was suggested that CXCL14 binds to its receptor on epithelial cells and enhances their proliferation, migration, and invasion and therefore acting as a paracrine factor. Although CXCL14 expression was found in myoepithelial cells of in situ breast carcinoma, in invasive breast carcinoma cells, the expression of CXCL14 was restricted to the tumor epithelial cells. The authors suggested that, during tumor progression, a paracrine factor may be converted into an autocrine one due to its up-regulation in the tumor epithelial cells (13).

In our study, CXCL14 transcript was markedly higher in PTC than in normal thyroid and thyroid lesions analyzed; this is the first report validating CXCL14 expression in PTC. Additionally, immunohistochemical analysis showed that its expression was specifically found in epithelial cancer cells. Although at this moment we can only guess, CXCL14, similar to other chemokines, is also likely to signal through a G protein–coupled receptor and contributes to progression of PTC.

SPP1 was first described in thyroid as a RET/PTC-induced transcriptional target gene (10). Later, it was shown that SPP1 is a feature of PTC cells with both RET/PTC or BRAF^V600E mutation (4). The authors showed that SPP1 prevalence and intensity of staining positively correlated with presence of lymph node metastases and tumor size (>1 cm). A low prevalence of SPP1 overexpression was found in follicular variant of PTC (16).

SPP1 was also significantly overexpressed in a variety of cancers, and recent studies strongly support the crucial role of SPP1 in tumor progression and metastasis (25). In this study, we assessed the expression of SPP1 using qPCR. Although SPP1 was overexpressed in PTC, the difference was not considered significant. Of note, the validation set included not only PTC classic variant not larger than 1 cm but also the follicular variant of PTC where SPP1 was found to be expressed at lower levels. Interestingly, SPP1 was found expressed at the same levels in FTC; this is the first analysis of SPP1 in FTC.

Most of the genes validated here were correlated previously with tumor progression (13, 15, 16). To determine whether CST6, DHRS3, CXCL14, and SPP1 could play a role in determining lymph node involvement, we analyzed their gene expression changes in a second set of samples that included metastatic and nonmetastatic PTC. CST6 and CXCL14 were associated with presence of lymph node metastases. SPP1 expression, however, was associated with presence of lymph node metastasis only for tumors larger than 1 cm. DHRS3 is expressed at lower levels in metastatic PTC than in nonmetastatic primary tumors. These findings corroborate with previous important studies, where it was suggested that there is a gene expression pattern associated with metastases that seems to be present in primary tumors and that this molecular program of metastasis is shared by multiple solid tumors (26, 27).

A common finding for the majority of PTC is activation of the RET/RAS/RAF/MEK/extracellular signal-regulated kinase pathway, which signals for cell growth and proliferation and prevents apoptosis and differentiation (1). Transcriptional profiling of either in vitro models or cohort of PTC harboring mutation in one of these oncogenes has produced some contradictory conclusions concerning gene expression altered by this pathway. In vitro models showed that RET/PTC3, HRAS, and BRAF induce changes in the expression of a widely overlapping set of genes; it is triggered by up-regulation of several CXC chemokines and their receptors (4). Other found that the mutational status is the principal determinant of the gene expression variation observed within PTC (6). One important aspect that emerged from this last study is the fact that these oncogenes may be able to signal through different pathways and therefore activate different targets. Additionally, different proteins may be coupled at different levels to this cascade. Therefore, a more complete understanding of the gene expression induced by the RET/RAS/BRAF/mitogen-activated protein kinase oncogenic cascade would be useful.

To investigate if there is an association between the presence of BRAF mutation and induction of transcripts expression, we compared the BRAF status in a set of metastatic and nonmetastatic PTC with gene expression changes. We found that CST6 and CXCL14 overexpression was positively correlated with the BRAF^V600E mutation in PTC. Our data are in agreement with one of the previously mentioned studies where there is a mutation-specific expression signature that reveals genes with expression altered by the mutation (6). Giordano et al. found 17-fold increased expression of CST6 and 14-fold increased expression of CXCL14 (mean expression values) in a set of BRAF mutant PTC compared with RAS mutants PTC (6). Additionally, their expression profile suggested that BRAF mutation has a unique role in initiating an immune response in PTC and that chemokines such as CXCL14 are involved in the regulation of inflammation status.

Lastly, we also show that CXCL14 was overexpressed in a thyroid carcinoma cell line with a homozygous BRAF mutation, whereas its expression was absent in the BRAF wild-type WRO. These results are consistent with the induction of CXCL14 expression by BRAF mutations in PTC.

In this study, we identified and validated new genes associated with the etiology of PTC. Indeed, CST6 and CXCL14 are associated with the presence of lymph node metastases and DHRS3 is underexpressed in metastatic PTC. Consequently, these transcripts might be prognostic markers for these tumors, but this hypothesis needs testing. Furthermore, CST6 and CXCL14 may be downstream targets of BRAF pathway mutational activation.

No potential conflicts of interest were disclosed.