Effect of BRAF mutational status on expression profiles in conventional papillary thyroid carcinomas

Background Whereas 40 % to 70 % of papillary thyroid carcinomas (PTCs) are characterized by a BRAF mutation (BRAFmut), unified biomarkers for the genetically heterogeneous group of BRAF wild type (BRAFwt) PTCs are not established yet. Using state-of-the-art technology we compared RNA expression profiles between conventional BRAFwt and BRAFmut PTCs. Methods Microarrays covering 36,079 reference sequences were used to generate whole transcript expression profiles in 11 BRAFwt PTCs including five micro PTCs, 14 BRAFmut PTCs, and 7 normal thyroid specimens. A p-value with a false discovery rate (FDR) < 0.05 and a fold change > 2 were used as a threshold of significance for differential expression. Network and pathway utilities were employed to interpret significance of expression data. BRAF mutational status was established by direct sequencing the hotspot region of exon 15. Results We identified 237 annotated genes that were significantly differentially expressed between BRAFwt and BRAFmut PTCs. Of these, 110 genes were down- and 127 were upregulated in BRAFwt compared to BRAFmut PTCs. A number of molecules involved in thyroid hormone metabolism including thyroid peroxidase (TPO) were differentially expressed between both groups. Among cancer-associated molecules were ERBB3 that was downregulated and ERBB4 that was upregulated in BRAFwt PTCs. Two microRNAs were significantly differentially expressed of which miR492 bears predicted functions relevant to thyroid-specific molecules. The protein kinase A (PKA) and the G protein-coupled receptor pathways were identified as significantly related signaling cascades to the gene set of 237 genes. Furthermore, a network of interacting molecules was predicted on basis of the differentially expressed gene set. Conclusions The expression study focusing on affected genes that are differentially expressed between BRAFwt and BRAFmut conventional PTCs identified a number of molecules which are connected in a network and affect important canonical pathways. The identified gene set adds to our understanding of the tumor biology of BRAFwt and BRAFmut PTCs and contains genes/biomarkers of interest.


Background
Over the last decades the incidence rate of thyroid cancer increases worldwide [1]. In Saudi Arabia, thyroid carcinoma (TC) is considered the second most common cancer in young women [2]. About 80% of all TCs are PTCs. The majority of PTCs are histologically classified as conventional PTCs. The follicular variant of PTC (FVPTC) represents the largest subtype and accounts for about 30% of all PTCs [3]. Minor and rare subtypes include Hurthle cell variant PTC and insular PTC which bears an aggressive clinical behavior [4]. Conventional PTCs are characterized on the molecular level by a moderate to high frequency of BRAF mutations (40 % -70 %) that distinguishes them from FVPTCs (10 % -20 %) [5].
BRAF is a cytoplasmic receptor serine/threonine kinase and a key molecule in the mitogen activated protein kinase (MAPK) pathway. BRAF is mutated in diverse human malignancies although frequency and clinical presentation varies considerably between different types of cancers [6]. Over 90 % of all BRAF mutations are a valine by glycine substitution at codon 600 (V600E) in exon 15. Other BRAF mutations affect commonly codons adjunct to codon 600. Although the impact of BRAF mutations in PTC is controversially discussed, many studies found an association of BRAF mut PTCs with unfavorable clinical features including larger tumor size, advanced tumor stage, vessel invasion, capsular invasion, tumor extension, higher risk for lymph node (LN) involvement, distant metastasis, and poor prognosis [7][8][9][10]. Consistent with this, patients with a BRAF mut PTC are considerably older than those with a BRAF wt PTC [5,10]. Whereas a BRAF mutation represents are valuable target for molecular therapy in advanced solid tumors including PTCs, molecular profiles of BRAF wt PTCs are less known and genetic screening for valuable target genes is primarily limited to research studies [11]. The major deregulated key genes in the BRAF wt group are RET and RAS. RAS consists of the highly related genes for HRAS, KRAS, and NRAS. Within the MAPK pathway the RAS molecules transmit signals to the downstream target BRAF. The most common RET/PTC fusions are paracentric fusions with the gene encoding coiled-coil domain containing 6 (CCDC6) contributing to~80 % and with the nuclear receptor coactivator 4 (NCOA4) contributing to~10 % of all known RET/PTC rearrangements. The frequency of RAS mutations and RET/PTC rearrangements differs between the populations studied and depends in part on the inclusion/exclusion criteria for the different histological PTC subtypes [12]. RAS mutations and RET rearrangements are unlikely to act as molecular drivers for onset of malignancy as they are already present in benign thyroid neoplasms [13,14]. This distinguishes them from BRAF mutations which are virtually absent in precursor lesions of PTCs. Until now only a few studies compared expression profiles between BRAF wt and BRAF mut PTCs using array technologies [15][16][17][18]. The relevance of expressional screening in PTC according to their BRAF mutational status is in part related to the different clinical behavior of both PTC groups with the necessity to identify appropriate biomarkers and in part related to the different tumor biology of both groups which is not thoroughly understood [7]. In our study we took advantage of current state-of-the-art technology to screen and analyze a case series of BRAF wt and BRAF mut PTCs for detecting new biomarkers which could become useful to distinguish both groups on the molecular level. We did not include FVPTCs and other smaller histological subtypes of PTC in our screening to minimize expressional bias which might be related to a different histology.

