Thyroid samples

We studied specimens from six FVPTCs, seven FAs, and seven normal thyroid (TN) samples that were derived from patients who were treated surgically in the period between February 2009 and April 2013 at the King Abdulaziz University Hospital (KAUH), Jeddah, and the King Faisal Specialist Hospital & Research Center (KFSH&RC), Jeddah. Normal thyroid specimens were derived from unaffected normal thyroid samples. Diagnosis was established by an experienced oncologic pathologist (JM) according to established criteria [18, 19]. DNA extraction of all samples and BRAF mutational screening for the 13 tumor samples was performed as described earlier [3]. This study was approved by the ethical review boards of KAUH (no. 358-10) and KFSH&RC (no. IRB2010-07).

RNA and array processing

Total RNA was extracted from fresh tissue specimens using the Qiagen RNeasy Mini Kit (Qiagen, Hilden, Germany). Extraction protocol included an on-column DNAse treatment as recommended by the manufacturer. Quality of purified RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). RNA integrity number for all evaluated samples was at least 5.0. The NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE) was utilized to determine RNA concentration. Samples each consisting of 250 ng 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. In following processing steps the Affymetrix GeneChip Hybridization, Wash and Stain Kit was utilized. The hybridization mixtures containing each 5500 ng of cDNA were hybridized for 17 hrs to Affymetrix Human Gene 1.0 ST GeneChip arrays in a hybridization oven at 45°C under rotation (60 rpm). This array type interrogates with a set of 764,885 probes 36.079 annotated reference sequences (NCBI build 36). Subsequent to wash and staining, the arrays were scanned in the GeneChip Scanner 3000 7G. Probe cell intensity data (CEL files) were generated using the GeneChip Command Console Software (AGCC).

Gene expression and alternative splicing analysis

