Abstract
Background/Aim: Extrathyroidal extension (ETE) indicates the invasiveness of primary thyroid tumor into the adjacent tissue, and its importance as a prognostic factor overrides tumor size in classifying the cancer stage. The aim of this study was to determine the molecular basis of ETE in papillary thyroid carcinomas (PTCs). Materials and Methods: We systematically defined genes and pathways regulated in ETE using mutation and gene expression profiles from The Cancer Genome Atlas, and examined the effect of BRAF and RAS mutations on ETE. The significance of these genes was further validated using public microarray data. Results: Genes related to extracellular matrix and immune response were significantly up-regulated in ETE and ion-transport genes were often down-regulated. Differentiation properties and WNT signaling were also found to be altered by ETE. BRAF and RAS mutations were shown to have distinct effects on genes associated with ETE. Specifically, PAX8 and its downstream target WNT4 were differentially expressed according to mutation status in addition to ETE, indicating their critical roles in cell motility. Conclusion: BRAF V600E mutation predisposes PTC cells toward invasive phenotypes, while RAS mutation confers resistance to ETE. The differential regulation appears to be mediated through WNT4.
Abbreviations: ANOVA: Analysis of variance; BRAF: B-Raf proto-oncogene, serine/threonine kinase; CDF: cumulative distribution function; DIO1: deiodinase, iodothyronine, type I; ETE: extrathyroidal extension; FOXE1: forkhead box E1; FV: follicular variant; GO: Gene Ontology; HHEX: hematopoietically expressed homeobox; NCBI: National Center for Biotechnology Information; NKX2-1: NK2 homeobox 1; OCV: oncocytic variant; OR: odds ratio; PAX8: paired box 8; PTC: papillary thyroid carcinoma; RET: ret proto-oncogene; RSEM: RNA-Seq by Expectation Maximization; SFRP: secreted frizzled-related protein; SLC5A5: solute carrier family 5 (sodium/iodide cotransporter), member 5; SLC26A4: solute carrier family 26 (anion exchanger), member 4; TCGA: The Cancer Genome Atlas; TCV: tall cell variant; THRA: thyroid hormone receptor alpha; THRB: thyroid hormone receptor, beta; TPO: thyroid peroxidase; WNT: wingless-type MMTV integration site family.
In recent years, thyroid cancer has received increasing global attention, as its incidence has increased by the fastest rate among all cancer types over the past two decades (1). Papillary thyroid carcinoma (PTC), which is derived from follicular thyroid cells, is the most common type of thyroid cancer and comprises of most newly identified cases. PTC is a differentiated cancers, and has an excellent prognosis, with a low mortality rate (1 per 200,000 population) (2). However, it often spreads to local lymph nodes, which increases the risk of disease recurrence (3). Unlike other types of solid cancer, age is the main prognostic factor for PTC, with the outcome being better for patients under the age of 45 years (4). In addition to age, the most important clinicopathological factor for tumor stage is extrathyroidal extension (ETE), which refers to the invasiveness of the primary tumor into adjacent tissues beyond the thyroid tissue (5). Macroscopic ETE is determined by the surgeon during thyroidectomy, while microscopic ETE can be confirmed by a pathologist through microscopic evaluation of biopsied tissue (6). There has been a marked increase in the rate of detection of microscopic ETE but its clinical implications are controversial (7, 8).
The Cancer Genome Atlas (TCGA) research network has produced a massive amount of cancer genomics data, including data on DNA variants, mRNA and microRNA expression, and methylation, using next-generation sequencing (9). The power of such data is in its enriched clinical annotation, which can improve the application of genomic data in clinical practice. A recent comprehensive analysis of thyroid cancer genomics revealed that molecular subgroups were differentially identified for cancer with BRAF V600E and RAS mutations (10). BRAF and RAS mutations are major driver mutations in thyroid cancer that have fundamentally different consequences (11, 12).
While many studies have investigated the clinical relevance of ETE (13, 14), its molecular properties have not been clarified, mainly due to the lack of clinical annotation of the genomic data. The aim of the present study was to define the genetic components associated with ETE by utilizing the abundance of clinical data in TCGA. The differential effects of main driver mutations (BRAF V600E and RAS mutation) on the development of ETE are described herein. The findings of this study contribute to the development of genetic markers for ETE, promoting their application for personalized diagnosis and treatment.
Materials and Methods
Data retrieval from TCGA. Genomic and clinical data in human PTC were retrieved from TCGA (http://cancergenome.nih.gov). Only patient samples with a primary tumor were used. The BRAF mutation analysis only used the V600E oncogene mutation, which is the most predominant form in PTC; three RAS oncogenes (HRAS, NRAS and KRAS) were included in the RAS mutation study. The accuracy of the pathological classification was improved by having the endocrine surgeon download the raw pathological reports and examine them against the corresponding clinical datasheet.
Validation using public microarray data. Gene regulation by driver mutations in TCGA data was validated using a public microarray data. The dataset contains transcriptional expression profiles of 51 PTCs with their morphology and RET/PTC, BRAF and RAS mutational status (GSE27155) (12). Microarray data normalization was carried out by the robust multi-chip analysis (RMA) method using Affymetrix® Expression Console™ Software (http://media.affymetrix.com/support/downloads/manuals/expression_console_userguide.pdf). We determined whether a gene was expressed using the MAS 5.0 method. Significant genes were selected by fold change greater than 1.2 and p-value of <0.05 (t-test).
Analysis of significant genes. Statistical analyses for differentially expressed genes were performed using Student's t-test or ANOVA test in the R program (15). Normalized values by RNA-Seq by Expectation Maximization (RSEM) method were used for mRNA expression analysis (16). The expression value for each experimental group was summarized using the median of measurements by group. k-means clustering was conducted to classify genes based on their expression pattern by MultiExperiment Viewer (MeV) (http://www.tm4.org/mev.html).
Gene-set analysis. Biological pathways associated with significant genes were analyzed by Gene Ontology (GO) analysis (http://www.geneontology.org). The mapping between genes and GO entries was obtained from the NCBI Gene database (http://www.ncbi.nlm.nih.gov/gene). The association of genes and pathways was assessed using the Fisher's exact test. Redundant concepts in GO results were removed by comparing overlapped portions between genes related to GO entries. The significance analysis of genes associated with a given GO group was also carried out using empirical cumulative distribution function (CDF) of log2 fold change in each dataset. The Kolmogorov-Smirnov test was used to determine if genes in the test set differed from the other genes. Network analysis of genes was carried out in cBioPortal (17).
Results
Evaluation of clinical features associated with ETE. The association between clinical features and ETE was evaluated to determine their roles in the development of ETE. The numbers of patients with PTC with clinical features according to the presence of ETE are given in Table I with the significant associations (assessed by the Fisher's exact test). Lymph node metastasis was found to occur more often in the presence of ETE [odds ratio (OR)=3.8, p=5.2×10−11]. Notably, ETE was associated with the histological type of PTC (p=9.5×10−09), being more common in the tall cell variant (TCV), which usually has an aggressive course, and comparably rare in follicular variant and oncocytic variant. The frequent occurrence of ETE in TCV is consistent with previous reports (18). The occurrence of ETE was positively correlated to BRAF V600E mutation (OR=2.5, p=3.8×10−06) but negatively correlated to RAS mutation including NRAS, HRAS, and KRAS (OR=0.3, p=1.2×10−03). The prevalence of ETE was not affected by maximal tumor size and recurrence.
The effect of clinical features on the gene expression response to ETE was evaluated. The patients were stratified into two groups according to the presence of ETE (i.e. ETE and no ETE). The relative expression changes according to ETE was quantified by log2 fold change of medians in the two sub-groups for each feature. The divergence of the gene-expression changes across clinical subgroups was measured by displaying Pearson correlation coefficients in a heatmap (Figure 1). The correlations tended to be strong across various clinical settings, suggesting that the genes associated with ETE were generally similar among clinically defined sub-groups. However, it is notable that RAS mutation was poorly correlated with other features (median correlation coefficient=0.06), indicating that samples with the RAS mutation were uniquely associated with ETE.
Genes correlated by the degree of ETE reflect significantly different regulation of ion transport and the tumor microenvironment. The genetic indication of the severity of ETE was identified by first excluding patient samples with RAS mutation (n=51) who exhibited a distinctive expression associated with ETE, and then grouping the remaining samples into three groups according to TCGA pathological data: no ETE (n=289), minimal ETE (n=125), and moderate/advanced ETE (n=21). ANOVA was used to select genes significantly differentially expressed among the three groups, and then Tukey's test was used to uncover which specific group pairs made a difference. In total, 326 genes were selected in association with ETE (relative fold change >2 and ANOVA p<0.00001), and their mRNA expression values were subject to k-means clustering. The results are shown in Figure 2A, where cluster 1 represents genes down-regulated with ETE, and cluster 2 and 3 represent those up-regulated with ETE.
In GO analysis, down-regulated genes (cluster 1) tended to change progressively with the extent of ETE, and were often involved in ion transmembrane transport. In cluster 2, genes associated with the immune response and the response to external stimuli were up-regulated with minimal ETE, but no further up-regulation was observed with the advance of ETE. The genes in cluster 3 were progressively up-regulated according to the degree of ETE. This cluster was highly enriched in genes related to extracellular matrix organization. The two types of patterns of up-regulation suggest that an activated immune response is more important for initiating and sustaining ETE, while extracellular matrix remodeling is related more to the progression of ETE. The regulation profiles of gene sets were also compared based on CDF (Figure 2B). The difference in two CDF curves represents the extent of differential regulation. Genes in “ion transmembrane transport” tended to be mildly down-regulated, while those in “extracellular matrix organization” tended to be highly up-regulated with ETE.
Network analysis of the interaction between the genes found in each cluster was conducted using cBioPortal (17) (Figure 3). Genes in clusters are marked by nodes with a thick border. In the network of cluster 1 and 2, genes were indirectly connected. Conversely, the network of cluster 3 exhibited strong association of genes encoding components of the extracellular matrix, such as fibronectin (FN1), decorin (DCN), collagens (e.g. COL1A1 and COL1A2), and thrombospondin 1 (THBS1). This suggested that ETE is associated with a coordinated regulation of genes involved in extracellular matrix organization.
ETE is associated with loss of differentiation and altered WNT signaling. Impairment of the ion-transport system is known to increase cell motility and invasiveness of PTC (19). In this category, two genes known to be important for iodide transport were also found to be down-regulated with ETE: SLC5A5 and SLC26A4 (Figure 4A). SLC5A5 encodes sodium-iodide symporter for iodide uptake in the basolateral membrane and SLC26A4 encodes pendrin (a Cl−/I− transporter) for iodide efflux across the apical membrane (20). Genes involved in the process of hormone biosynthesis were also down-regulated, reflecting the reduced production of thyroid hormone in PTC.
Since the level of cellular differentiation is associated with tumor progression, the genes involved in the regulation of thyroid function were evaluated, as previously seen in a study of TCGA thyroid cancer (Figure 4B) (10). During ETE progression, a significant decrease in the mRNA expression of many genes was observed, indicating that PTC undergoes loss of function with ETE progression. In particular, genes involved in thyroid hormone production (e.g. DIO1 and TPO) tended to be down-regulated, but those encoding thyroid hormone receptors (e.g. THRA, THRB) were unchanged. Since genes encoding thyroid transcription factors regulate expression of genes involved in thyroid function, their mRNA regulation was also examined according to ETE (21) (gene names marked in green in Figure 4B). NKX2-1, FOXE1, PAX8 and HHEX were found to be expressed together in epithelial thyroid follicular cells. The mRNA of thyroid transcription factors tended to be down-regulated in the progression of ETE even though the extent of regulation was relatively small (log2 fold change <0.5).
The aberrant expression of WNT signaling components can predict patient outcome but the relationship is often context-dependent (22, 23). In the present study, genes related to the WNT signaling pathway were often de-regulated in association with ETE (Figure 4C). In particular, genes encoding proteins secreted in the extracellular region tended to be highly de-regulated. Among them, those encoding SFRP2, SFRP4, WNT7A, and WNT2 were significantly up-regulated in ETE but the gene encoding WNT4 was down-regulated. SFRP4 is reportedly a WNT4 antagonist, by binding directly to WNT4 proteins and preventing their interactions with Frizzled receptors in ovarian biology (24). Although the role of SFRP4 in PTC is still not clear, the opposing regulation of SFRP4 and WNT4 appears to coordinate an ETE-promoting response.
BRAF and RAS mutations have distinct effects on genes relevant to ETE. Point mutations in BRAF and RAS genes are common in PTC and are mutually exclusive. Both genes are active in the mitogen-activated protein kinase (MAPK) signaling pathway but the phenotypic consequences of the two oncogenic mutations are quite different in PTC (10). As demonstrated by the data in Table I, RAS mutation was related to a lower incidence of ETE, while BRAF V600E was more likely with ETE. A recent TCGA genomic study of PTC found that RAS-like tumors are more differentiated than are BRAF V600E-like tumors (10). Therefore, whether the two mutations can affect the incidence of ETE by regulating thyroid cell differentiation was herein assessed.
To examine whether these PTC driver mutations render cells susceptible or resistant to ETE, the relative gene expressions according to ETE (ETE vs. no ETE) and the two mutations (BRAF vs. no BRAF, RAS vs. no RAS) were compared (Figure 5A). The expression changes according to mutation were strongly correlated with ETE, but the direction of the regulation varied with mutation type. In addition, many genes deregulated in association with ETE also tended to be re-gulated by BRAF and RAS mutation in opposing directions (Figure 5B). The relevant GO terms associated with the shared genes were very similar to those associated with ETE (Figure 2). The relationship and GO entries were also confirmed on microarray data even though number of genes shared between ETE and mutation types was reduced (Figure 5C).
Since genes encoding thyroid transcription factors tended to be down-regulated with ETE (Figure 4B), their expression levels were examined by mutation status (Figure 6A). Notably, their basal levels were low in tumors with BRAF mutation compared to those with RAS mutation. This suggests that two major mutations determine the differentiation level of PTC by regulating mRNA expression of transcription regulators. Their differential regulation was reproducible on a public microarray data (GSE27155) (Figure 6B).
In the direct comparison between BRAF and RAS mutations, WNT4 was one of the most significantly de-regulated genes of the WNT pathway (Figure 7A), being the most down-regulated gene associated with ETE of the WNT pathway (Figure 4C) and oppositely regulated by two driver mutations (red dots in Figure 5A). The basal level of WNT4 was highly reduced in tumors with BRAF mutation but was maintained in those with RAS mutation compared to normal samples (Figure 7B). The regulation pattern of WNT4 according to the two mutations was consistent with microarray data (Figure 7C). WNT4 has recently emerged as an important gene for maintaining the epithelial phenotype of thyroid cells, as evidenced by its down-regulation in human anaplastic thyroid carcinomas (25). BRAF and RAS mutations appear to act as negative and positive regulators of WNT4, respectively. The expression level of WNT4 may be a molecular signal reflecting the extent of ETE progression.
Discussion
ETE has been identified as an important pathological feature reflecting adverse outcomes, but there are relatively few studies on its molecular mechanism. In this study, significant genetic components of PTC in association with ETE were identified by combined analysis of genomic and clinical data. The genes up-regulated in association with ETE involve extracellular matrix remodeling and activation of the immune response, while the down-regulated genes indicate impairment of the ion-transport system in ETE. Mutation status was shown to influence genes associated with ETE by regulating molecular signaling leading to ETE. The BRAF and RAS mutations were found to be positively and negatively correlated to the development of ETE, respectively. The association between BRAF mutation and ETE has been reported elsewhere (26, 27), but as far as we are aware, this is the first study to describe the impact of RAS mutation in the development of ETE.
The deregulation of thyroid transcription factors leads to failure to maintain the differentiated state of the thyroid cells and is often linked to the development of cancer (28). In the present study, all the genes encoding thyroid transcription factors were significantly down-regulated in the presence of ETE and BRAF V600E mutation but their expressions were not reduced in association with RAS mutation (Figure 6). It was recently reported that the level of differentiation was different between tumors with BRAF and those with RAS mutation (10). Our study further showed that the two driver mutations differentially regulate the mRNA level of thyroid transcription factors, resulting in different pathological phenotypes including ETE. The distinct effect of the mutations on transcriptional upstream regulators may provide clues on controlling the transcription program as a therapeutic application.
The reduced expression of WNT4 was found in association with ETE and BRAF V600E mutation, but not RAS mutation (Figure 4C and Figure 7). WNT4 repression has been reported to correlate to malignant cancer types and is involved in the alteration of cell motility and invasion (29-31). A previous study showed that WNT4 mRNA expression is very low in anaplastic thyroid cancer but its expression is highly variable in PTC samples (25). However, the current study also identified that WNT4 expression in PTC is dependent on driver mutation types, suggesting the critical role of WNT4 in PTC. Recently, WNT4 was discovered to be a transcriptional target of PAX8 (32), which was also down-regulated by BRAF mutation in this study. Therefore, a molecular scheme is proposed with PAX8 and WNT4 as downstream targets of oncogenic BRAF signaling, thus affecting cell movement (Figure 7D). This model conflicts with a previous study that found WNT4 activity to be repressed by RAS oncogenic signaling (25). This conflict may indicate the presence of various regulatory modes of WNT4.
In conclusion, BRAF and RAS driver mutations in the activation of MAPK signaling significantly differ in their expression profiles relevant to ETE. BRAF V600E mutation appears to predispose PTC cells to ETE, while RAS mutation appears to protect against the development of ETE. Their differential consequence on ETE appears to be mediated, at least in part, by the regulation of cell motility through WNT4. These results can be used to define the molecular characteristics of ETE at an early stage, and enable more personalized treatments of patients with PTC.
Acknowledgements
The Authors thank all members of Seoul National University Biomedical Informatics for their helpful discussion. Special thanks go to Jeremy Liu (Yale University Computer Science) for the careful reading of the manuscript. This work was partly supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (2012-0000994), a grant of the Korean Health Technology R&D Project, Ministry of Health and Welfare (HI13C2164), and the ICT R&D program of MSIP/IITP (B0101-15-247, Development of open ICT healing platform using personal health data).
- Received November 13, 2015.
- Revision received December 15, 2015.
- Accepted December 18, 2015.
- Copyright© 2016, International Institute of Anticancer Research (Dr. John G. Delinasios), All rights reserved