Thyroid samples
We examined 25 specimens from PTCs and seven normal thyroid samples (TN) from patients which were treated surgically in the period between November 2008 and February 2013 at the King Abdulaziz University Hospital, Jeddah, and the King Faisal Specialist Hospital & Research Center (KFSH&RC), Jeddah. In two BRAF mut PTCs, specimens were derived from a recurrence or a local metastasis and in one BRAF wt PTC from an LN metastasis. Normal thyroid specimens were derived from histopathologically unaffected normal thyroid tissue in the course of lobo-or thyroidetomies of thyroid lesions (4 goiters, 1 hyperthyroiditis, 1 PTC, and 1 FVPTC). Histopathological diagnosis and staging of thyroid lesions was performed by an experienced oncologic pathologist (JM) according to established criteria [19,20]. Five BRAF wt PTCs were classified as micro PTCs (≤ 1 cm). This study was approved by the Research Ethics Committee of the King Abdulaziz University, Faculty of Medicine, #358-10, and the Institutional Review Board of the KFSH&RC, #IRB2010-07, and included written informed consent provided by the participants.

DNA extraction and BRAF mutational screening
Genomic DNA from fresh-frozen samples was extracted using the QIAmp DNA tissue kit (Qiagen, Hilden, Germany). DNA concentration was measured with the Nanodrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE). Screening of the BRAF mutational hotspot region of exon 15 was performed as described earlier involving direct sequencing of PCR products spanning the region [5].

RNA sample and array processing
Total RNA was extracted from freshly preserved thyroid tissue specimens using the Qiagen RNeasy Mini Kit (Qiagen, Hilden, Germany) including an on-column DNAse treatment according to manufacturer's recommendations. Quality of the purified RNA was verified on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). RNA integrity number for all evaluated samples was at least 5.0. RNA concentrations were determined using a NanoDrop ND-1000 spectrophotometer. Samples containing each 250 ng of RNA were processed using the Ambion WT Expression Kit (Life Technologies, Austin, TX) and the GeneChip WT Terminal Labeling and Controls Kit (Affymetrix, Santa Clara, CA) according to the manufacturers`recommendations. Fragmentation and endlabeling of samples were monitored by electrophoresis on 3 % agarose gels. Affymetrix GeneChip hybridization, wash and stain kits were utilized in subsequent processing steps. Hybridization mixtures containing each 5500 ng of cDNA were hybridized at 45°C for 17 hrs and 60 rpm to Affymetrix Human Gene 1.0 ST GeneChip arrays. Subsequent to wash and staining at the GeneChip Fluidics Station 450, the arrays were scanned with the GeneChip Scanner 3000 7G. Probe cell intensity data (CEL files) were generated using the GeneChip Command Console Software (AGCC). Human Gene 1.0 ST GeneChip arrays interrogate in total with a set of 764,885 probes 36.079 reference sequences (NCBI build 36).

Gene Expression Analysis
CEL files were imported to Partek Genomics Suite version 6.6 (Partek Inc., MO) and a log-transformed data set of robust multi-array averaged (RMA), background-adjusted, and normalized values was generated. Non-annotated genes and multiple transcripts generated from the same gene were excluded from further analysis. Principal component analysis (PCA) was performed to assess quality as well as overall variance in gene expression between groups of samples. Analysis of Variance (ANOVA) was applied to generate a list of differentially expressed genes using a p-value with a false discovery rate (FDR) (Step up method) < 0.05 and a fold change > 2.0. Two dimensional average linkage hierarchical clustering was performed using Spearman's correlation as a similarity matrix. The generated array data set complies with MIAME [21] and was submitted to NCBI's Gene Expression Omnibus (GEO), accession number GSE54958.