CEL files were imported to Partek Genomics Suite version 6.6 (Partek Inc., MO) and a log-transformed data set of robust multi-chip averaged (RMA), background-adjusted, and normalized values generated using the RMA defaults settings as configured in the Advanced Import Option Tab. Quality of experiments was assessed on the basis of the QC metrics table and QC graphical report. Principal component analysis (PCA) was performed to assess quality as well as overall variance in gene expression between sample groups. Analysis of Variance (ANOVA) was applied to generate a list of differentially expressed genes using a false discovery rate (FDR) < 0.05 and a fold change cut-off >2. FDR is the proportion of false positives among all positives. In Genomics Suite the step up method was implemented to control the FDR. Two dimensional average linkage hierarchical clustering was generated utilizing Spearman’s correlation as a similarity matrix on the set of differentially expressed genes (FVPTC vs. FA). To identify significant genes that intersect or non-intersect between differentially expressed groups, a Venn diagram was created from the Venn diagram menu tab by selecting the corresponding gene lists. For alternative splicing analysis, ANOVA was applied to identify probes with differentially expressed exons using FDR < 0.05 and a fold change cut-off > 2. The performance of Human Gene 1.0 ST arrays to detect alternative splicing events has been assessed [20]. The Partek Gene Ontology (GO) enrichment tool was used in the gene expression workflow to group significantly expressed genes (FVPTC vs. FA) into functional categories in the domains cellular components, molecular functions, and associated biological processes [21]. The GO enrichment score indicates the level of differential expression of genes in a functional category. The score is determined by utilizing a Fisher`s Exact test to compare within a GO domain the proportion of genes represented in a functional category to the proportion not presented in the functional category. The GO data resources are curated by the GO consortium in a controlled and structured vocabulary. The generated array data set complies with MIAME [22] and has been deposited at the NCBI’s Gene Expression Omnibus (GEO) under accession number GSE54958.

Functional network and pathway analysis

The Ingenuity Pathways Analysis software (IPA) (Ingenuity Systems, Redwood City, CA) was employed to interpret our expression data in more detail. The main features of IPA are to provide biological insights into a data set by (i) mapping data to the Ingenuity knowledge base that contains IPA curated and GO information, (ii) translating data into diseases and biological functions, (iii) creating de novo molecular networks using algorithms based on biological information, (iv) and determining canonical pathways that are related to the data set. IPA results are sorted and ranked based on statistic scores. We initiated IPA by uploading the gene symbols as clone identifiers, p-values, and fold changes of the set of our differentially expressed genes. IPA workflow comprised core and functional analysis and structural categorization of the processed data. Network assembling was based on the connectivity of the differential expressed genes under the assumption the higher the connectivity of a gene the more significant it is. Non-connective areas were filled with genes from the IPA knowledge base. Canonical pathway analysis included an overlay analysis of predefined pathways with the gene list.

Expression of FVPTC vs. FA

Three-D presentation of the PCA demonstrated that FAs and FVPTCs are clustering closely together as indicated by the bounding ellipsoids (Figure 1). We identified 55 transcripts that were differentially expressed between FAs and FVPTCs (p-value with FDR < 0.05 and fold change cut-off > 2) (Table 1). Forty genes were upregulated and 15 genes downregulated in FVPTCs compared to FAs. The most upregulated genes in FVPTCs include GABA acid B receptor, 2 (GABBR2), neuronal cell adhesion molecule (NRCAM), extracellular matrix protein 1 (ECM1), heparan sulfate 6-O-sulfotransferase 2 (HS6ST2), and retinoid X receptor, gamma (RXRG). The most significantly downregulated genes in FVPTCs include interaction protein for cytohesin exchange factors 1 (IPCEF1), G protein-coupled receptor 155 (GPR155), Purkinje cell protein 4 (PCP4), chondroitin sulfate N-acetylgalactosaminyltransferase 1 (CSGALNACT1), and glutamate receptor interacting protein 1 (GRIP1). In the top five genes that were overexpressed in FVPTCs, the fold change (range, 4.09-18.36) was comparably higher than in the top five genes that were overexpressed in FAs (range, 2.47 - 6.48). Other genes, all overexpressed in FVPTCs, with a fold change ≥ 10 include leucine-rich repeat (kinase 2) (LRRK2), odz, odd Oz/ten-m homolog 1 (Drosophila) (ODZ1), lipase, member H (LIPH), and proteolipid protein 1 (PLP1).

Hierarchical cluster analysis displays the diverging expression profiles of the 55 differentially expressed genes between FAs and FVPTCs (Figure 2). TN samples were included to provide expression values of normal thyroid tissue. Comparison of either FVPTCs or FAs with TN samples revealed 12280 and 11706 differentially expressed transcripts, respectively. The number of intersecting and non-intersecting genes between the three sets of differentially expressed genes is projected into a Venn diagram illustrating the molecular relationships between the groups (Additional file 2).

Enrichment score analysis

An enrichment analysis for functional categories in the GO domains cellular component, molecular process, and biological functions was performed on the set of 55 differentially expressed genes (Additional file 3). The most dominant functional categories in the cellular component domain were membrane part (enrichment score 4.90), synapse part (5.73), and membrane (6.16). In the molecular function domain, the prevalent categories consisted of molecular transducer activity (2.72), binding (2.72), and nucleic acid binding transcription factor activity (2.87) whereas in the biological process domain, the dominant categories included localization (5.80), single-organism process (6.03), and developmental process (10.23).

Exon splicing

Statistical exon splicing analysis identified 62 genes for which at least one exon was significantly upregulated and 25 genes for which at least one exon was significantly downregulated in FVPTCs compared to FAs (Additional file 4). This gene list highly overlaps with the list of differentially expressed genes (Table 1) as 52 of the genes were identical in both sets. An example for an exon splicing event is provided for HS6ST2 (Additional file 5).

Network and pathway analysis

Using IPA, a network was assembled based on those connective relationships that were identified in our differentially expressed genes (Figure 3). Upregulated network molecules in FVPTCs include molecules such as ECM1, vitamin D (1,25-dihydroxyvitamin D3) receptor (VDR), NRCAM, and RXRG and downregulated molecules in FVPTCs include SMAD family member 9 (SMAD9), v-kit Hardy-Zuckerman 4 feline sarcoma viral oncogene homolog (KIT), and GRIP1. The most significantly associated canonical pathway related to the differentially expressed gene set was thyroid cancer signaling (p = 5.58x10^-3).

Top genes upregulated in FVPTC

Among the most overexpressed genes in FVPTCs in our study were GABBR2, NRCAM, ECM1, HS6ST2, RXRG, LRRK2, ODZ1, LIPH, and PLP1. GABBR2 is a component of the heterodimeric G protein-coupled receptor for GABA. Overexpression of GABBR2 has been identified in pediatric post-Chernobyl PTCs compared to sporadic adult PTCs [25]. NRCAM is a member of the neuronal cell adhesion molecules. Prominent overexpression of NRCAM has been detected in PTCs compared to matched normal tissue by RT-PCR and immunohistochemistry which led to the suggestion that the molecule is involved in PTC pathogenesis [26]. ECM1 is a soluble molecule that is implicated in angiogenesis and endochondral bone development. Upregulation of ECM1 has been revealed in malignant compared to benign thyroid neoplasms and it has been proven to be useful as a diagnostic marker for thyroid malignancy in fine needle aspiration biopsies [27]. HS6ST2 is a membrane and extracellular matrix associated molecule. Gene silencing of HS6ST2 in pancreatic cells resulted in inhibition of cell progression [28]. RXRG is a member of the retinoid X receptors and forms heterodimers with retinoic acid, VDR, and the thyroid hormone receptor. Significant overexpression of RXRG has been detected in PTCs compared normal thyroid tissue and the receptor proved to be useful in a classifier set consisting of 19 genes to accurately discriminate PTCs from normal thyroid tissue [29]. A biostatistics and microarray survey in PTC and papillary renal cell carcinomas detected concomitant overxpression of LRRK2 and the MET proto-oncogene [30]. Knockdown of LRRK2 in cell culture experiments minimized MET phosphorylation and impaired MET downstream signaling to STAT3 and mTOR resulting in cell cycle arrest. LRRK2 was considerably higher expressed in our FVPTCs compared to FAs (fold change, 20.97) underscoring the known relevance of the MET pathway in thyroid malignancies [31]. ODZ1 belongs to a conserved group of signaling molecules acting as transmembrane receptors and with their intracellular domain as transcriptional factors [32]. An oligonucleotide array study detected upregulation of ODZ1 in PTCs compared to matched normal thyroid samples and the results were confirmed by RT-PCR [33]. The PLP1 protein is a major myelin component and a marker for oligodendrocytes. Functional involvement of PLP1 in thyroid cancer is virtually unknown; however, in glioblastoma multiforme, PLP1 expression correlated with the extent of the oligodendroma component [34]. LIPH is a secreted enzyme that hydrolysis phosphatidic acid into free fatty acid and lysophosphatidic acid which is known as a ligand for a number of G protein-coupled receptors and implicated in receptor signaling in cancer [35]. Based on a microarray data survey, LIPH has been recently discovered as a putative biomarker for lung adenocarcinomas and bronchioloalveolar carcinomas [36]. LIPH was also found to be higher expressed in the serum of lung cancer patients.

Top genes downregulated in FVPTC

Among the most downregulated genes in FVPTCs compared to FAs were IPCEF1, GPR155, PCP4, CSGALNACT1, and GRIP1. IPCEF1 has been identified as the C-terminal part of a protein in which CNK3 represents the N-terminal part [37]. In vitro experiments have demonstrated that knockdown of CNK3/IPEF1 leads to inhibition of hepatocyte growth factor induced cell migration. GPR155 is an integral membrane G protein-coupled receptor harboring an auxin efflux carrier domain, a pleckstrin/G protein-interacting region, and a winged helix repressor DNA-binding domain; however, the function of the receptor is barely known [38]. PCP4 codes for a sodium-dependent neurotransmitter. Upregulation of PCP4 was detected in cerebral tissue from rats that were hyperthyroidic through triiodothyronine application [39]. CSGALNACT1 is a critical glycosyltransferase of chondroitin sulfate biosynthesis. Highest RNA expression among different tissue types tested was identified in the thyroid and placenta and statistical network reconstruction on expression data of multiple myeloma identified that high expression of CSGALNACT1 was associated with favorable prognosis [40, 41]. GRIP1, also known as nuclear receptor coactivator 2 (NCoA-2) encodes a transcriptional co-activator of steroid, thyroid, retinoid, and vitamin D receptors including thyroid hormone receptor alpha [42, 43]. GRIP1 assists these nuclear receptors in promoter-driven DNA expression.

Network genes

Network analysis displayed interactions wherein approximately one third of our set of 55 differentially expressed genes participated either in the extracellular space, transmembrane domain or intracellular compartments (Figure 3). Involved in multiple interactions were a number of molecules including ECM1, VDR, SMAD9, KIT, NRCAM, and RXRG. VDR is involved in gene activation through heterodimerization with the retinoid X receptor and subsequent binding to vitamin D responsive elements [44]. The functional role of the retinoid X receptor as a heterodimerization partner of the thyroid hormone receptor in hormone signaling pathways depends on the cellular context [45]. Further investigations should be envisaged to address the functional role of RXRG, VDR, and the thyroid hormone receptor in thyroid cancer, especially in FVPTCs. SMAD9 is a member of the receptor SMAD family that is involved in the canonical bone morphogenetic pathway. SMAD9 is known to interact with high-mobility group AT-hook 2 (HMGA2) for which a higher expression has been detected in a quantitative RT-PCR assay in FVPTCs and PTCs compared to FAs [46]. KIT encodes a receptor tyrosine kinase type III and activation mutations are associated with a number of tumor types; however, decreased KIT expression has been identified in malignant thyroid lesions compared to benign lesions [47] which is consistent with our findings.

Nucleic acid binding factors

IGF2BP2 is a member of the insulin-like growth factor pathway and it has been shown that its expression correlates with tumor progression in head and neck squamous cell carcinoma patients [48]. Knockdown of HGMA2 in embryonic rhabdomyosarcoma cells resulted in severe downregulation of IGF2BP2 and knockdown of either HGMA2 or IGF2BP2 led to reduced protein levels encoded by the neuroblastoma RAS viral (v-ras) oncogene homolog (NRAS) [49]. Furthermore, binding of IGF2BP2 to the 3`UTR of NRAS resulted in increased stability of the NRAS mRNA and enhanced NRAS protein levels. As gain-of-function mutations in NRAS are found in a minority of FVPTCs [5], overexpression of IGF2BP2 would be an intriguing alternative that leads to enhanced NRAS signalling in thyroid cancer with follicular histology. CRABP2 is known as a co-activator of retinoid acid and retinoid acid receptor signaling [50]. Hypermethylation of CpG islands in the CRABP2 promoter was associated with its downregulation in PTC compared to normal thyroid tissue [51]. We observed a higher expression of CRABP2 in FVPTCs compared to FAs but of notice, CRABP2 was lower expressed (fold change, 1.1) in our FAs compared to TN samples. TFCP2L1 is a crystallin, alpha B (CRYAB) gene promoter interacting transcription factor and lower expression of TFCP2L1 and CRYAB was revealed in anaplastic thyroid carcinomas compared to benign goiters [52]. Lower expression of TFCP2L1 was also detected in PTCs compared to matched normal thyroid samples in a combined microarray/RT-PCR study [53]. The zinc finger factor PLAG1 is a transcriptional activator of a number of genes including insulin-like growth factor II and the cytokine-like factor 1[54]. An RT-PCR analysis of a number of candidate genes for differentiating benign and malignant thyroid lesions detected concordant higher expression of PLAG1 and HMGA2 in FVPTCs, and to a lesser extent in FTCs, compared to other thyroid malignancies (PTCs and a Hurthle cell carcinoma) or benign forms of proliferative thyroid lesions [15]. A strong correlation of PLAG1 and HMGA2 expression was also detected in a number of thyroid carcinomas and uterine leiomyomas and cell culture experiments led to the conclusion that HMGA2 is an upstream activator of PLAG1[46]. Expression of zinc finger transcription factor ZFN521 is known to be required for embryonic stem cell differentiation from the epiblast state into neuroectodermal progenitor cells [55]. Furthermore, cell culture assays indicated that enhanced expression of ZFN521 is associated with reduced recognition of tumor cells by natural killer cells [56]. Downregulation of LMO3 in radiation-related thyroid carcinomas in comparison to matched normal thyroid samples has been detected in a microarray/RT-PCR study conducted on individuals who received iodine-131 doses as a consequence of the Chernobyl accident [57]. In contrast, overexpression of LMO3 has been demonstrated as a marker for development and progression of neuroblastoma and furthermore, cell culture experiments indicated that LMO3 is a p53 co-suppressor [58]. In lung adenocarcinomas bearing a NK2 homeobox 1 (NKX2-1) amplification, the oncogenic capacity of LMO3 as a transcription regulator downstream of NKX2-1 has been demonstrated [59]. The neuro-oncological ventral antigen 1 (NOVA1) encodes an alternative splicing factor and is transcriptionally subject to A-to-I editing [60]. High expression of the NOVA1 RNA-binding protein has been found to be associated with unfavorable prognosis in hepatocellular carcinoma [61].

In summary, our expression profile study identified a set of candidate genes which are implicated in molecular mechanisms that differentiate FAs from FVPTCs and thus, add to our knowledge of the tumor biology of these histopathologically related, but clinically distinct entities. Further studies have to address the question to what extent the identified genes are solely discriminators of histopathological interest or have functional implication in a not yet defined transition process from FA to FVPTC.