Functional network and pathway analysis
To define molecular networks and canonical pathways in differentially regulated gene sets, pathway analyses were performed by using Ingenuity Pathways Analysis (IPA) software (Ingenuity Systems, Redwood City, CA). For this purpose, statistically differentially expressed genes and their corresponding probe set ID, gene symbol as clone identifier, p-value and fold change were imported into IPA. The program identifies with its functional algorithms those biological functions, interacting drugs and/or diseases that are most significantly related to a differentially expressed gene set. The canonical pathway analysis identifies pathways that are most significantly related to the data set.

Results
Demographic data of patients and histopathological criteria of BRAF wt and BRAF mut PTC are listed in Table 1. A gender shift towards females was observed in the BRAF wt group and mean age was considerable lower in the BRAF wt than in the BRAF mut group (30.9 years vs. 45.9 years). Histopathological criteria including tumor size, LN involvement, tumor focality and tumor stage were comparably more unfavorable in BRAF mut PTCs (Table 1).

Expression BRAF wt vs. BRAF mut PTCs
Employing whole-transcript microarrays (HuGene 1.0 ST) we compared expression profiles of 11 BRAF wt PTCs with 14 BRAF mut PTCs. Seven TN specimens were used as a reference for normal thyroid tissue expression which allowed us to identify differentially expressed genes between both PTC groups and TN samples. Three-D presentation of the PCA displays clustering of BRAF wt PTCs, BRAF mut PTCs and TN samples ( Figure 1). We identified, after excluding nonannotated genes and multiple transcriptional isoforms, 237 annotated genes that were significantly differentially expressed (p-value with FDR <0.05 and a fold change >2) between BRAF wt and BRAF mut PTCs (Additional file 1). Of these genes, 127 were up-, and 110 were downregulated in BRAF wt compared to BRAF mut PTCs. The most significantly upregulated genes in BRAF wt were inositol 1,4,5-triphosphate receptor, type 1 (ITPR1), hepatic leukemia factor (HLF), potassium voltage-gated channel, shaker-related subfamily, beta member 1

Expression PTCs vs. TN samples
Comparison of genes differentially expressed between BRAF wt and TN samples revealed a set of 8249 genes and between BRAF mut and TN samples a set of 8836 genes. To identify genes of interests, e.g. conceivable immunohistochemistry markers, we selected those genes which are differentially expressed between all three groups ( Table 2). In sum, 32 genes met the statistical criteria and the vast majority revealed the comparatively highest expression in BRAF mut PTCs.

Gene networks and canonical pathways
Most significant network functions identified by IPA algorithms and associated with the set of 237 differentially expressed genes were involved in cell signaling, cancer, and cellular development ( Figure 3). This comprehensive network includes 12 molecules which are overexpressed and 16 molecules which are underexpressed in BRAF wt compared to BRAF mut PTCs. The most significantly associated canonical pathway related to the differentially expressed set of 237 genes is the proteinkinase A (PKA) signaling cascade (Additional file 2). The PKA pathway is involved in second messenger signaling and it is stimulated by upstream cascades including the G protein-coupled receptor pathway (Figure 4). The most significant pathways with thyroid specific functions were the thyroid hormone metabolic and thyronamine/iodothyronamine metabolic pathways.

Discussion
We performed one of the first studies using whole transcript microarrays to compare expression profiles solely in conventional BRAF wt and BRAF mut PTCs. Comparison to studies including different histological subtypes of PTCs may result in detecting a lower number of common genes. An expression array study including histological variants of PTC identified, on basis of enhanced stringent threshold criteria, over 80 up-and downregulated genes in the BRAF mut group in comparison to PTCs with either a RAS mutation or a RET/PTC rearrangement [15] and the 40 most up-and downregulated genes have an overlap of~40 % to our list of 237 differentially expressed genes.
Selected differentially expressed molecules DCSTAMP, also known as TM7SF4, has been previously identified as one of the most overexpressed genes in BRAF mut PTC compared to BRAF wt PTC as well in PTC with undetermined mutational status compared to normal thyroid tissue [15,22,23]. It has been supposed that DCSTAMP expression is an immune response related to BRAF mut tumors [15]. DCSTAMP contains signature motifs owned by the family of transmembrane serine proteases and it exhibits degradation activity against extracellular matrix proteins. One of the function of the hepatic leukemia factor (HLF) as a transcription factor is to mediate thyroid hormone activation from the thyroid hormone receptor/retinoid X receptor heterodimer to the hypoxiainducible factor (HIF-1alpha) [24]. The function of C19orf33, also known as H2RSP (Hepatocyte growth factor activator inhibitor type 2-related small protein) is virtually unknown. An enhanced expression with higher levels in LN positive tumors has been observed for C19orf33 in colorectal adenocarcinoma cells at the invasive front [25]. The voltage gated channel molecule KCNAB1 exhibits diverse functions including electrolyte transport, and insulin secretion. Downregulation of KCNAB1 expression has been identified in follicular thyroid carcinomas compared to benign follicular adenomas [26]. ELMO1 functionally interacts with dedicator of cytokinesis 1 (DOCK1) that promotes Rac guanine exchange factor (GEF) activity for Rac proteins of the Rho GTPases family. GTP-loaded Rac proteins initiate downstream pathways that promote cell elongation, migration, and cytoskeleton remodeling [27,28]. The active ELMO1/ DOCK1 complex is anchored via phosphoinositides to the membrane. In our study, the DOCK1 related DOCK5 gene was identified as a 2-fold upregulated molecule in BRAF wt compared to BRAF mut PTCs. The leucine-rich repeat and immunoglobulin-like domain-containing nogo receptorinteracting protein 2 (LINGO2) encodes a single-pass type I membrane protein that is primarily expressed during development in cells adjunct to the epithelial lining of the olfactory pit and in adult brain [29]. In our study, LINGO2 was downregulated in BRAF wt PTCs and even more in BRAF mut PTCs compared to TN samples which may imply a tumor suppressor function for this molecule.
Higher expression levels of the ERK1/2-specific cytoplasmic dual specificity phosphatase 6 (DUSP6) in comparison to benign and normal thyroid cells has been previously associated with PTC, especially with advanced thyroid carcinomas [30,31]. In our study, DUSP6 was 2.1 times higher expressed in BRAF mut than in BRAF wt PTCs assuming that the ERK1/2 related pathway is frequently more utilized in the BRAF mut group. MET overexpression in thyroid cancer has been identified in a number of studies and this molecule was 2.3-fold higher expressed in our study in BRAFmut compared to BRAF wt PTCs which is in accordance with another survey [15] ( Table 2).

Second messenger molecules
ITPR1 is an intracellular receptor for inositol 1,4,5-trisphosphate (IP3) and implicated in the thyroid hormone synthesis pathway. The receptor mediates calcium release from the endoplasmic reticulum upon stimulation by IP3.
Downregulation of ITPR1 has been demonstrated in thyroid cancer in comparison to non-malignant thyroid tissue in a number of studies [32][33][34][35]. Among the members of the phospholipase C (PLC) family, PLCH1 was upregulated in BRAF wt PTCs whereas PLCD3 was downregulated. PLC molecules hydrolyze phosphatidylinositol 4,5-bisphosphate to generate the second messenger diacylglycerol and IP3.
Thyroid hormone pathway molecules LAD1 has been identified as an overexpressed gene in BRAF mut thyroid carcinomas compared to those with a RET/PTC rearrangement [18]. Lad1 expression is regulated by the glucocorticoid receptor (GR) and requires for induction the GR coactivators thyroid hormone receptor associated protein 220 (MED1) and the thyroid hormone receptor associated protein 170 (MED14) [36]. Downregulation of TPO has been considered in a number of studies comparing thyroid carcinoma with thyroid carcinomas with benign thyroid tumors or normal thyroid samples [32,33,37,38]. However, downregulation of TPO and the sodium iodine symporter genes as been previously associated with BRAF mut PTCs in comparison to PTCs with a RET/PTC rearrangement [39]. This is in line with observations that BRAF mut tumors are refractory for radioactive iodine ablation due to downregulation of thyroid hormone biosynthesis pathways [40]. TPO was in our study 13-fold downregulated in BRAF mut compared to BRAF wt PTCs. In addition, other genes involved in thyroid hormone biosynthesis pathway including the solute carrier family-5 member-8 (SLC5A8), solute carrier family 26 (anion exchanger), member 4 (SLC26A4), deiodinase, iodothyronine, type I (DIO1) and deiodinase, iodothyronine, type II (DIO2) were downregulated in our set of BRAF mut PTCs.
Other members of the SLC family that were downregulated in BRAF mut PTCs were SLC1A3, SLC4A4, SLC16A2, and SLC26A7 whereas SLC22A31, SLC30A2, and SLC34A2 were downregulated in BRAF wt PTCs.

ERBB3 and ERRB4 genes
Two of the four structurally related receptor tyrosine kinases ErbB, namely ERBB3 and ERBB4, were differentially expressed in the two PTC groups, i.e. ERBB3 was downregulated and ERBB4 upregulated in BRAF wt PTCs (Additional file 1). Of notice, ERBB4 was also lower expressed in BRAF mut PTCs compared to TN samples (p = 0.0035) (data not shown). Oncogenic functions of ERBB3 and ERBB4, which can form heterodimers and signal through the PI3K/AKT signaling pathway have not been elucidated yet in detail in relation to the BRAF mutational status in thyroid cancer [41][42][43]; however, decreased expression of ERBB4 in PTC vs. normal thyroid tissue has been demonstrated in a RT-PCR study whereas ERBB2 and ERBB3 expression was shown to be increased [44]. On the protein level, a tissue microarray study in proliferative thyroid lesions found a correlation of ERBB3 expression with LN metastasis whereas ERBB4 expression correlated with lower tumor stage [45]. A possible link of ERBB3 to apoptosis can be deduced from a functional in vivo study wherein deletion of ERBB3 in mouse intestinal epithelium induced tumor-specific cell death [46].

MicroRNAs
The most significantly downregulated microRNAs in BRAF wt compared to BRAF mut PTCs was mir492 and the most upregulated one was mir32. Overexpression of mir492 has been previously linked to progressive hepatoblastoma and tumorigenesis of retinoblastoma [47,48]. Of notice, mir492 is processed from the KRT19 gene [47] and both are higher expressed. A possible target for miR492 is the 3`UTR of KRT19 [49] which is 4.4-fold higher expressed in BRAF mut PTCs compared to BRAF wt PTCs suggesting an accumulation of the KRT19 mRNA. Another possible target of miR492 is the thyroid hormone receptor-associated protein 3 (THRAP3) mRNA that harbors two predicted target sites for miR492 [49] and that is significantly downregulated in our case series by approximately 3-fold in both, BRAF wt and BRAF mut PTCs compared to TN samples. Upregulation of mir32 in thyroid cancer vs. benign thyroid tumors has been detected in a microarray study; however functional implication of this microRNA in PTC is not known yet [50]. In summary, our microarray expression study provides a detailed overview of differentially expressed genes, networks, and pathways between BRAF wt and BRAF mut PTCs that gain interest for basic molecular genetics and translational studies in PTCs.

Additional material
Additional file 1: The 237 most differentially expressed genes in BRAF wt vs. BRAF mut papillary thyroid carcinomas.
Additional file 2: The canonical PKA pathway is a second messenger cascade and involved in diverse functions as growth, development, metabolism, DNA replication/recombination, DNA repair and cellular organization. A number of molecules including members of PKA, ryanodine receptors (RYR), inositol trisphosphate receptors (IP3R), and lymphoid enhancing factors/T-cell factors (TCF/LEF) are upregulated (red) and number of molecules including members of phospholipases C (PLC), 14-3-3 proteins, and protein tyrosine phosphatases (PTP) are downregulated (green) in BRAF wt compared to BRAF mut PTCs.

Competing interests
The authors declare that they have no competing interests.
Authors' contributions JM, EH, MG, and MHQ made substantial contributions to conception and design of the study. RA, AA, and MA processed expression arrays, performed BRAF mutational analysis, and were involved in data interpretation. KG, FM, and OAH were responsible for surgeries, oversight of clinical databases and contributed to the conception and design of the study. JM performed histological examinations. SK and HJS performed data analysis. HJS had general oversight of the study. HJS, JM, and MHQ interpreted data and drafted the manuscript. All authors read and approved the final manuscript.