In silico Analysis of Publicly Available Transcriptomics Data Identifies Putative Prognostic and Therapeutic Molecular Targets for Papillary Thyroid Carcinoma.

Background: Thyroid cancer is the most common endocrine malignancy. However, the molecular mechanism involved in its pathogenesis is not well characterized. Purpose: The objective of this study is to identify key cellular pathways and differentially expressed genes along the thyroid cancer pathogenesis sequence as well as to identify potential prognostic and therapeutic targets.
Methods: Publicly available transcriptomics data comprising a total of 95 samples consisting of 41 normal, 28 non-aggressive and 26 metastatic papillary thyroid carcinoma (PTC) cases were used. Transcriptomics data were normalized and filtered identifying 9394 differentially expressed genes. The genes identified were subjected to pathway analysis using absGSEA identifying PTC related pathways. Three of the genes identified were validated on 508 thyroid cancer biopsies using RNAseq and TNMplot.
Results: Pathway analysis revealed a total of 2193 differential pathways among non-aggressive samples and 1969 among metastatic samples compared to normal tissue. Pathways for non-aggressive PTC include calcium and potassium ion transport, hormone signaling, protein tyrosine phosphatase activity and protein tyrosine kinase activity. Metastatic pathways include growth, apoptosis, activation of MAPK and regulation of serine threonine kinase activity. Genes for non-aggressive are KCNQ1, CACNA1D, KCNN4, BCL2, and PTK2B and metastatic PTC are EGFR, PTK2B, KCNN4 and BCL2. Three of the genes identified were validated using clinical biopsies showing significant overexpression in aggressive compared to non-aggressive PTC; EGFR (p < 0.05), KCNN4 (p < 0.001) and PTK2B (p < 0.001). DrugBank database search identified several FDA approved drug targets including anti-EGFR Vandetanib used to treat thyroid cancer in addition to others that may prove useful in treating PTC.
Conclusion: Transcriptomics analysis identified putative prognostic targets including EGFR, PTK2B, BCL2, KCNQ1, KCNN4 and CACNA1D. EGFR, PTK2B and KCN44 were validated using thyroid cancer clinical biopsies. The drug search identified FDA approved drugs including Vandetanib in addition to others that may prove useful in treating the disease.

Introduction

Thyroid cancer was ranked as the most common endocrine malignancy.1 Globally, thyroid cancer incidence has been on the rise over the past three decades. Between 2006 and 2012, the annual incidence rate was 6.5% in women and 5.4 in men.2,3 In the United States between 2000 and 2009, thyroid cancer incidence rate was the highest among all cancers.4 The mortality rate of thyroid cancer is considered to be low, whilst the reoccurrence and persistence of the disease is still considered high.5

Morphologically, thyroid cancers are classified into different cellular subtypes such as papillary, follicular, medullary and anaplastic. Differentiated papillary thyroid carcinoma (PTC) form is the most common type comprising more than 80% of all thyroid cases as shown in Table 1. Genetic mutations have been associated with PTC.6

Table 1

List of Subtypes of Thyroid Carcinoma and the Current Treatment Provided

Tumor Subtype88	Origin89	% of Other Subtypes90	Survival91	Treatment92,93
Papillary	Follicular thyroid cells	80–90	10-year survival: 74–93%	Total thyroidectomy/131I administration/Thyroid-stimulating hormone suppression with thyroxine
Follicular	Follicular thyroid cells	10–15	10-year survival 43–94%	Total thyroidectomy/131I administration/Thyroid-stimulating hormone suppression with thyroxine
Medullary	Parafollicular thyroid cells- C cells	2–3	65–89%94	Total Thyroidectomy/palliative chemotherapy/ teleradiotherapy and substitutive doses of L-thyroxine95
Anaplastic	Follicular thyroid cells	2–3	4–5 months from diagnosis	Surgery: tracheostomy/Chemotherapy
Follicular Thyroid Adenoma	Follicular thyroid cells	Benign	-	Thyroid lobectomy and isthmusectomy
Poorly differentiated thyroid cancer (PDTC)	Follicular thyroid cells	5–10	-	Surgery, radioactive iodine and/or radiation therapy
Thyroid Primary Lymphoma	Lymphocytes	<1	82%88	Chemotherapy/radiation therapy
Metastasis to Thyroid gland from other organs	Non thyroid cells	<1	-	Total thyroidectomy and substitutive doses of L-thyroxine95

Whilst many genomic mutational screening studies were carried out on thyroid cancer in general and PTC in particular, only few have identified mutated genes that are correlated with progression of PTC including TP53 and KRAS/BRAF7. However, although such studies suggested that thyroid cancer has high degree of intra-tumoral heterogeneity,8 the mutations identified did not provide clear insights into the molecular mechanism of thyroid cancer phenotypes and progression. Thus, for better clinical outcomes, there is a compelling need to actively study alterations in cellular pathways linked to the underlying mechanism of thyroid cancer initiation and progression.

Few transcriptomic analyses were carried out on PTC identifying some of the cellular pathways involved in its pathogenesis9. However, such studies were generally carried out on small number of patients using standard bioinformatics analysis focusing on list of differentially expressed genes. This provided limited insights into the molecular basis of PTC without clear association to diagnostic, prognostic and therapeutic targets.

In this study, we carried out comprehensive and systematic in silico pathway analysis of PTC using in-house bioinformatics pipeline that has shown good ability to identify the transcriptomic profiles and related differentially expressed genes between different subtypes of the same disease.10 The aim of this study is to attempt to identify the key transcriptomic signatures that drive non-aggressive and metastatic PTC as well as using such signature to identify putative drug targets for PTC. Such approach can provide insights into some of the molecular mechanisms involved in PTC progression and facilitate the identification of key prognostic and therapeutic targets that might help provide better ways for patient management of PTC.

Methods

Publicly Available Data Sets for Papillary Thyroid Carcinoma

Discovery Set

In order to identify the cellular pathways and differentially expressed genes related to papillary thyroid carcinoma, PTC gene sets were searched and retrieved from gene expression omnibus (GEO). Datasets inclusive of patient’s matched normal thyroid tissue transcriptome were considered for analysis. In order to eliminate platform bias, the gene sets obtained were from the Affymetrix Human Genome U133 Plus 2.0 Array platform. Three gene sets that met such criteria were downloaded. Those were GSE6004, GSE60542, and GSE3678 (Table 2). In total 95 cases were identified and the raw CEL files corresponding to these gene sets were extracted and further processed for Gene Set Enrichment Analysis (GSEA).

Table 2

List of Gene Sets Included in the Study

S No.	Gene Set ID	Population	Type of Sample
			Normal	Non-Aggressive	Metastatic
1	GSE600496	Ukraine	4	7	7	Discovery set
2	GSE6054297	Belgium and France	30	14	19
3	GSE3678	USA	7	7	0
		Total	41	28	26	Grand Total = 95
4	GSE35570	Ukraine	51	32		Validation set
5	GSE50901	Brazil	4	61
6	GSE129562	South Korea	8	8

Validation Set

In order to validate the pathways and genes identified from the discovery set, an independent validation set was constructed from 3 independent gene sets from different populations; Ukraine GSE35570 with 51 normal and 32 thyroid cancer tissue biopsies, Brazil GSE50901 with 4 matched normal thyroid and tumor samples and 57 unmatched thyroid tumor biopsies and South Korea GSE129562 with 8 matched normal and thyroid tumor samples (Table 2). The analysis for this study was approved by the Research Ethics Committee of University Hospital Sharjah (UHS); the ethical approval number of the study is UHS-HERC- 011-10062019.

Raw Microarray Normalization and Adaptive Filtering

Each Affymetrix microarray consists of > 54,000 probes. The raw CEL files for the 95 PTC patients obtained from the GEO for normal, non-aggressive and metastatic thyroid samples were normalized using in house R script as described previously.10 Briefly, Affymetrix microarray suite 5 (MAS5) and Gene Chip Robust Multiarray Averaging (GCRMA) packages in R software were applied to normalize and remove the background noise. The invariant probes were removed from the transcripts list, and non-specific filtering was performed to obtain the common set of variant probes. Adaptive filtering was carried out using R script. Probes with MAS5 value >50 and coefficient of variation (CV) 10–100% in GCRMA across all cases were generated and intersected to obtain probes with common variant probes set. The filtered probes from all the samples were then mapped to gene list using Broad Institute software ().11 The probes with maximum expression for each gene were chosen as the expression value for the gene. Probes corresponding to housekeeping genes or not assigned to any gene were excluded.

Pathway Analysis Using Gene Set Enrichment Analysis

The mapped gene expression list was subjected to Gene set enrichment analysis (GSEA) to identify the activated and enriched cellular pathways in non-aggressive (NAG) and metastatic papillary thyroid carcinoma (PTC) samples in comparison to normal tissue. Absolute GSEA search was carried out on the expression data using around 20,500 annotated cellular pathways obtained across seven well annotated gene sets C1 to C7 obtained from the Broad Institute’s database (). The significantly activated pathways in different types of PTC samples were selected based on p < 0.05 and FDR < 0.25 as previously described.10,12 The selected pathways were further processed to identify differentially enriched genes between the normal versus non-aggressive and normal versus metastatic PTC cases. This was followed by reducing the set of available genes by identifying the frequency of gene occurrence across differentially activated cellular pathways.

Differential Gene Expression in PTC Samples Compared to Normal Thyroid Tissue

The differential gene expression analysis was carried out using two approaches in order to obtain information based on pathway enrichment as well as microarray gene expression. Firstly, the significantly enriched pathways for each sample set were used to obtain genes occurring frequently across all the enriched pathways using R script as described previously.10 Using statistical analysis, the 95-percentile cut-off was calculated for each sample. Secondly, the differentially expressed genes in both non-aggressive and metastatic samples were obtained by calculating the average expression value across each sample set for each gene and a fold change value based on normal tissue expression was determined. Genes with fold change >1.5 were considered as upregulated and fold change <0.5 as downregulated.

In vitro Validation of the Pathways and Genes Identified by GSEA in Independent Cohort

Metascape Analysis

In order to validate the pathways identified by GSEA, the most frequent genes from non-aggressive and metastatic samples were considered. The commonly occurring genes with high frequency amongst both groups were inputted in the Metascape software ()13 to identify significantly activated cellular pathways.

Drug Bank Database Search

The genes differentially expressed and enriched with high frequency in GSEA in NAG and metastatic PTC were used to search in drug bank database to identify the potential drug targets for papillary thyroid carcinoma. Pharmacoinformatics search using the differentially expressed genes identified as targets to search for matching drug using DrugBank repository14 was carried out. Among these, the approved drugs used to treat thyroid cancer were sorted and novel drug targets were listed for the ones not prescribed.

In order to determine the putative therapeutic targets based on different populations, the most upregulated unique genes from each population were used to search DrugBank for associated drugs.

In vivo Validation from Early and Late Thyroid Cancer Tissue Biopsies

Sample Details

Six well characterized United Arab Emirates (UAE) patients biopsies from early and late thyroid cancer were recruited for the study (Table 3). The formalin fixed paraffin embedded (FFPE) tissue biopsies from those cases were subjected to microdissection to enrich the tumour content followed by RNA extraction using modified Recover All protocol as previously described.10 The transcriptomic analysis for this study was approved by the Research Ethics Committee of University Hospital Sharjah (UHS); the ethical approval number of the study is UHS-HERC- 011-10062019.

Table 3

Patient Characteristics for the Six Biopsies Collected from Thyroid Cancer Patients in UAE

S No	Gender	Age	Nationality	Subtype
1	Female	43	Egyptian	Early Thyroid cancer
2	Male	65	UAE	Early Thyroid cancer
3	Female	60	UAE	Early Thyroid cancer
4	Female	33	Tunisian	Late Thyroid Cancer
5	Male	43	Egyptian	Late Thyroid Cancer
6	Female	33	Philippines	Late Thyroid Cancer

RNA Sequencing

Next Generation Sequencing (NGS) RNAseq was applied to the RNA extracted from the microdissected FFPE thyroid samples using AmpliSeq Whole Transcriptome on S5 System (ThermoFisher) as previously described.15 Briefly, the targeted RNA-seq library was prepared using Ion AmpliSeq Transcriptome Human Gene Expression Kit (Thermo Fisher Scientific) which is designed to profile over 21,000 distinct human RNA targets. The prepared template libraries were then sequenced on the Ion S5 XL Semiconductor sequencer using the Ion 540 Chip (Life Technologies Corporation, Carlsbad, CA).

Bioinformatic Analysis

RNAseq data were analyzed using the Ion Torrent Software Suite version 5.4. Alignment was carried out using the Torrent Mapping Alignment Program (TMAP) optimized for Ion Torrent sequencing data for aligning the raw sequencing reads against reference sequence derived from hg19 (GRCh37) assembly. Differential gene expression (DGE) analysis was performed using R/Bioconductor package DESeq230 with raw read counts from RNA-seq data. Read count genes with less than ten normalized read counts were excluded from further analysis. Differentially expressed genes were selected at significance of p<0.05.

Cross Validation of the Molecular Targets on Large Cohort of Cases

Additional validation for the differentially expressed genes from in silico analysis was performed on a larger independent cohort for thyroid cancer RNA-seq data obtained from The Cancer Genome Atlas Program (TCGA). The cohort comprises of 502 non-metastatic thyroid tumor samples, 8 metastatic cases and 58 normal thyroid tissue. The analysis was carried out using TNM plotter ()16 and Kruskal–Wallis test was used for statistical comparison. p<0.05 was considered to be statistically significant.

Results

Normalization and Filtration of the Transcriptome Data for Papillary Thyroid Carcinoma

The flow chart for the process of normalization and filtration is shown in Figure 1. From the total number of 54,675 probes in the Affymetrix Human Genome U133 Plus 2.0 Array, following MAS5 and GCRMA filter 15,801 probes were extracted. These filtered probes were mapped to 9394 genes in GSEA as described in the methods section.

Figure 1

Flow chart of transcriptomics data normalisation and gene set enrichment analysis

Gene Set Enrichment Analysis Identifies the Activated Cellular Pathways in Non-Aggressive and Metastatic PTC Compared to Normal Tissue

The three groups; normal, non-aggressive and metastatic papillary thyroid cancer samples were processed using absolute GSEA. The differentially activated significant pathways across the three different samples were identified by comparing the cancer samples with normal tissue. Significantly differentially activated pathways were obtained based on p < 0.05 as well as false discovery rate (FDR) < 0.05 cutoffs. The results identified around 1795 significantly differentially activated pathways from the molecular functions and biological processes ontology gene sets (Table 4). The most significantly enriched pathways include transforming growth factor beta receptor binding, phosphatase regulator activity, protein tyrosine phosphatase activity, protein kinase activity and calcium dependent protein kinase activity in normal versus non-aggressive set (Table 5). The complete list of pathways enriched can be seen in .

Table 4

List of Number of Significant Pathways Enriched in Non-Aggressive and Metastatic PTC Compared to Normal Thyroid Tissue in Absolute GSEA

Gene Set Analyzed	Description	Total Number of Pathways	Significant Pathways from Absolute GSEA
			NAG	MET
C2	Curated gene sets eg KEGG REACTOME	6229	447	294
C5.bp	Ontology Gene set: biological processes	7573	860	728
C5.mf	Ontology Gene set: molecular functions	1697	107	100
C6	Oncogenic signature	189	78	117
C7	Immunologic signature	4872	701	730

Table 5

List of the Pathways Activated in Non-Aggressive Samples in Comparison to Normal Thyroid Tissue Analyzed by GSEA

Gene Set	Size	ES	NES	NOM p-val	FDR q-val	FWER p-val	Tag %	Gene %	Signal	glob.p.val
Go_regulation_of_ion_transport	298	0.489	2.09	<0.0001	0.004	0.048	0.292	0.184	0.246	0
Go_positive_regulation_of_nervous_system_development	287	0.478	2.147	<0.0001	0.002	0.026	0.401	0.287	0.295	0
Go_regulation_of_hormone_levels	240	0.517	2.272	<0.0001	0	0.002	0.354	0.209	0.288	0
Go_regulation_of_developmental_growth	179	0.471	1.995	<0.0001	0.008	0.141	0.413	0.3	0.295	0.001
Go_regulation_of_membrane_potential	173	0.542	2.333	<0.0001	0.001	0.001	0.318	0.159	0.273	0
Go_organic_acid_transport	168	0.478	2.118	<0.0001	0.003	0.032	0.327	0.214	0.262	0
Go_intracellular_receptor_signaling_pathway	161	0.416	1.835	<0.0001	0.02	0.449	0.366	0.294	0.263	0.002
Go_hormone_transport	156	0.469	2.056	<0.0001	0.005	0.078	0.308	0.209	0.248	0.001
Go_positive_regulation_of_growth	146	0.481	1.963	<0.0001	0.01	0.189	0.349	0.243	0.269	0.001
Go_regulation_of_blood_circulation	128	0.532	2.138	<0.0001	0.002	0.029	0.352	0.197	0.286	0
Go_peptide_hormone_secretion	128	0.489	2.084	<0.0001	0.004	0.055	0.312	0.197	0.254	0
Go_regulation_of_hormone_secretion	126	0.482	2.088	<0.0001	0.004	0.05	0.317	0.209	0.255	0
Go_insulin_secretion	110	0.502	2.149	<0.0001	0.002	0.026	0.327	0.197	0.266	0
Go_regulation_of_peptide_hormone_secretion	105	0.485	2.074	<0.0001	0.005	0.061	0.314	0.197	0.255	0
Go_cell_substrate_adhesion	236	0.455	1.857	0.002	0.017	0.4	0.297	0.198	0.244	0.001
Go_g_protein_coupled_receptor_signaling_pathway	345	0.477	1.933	0.002	0.011	0.233	0.31	0.203	0.257	0.001
Go_regulation_of_wnt_signaling_pathway	221	0.421	1.73	0.004	0.031	0.674	0.281	0.226	0.222	0.001
Go_transmembrane_receptor_protein_serine_threonine_kinase_signaling	198	0.471	1.846	0.004	0.018	0.424	0.379	0.261	0.286	0.002
Go_response_to_transforming_growth_factor_beta	161	0.437	1.72	0.006	0.032	0.69	0.466	0.352	0.307	0.001
Go_positive_regulation_of_apoptotic_signaling_pathway	130	0.406	1.694	0.008	0.036	0.741	0.308	0.249	0.234	0.001
Go_positive_regulation_of_map_kinase_activity	169	0.451	1.725	0.01	0.032	0.681	0.308	0.216	0.246	0.001
Go_positive_regulation_of_peptidyl_tyrosine_phosphorylation	106	0.487	1.696	0.012	0.036	0.738	0.34	0.199	0.275	0.001
Go_regulation_of_protein_serine_threonine_kinase_activity	330	0.401	1.601	0.012	0.055	0.866	0.264	0.212	0.215	0
Go_cell_cycle_arrest	141	0.374	1.63	0.014	0.049	0.837	0.348	0.331	0.236	0
Go_regulation_of_apoptotic_signaling_pathway	256	0.375	1.605	0.018	0.054	0.861	0.285	0.25	0.22	0
Go_positive_regulation_of_erk1_and_erk2_cascade	109	0.491	1.624	0.028	0.05	0.846	0.404	0.243	0.309	0

Abbreviations: ES, enrichment score; NES, normalized ES; NOM, nominal; FDR, false discovery rate; FWER, family-wise error rate; Tag%, the percentage of gene tags before (for positive ES) of after (for negative ES) the peak in the running enrichment score; gene %, the percentage of genes in the gene list before (for positive ES) of after (for negative ES) the peak in the running enrichment score; GO, gene ontology.

Amongst the normal versus metastatic set, negative regulation of peptide hormone secretion, insulin like growth factor receptor signaling pathway, activation of MAPK activity, regulation of MAPK cascade, regulation of protein serine threonine kinase activity, and transmembrane receptor protein tyrosine kinase signaling pathway were among the significantly enriched pathways (Table 6). Example representation of the output from the gene set analysis for each data set is shown in Figure 2.

Table 6

List of the Pathways Activated in Metastatic Samples in Comparison to Normal Thyroid Tissue Analyzed by GSEA

Gene Set	Size	ES	NES	NOM p-val	FDR q-val	FWER p-val	Tag %	Gene %	Signal	FDR (Median)	glob.p.val
Go_growth	563	0.429	1.893	<0.0001	0.014	0.3	0.329	0.26	0.259	0	0.001
Go_regulation_of_cell_development	525	0.447	1.948	<0.0001	0.01	0.19	0.417	0.31	0.305	0	0.001
Go_positive_regulation_of_transport	515	0.41	1.698	<0.0001	0.037	0.792	0.357	0.297	0.266	0.014	0.001
Go_cation_transport	511	0.439	2.008	<0.0001	0.008	0.1	0.399	0.306	0.293	0	0
Go_ion_transmembrane_transport	510	0.448	2.129	<0.0001	0.005	0.014	0.445	0.335	0.313	0	0.001
Go_g_protein_coupled_receptor_signaling_pathway	345	0.474	1.859	<0.0001	0.017	0.397	0.441	0.305	0.318	0	0.001
Go_cell_cell_signaling_by_wnt	311	0.385	1.635	<0.0001	0.05	0.882	0.334	0.284	0.248	0.023	0.001
Go_anion_transport	307	0.444	1.976	<0.0001	0.009	0.149	0.423	0.305	0.304	0	0.001
Go_regulation_of_ion_transport	298	0.469	1.951	<0.0001	0.01	0.187	0.436	0.306	0.313	0	0.001
Go_response_to_extracellular_stimulus	280	0.404	1.737	<0.0001	0.031	0.712	0.443	0.36	0.292	0.011	0.001
Go_regulation_of_transmembrane_transport	266	0.467	1.963	<0.0001	0.01	0.166	0.466	0.336	0.319	0	0.001
Go_organic_anion_transport	240	0.444	1.948	<0.0001	0.01	0.19	0.438	0.305	0.312	0	0.001
Go_cell_substrate_adhesion	236	0.483	1.872	<0.0001	0.016	0.366	0.39	0.259	0.296	0	0.001
Go_regulation_of_wnt_signaling_pathway	221	0.42	1.684	<0.0001	0.039	0.81	0.353	0.284	0.259	0.016	0.001
Go_positive_regulation_of_neuron_differentiation	215	0.448	1.904	<0.0001	0.013	0.28	0.433	0.319	0.301	0	0.001
Go_regulation_of_ion_transmembrane_transport	211	0.478	1.984	<0.0001	0.009	0.135	0.441	0.304	0.314	0	0.001
Go_canonical_wnt_signaling_pathway	197	0.427	1.689	<0.0001	0.038	0.803	0.365	0.284	0.267	0.015	0.001
Go_negative_regulation_of_cell_development	179	0.456	1.895	<0.0001	0.014	0.299	0.419	0.299	0.3	0	0.001
Go_regulation_of_membrane_potential	173	0.541	2.343	<0.0001	0	0	0.347	0.178	0.29	0	0
Go_regulation_of_cation_transmembrane_transport	165	0.489	1.952	<0.0001	0.01	0.187	0.467	0.304	0.33	0	0.001
Go_intracellular_receptor_signaling_pathway	161	0.418	1.797	<0.0001	0.023	0.555	0.292	0.223	0.231	0.006	0.001
Go_hormone_transport	156	0.46	2	<0.0001	0.007	0.112	0.41	0.307	0.289	0	0.001
Go_positive_regulation_of_growth	146	0.425	1.708	<0.0001	0.035	0.766	0.199	0.113	0.179	0.013	0.001
Go_calcium_ion_transmembrane_transport	153	0.424	1.813	0.002	0.021	0.522	0.366	0.294	0.263	0.005	0.001
Go_regulation_of_protein_localization_to_membrane	133	0.445	1.678	0.002	0.041	0.819	0.429	0.32	0.296	0.016	0.001
Go_transmembrane_receptor_protein_tyrosine_kinase_signaling	452	0.416	1.723	0.002	0.033	0.734	0.358	0.289	0.268	0.012	0.001
Go_positive_regulation_of_protein_serine_threonine_kinase	218	0.461	1.7	0.004	0.036	0.79	0.394	0.289	0.287	0.014	0.001
Go_regulation_of_mapk_cascade	434	0.438	1.708	0.008	0.035	0.767	0.366	0.27	0.28	0.013	0.001
Go_positive_regulation_of_map_kinase_activity	169	0.462	1.647	0.008	0.047	0.871	0.402	0.289	0.291	0.021	0.001
Go_regulation_of_peptidyl_tyrosine_phosphorylation	142	0.483	1.67	0.01	0.043	0.827	0.423	0.287	0.306	0.018	0.001
Go_response_to_wounding	381	0.427	1.655	0.012	0.046	0.861	0.399	0.312	0.286	0.02	0.001
Go_regulation_of_apoptotic_signaling_pathway	256	0.39	1.646	0.018	0.048	0.872	0.387	0.308	0.275	0.021	0.001
Go_extracellular_structure_organization	236	0.505	1.678	0.021	0.041	0.819	0.458	0.27	0.343	0.016	0.001

Abbreviations: ES, enrichment score; NES, normalized ES; NOM, nominal; FDR, false discovery rate; FWER, family-wise error rate; Tag%, the percentage of gene tags before (for positive ES) of after (for negative ES) the peak in the running enrichment score; gene %, the percentage of genes in the gene list before (for positive ES) of after (for negative ES) the peak in the running enrichment score; GO, gene ontology.

Figure 2

Representation of heatmaps and graphs for GSEA for significant pathways with enrichment scores. (A) The result file for normal and non-aggressive dataset is presented here with graph for enrichment score. (B) Graphical representation for the GSEA for normal versus metastatic data

Genes Differentially Expressed Among Non-Aggressive and Metastatic PTC in Comparison to Normal Thyroid Tissue

The enriched pathways from GSEA were subjected to gene frequency cutoff using the 95-percentile as a cut-off. Gene frequency can be defined as the number of times a gene occurs across all the enriched gene component from the significantly activated cellular pathways. This type of analysis showed the value for the frequency for non-aggressive (NAG) to be 13 and metastatic (MET) to be 10. Based on those frequency cutoff values, the number of genes with frequency higher than the cutoff in NAG was 355 and in MET was 280. The top 40 genes based on frequency cutoff were shown in Tables 7 and 8

Table 7

List of the Top 40 Genes Based on Frequency in Normal versus NAG Set

Gene	Frequency	Gene	Frequency
KCNQ1	38	CACNA1A	29
CACNA1D	37	EDN3	29
PTK2B	35	EGFR	29
EDN1	34	KCNAB1	29
SFRP1	33	KCNE4	29
ABAT	32	RYR2	29
KCNJ2	32	KCNE3	27
KCNJ5	32	KCNJ8	27
KCNS3	32	KCNK1	27
ADRA2A	31	KCNMA1	27
ANO1	31	KCNQ3	27
BCL-2	31	SCN4B	27
CACNA2D2	31	ADORA1	26
FKBP1B	31	CXCL12	26
GPER1	31	GRIN2C	26
AGT	30	PTEN	26
CACNB3	30	RGS4	26
HCN4	30	AKR1C3	25
ITPR1	30	GRIK2	25
KIT	30	ITPR3	25

Table 8

List of the Top 40 Genes Based on Frequency in Normal versus MET Set

Gene	Frequency	Gene	Frequency	Gene	Frequency
EGFR	26	ITPR1	20	ANXA6	17
PTK2B	25	RGS2	20	CACNA2D1	17
RYR2	24	SRC	20	FGF13	17
BCL-2	23	ANK2	19	KCNJ5	17
CACNA1D	23	CRABP2	19
SFRP1	23	FYN	19
CXCL12	22	AGT	18
GPER1	22	AKT1	18
KCNJ2	22	CACNA1A	18
KCNQ1	22	CAV1	18
RYR1	22	CX3CL1	18
ABL1	21	EFEMP1	18
ADRA2A	21	FGFR3	18
SLC8A1	21	HBEGF	18
CDK5	20	INHBB	18
DMD	20	KIT	18
EDN1	20	PSEN1	18
FKBP1B	20	ADRB2	17

Based on fold change method, 144 genes upregulated in non-aggressive samples 27 genes down regulated. Among metastatic PTC 138 genes were upregulated and 20 genes down regulated (). The intersection of genes upregulated between both the NAG and MET samples were determined using InteractiVenn17 (). Around 114 genes were seen commonly upregulated in both the sets. The genes unique to NAG set were 30 and for MET was 24 (Figure 3). The list of commonly upregulated genes in both the sets given in Table 9.

Figure 3

Intersection of DEGs among non-aggressive and metastatic set compared to normal samples

Table 9

List of the 114 Genes Commonly Upregulated in Both the Types of PTC

DCSTAMP	LEMD1	LINC02555	MIR31HG	AGR2	ABTB2	CRLF1
KLK10	FAXC	ABCC3	TMEM163	GLT1D1	MIR100HG	CLDN1
GABRB2	FAM230B	KCNJ2	SPTBN2	GALE	TUSC3	KRT19
RXRG	SYT12	KCNN4	SLC34A2	CLDN16	LRRK2	NAT8L
SYTL5	GOLT1A	EGFEM1P	ADTRP	HLA-DQB2	TMEM79	IL17RD
CLDN10	LAMP5	RAB27B	ADORA1	NOD1	NOX4	TNFRSF12A
PRSS2	ZCCHC12	NMU	THRSP	NR2F1-AS1	DOCK9-DT
HMGA2	KLHDC8A	TRDC	ALOX15B	DPP4	B3GNT3
PRR15	CITED1	CD1A	CHI3L1	LPAR5	CORO2A
LRRC52-AS1	NGEF	BRINP1	GLDN	ULBP2	HPCAL4
PDZK1IP1	LRRK2-DT	LIPH	STK32A	MMP16	ECM1
ARHGAP36	GDF15	FAM20A	CTXND1	KISS1R	NRCAM
TMPRSS4	RIMS2	TENM1	ALDH1A3	EVA1A	PLAU
AHNAK2	KCNQ3	KLK11	TIAM1	NFE2L3	TACSTD2
ST6GALNAC5	SCEL	PDZRN4	SYT1	CCL13	PCSK1N
GAP43	LCN2	CDKN2B	COMP	MAMLD1	LINC00891
LAMB3	CDH3	RYR1	SHROOM4	CYP1B1	NHSL2
METTL7B	SLC27A6	LRP4	CEACAM6	IGSF1	INAVA

The fold change in expression for the most frequent genes were retrieved from the microarray data and plotted to compare the differential expression pattern among NAG (n=28), MET (n=26) patients’ samples in comparison to healthy thyroid tissue (n=41). Three genes showed significant differential expression between healthy and thyroid cancer; EGFR (p < 0.05), PTK2B (p < 0.001), KCNN4 (p < 0.001). The 3 genes showed significantly higher expression in NAG and MET samples compared to healthy thyroid tissue (Figure 4).

Figure 4

Box plots for log fold expression from microarray data for the three differentially expressed genes identified from in silico analysis between healthy, non-aggressive and metastatic groups. (A) differential expression of EGFR, (B) differential expression of PTK2B and (C) differential expression of KCNN4. *p < 0.05, ***p < 0.01

In silico Validation of GSEA Using Metascape Analysis

The most frequently present genes across the enriched pathways identified using the absolute GSEA were used to validate the significantly activated cellular pathways between non-aggressive and metastatic samples in comparison to normal samples. The validation was carried out using Metascape relying on large and well annotated cellular pathways derived from gene ontology.18,19 The analysis showed that calcium ion transport, positive regulation of protein phosphorylation and signaling by receptor tyrosine kinase were enriched in non-aggressive PTC (Figure 5A). In the metastatic PTC, significantly activated pathways included positive regulation of protein phosphorylation, MAPK cascade, apoptotic and growth signaling pathways (Figure 5B). Interestingly, although the MAPK pathway activation was present in both NAG and metastatic thyroid the data showed that MAPK pathway came up 3 times in the metastatic set.

Figure 5

Metascape analysis for the high frequent genes from (A) normal versus non-aggressive set and (B) normal versus metastatic set

Similarly, when the commonly upregulated genes were input in Metascape, positive regulation of protein phosphorylation, extracellular matrix organization and cellular response to transforming growth factor beta stimulus pathways were identified (Figure 6).

Figure 6

Metascape for DEGs commonly upregulated in both non-aggressive and metastatic PTC

Analysis of the immune component using the enriched genes from both the NAG and metastatic PTC revealed that NAG has less inflammatory component than the metastatic PTC as shown by imbalance in the M1/M2 ratio as well as the decrease in the NK fraction in the metastatic compared to the non-aggressive PTC. In addition, increase in memory:naïve B-cell ratio was observed in NAG set (Figure 7).

Figure 7

Immune cells enriched in non-aggressive and metastatic PTC in comparison to normal thyroid tissue

DEGs and Enriched Pathways of Thyroid Cancer Across Different Populations

Differentially expressed genes from the microarray data available from other populations such as Ukraine, Brazil and South Korea were analyzed and the top genes upregulated in thyroid cancer in each population were subjected to pathway analysis using Metascape. The results identified unique set of pathways activated for each population. However, key pathways known to be affected in thyroid cancer such as PI3 kinase, MAP kinase and tyrosine metabolism were identified across the various populations (Figures 8–10) indicating that MAPK pathway is probably commonly activated in thyroid cancer across different population cohort. Interestingly, the study identified response to steroid hormone and hormone metabolism activated more in the Ukrainian patients whereas response to inorganic substance and small molecule metabolism was detected in Brazilian thyroid cancer patients. In case of South Korean patients, viral entry and wound healing pathways were observed.

Figure 8

Pathway analysis using Metascape on Ukrainian thyroid cancer samples

Figure 9

Pathway analysis using Metascape on Brazilian thyroid cancer samples

Figure 10

Pathway analysis using Metascape on South Korean thyroid cancer samples

Drugs Associated with Genes Identified from Bioinformatics Analysis

The pharmacoinformatics search using the differentially enriched genes to search the DrugBank database identified many FDA approved drug targets that are currently used in treating metastatic thyroid cancer. Most of those drugs targeted EGFR and vascular EGFR (Table 10). In addition, the search also identified other drugs that target molecules linked to thyroid cancer (Table 11) some of which might be useful as target for pre-clinical trials to treat thyroid cancer pending future studies.

Table 10

List of Drugs Approved by FDA to Treat Thyroid Cancer

Drugs Approved for Thyroid Cancer Treatment	Target Known	Stage of Thyroid Cancer
Cometriq (Cabozantinib-S-Malate)49–51,98	VEGFR	Differentiated and spread; metastasized
Vandetanib99–103	VEGFR and EGFR inhibitor	Metastasized

Table 11

List of Drugs Related to Other Genes Possibly Involved in Thyroid Cancer

Gene Symbol	Drugs Known to Target the Gene	Conditions Associated	Mechanism
KCNQ1	Enflurane	General Anesthesia104,105	Voltage-gated Potassium Channels inhibitor106
	Promethazine	Sedative therapy, Allergic conjunctivitis	Voltage-gated Potassium Channels inducer
CACNA1D107,108	Isradipine35	Hypertension	Calcium channel blocker
PTK2B109	Genistein63,64	Calcium deficiency	Unknown
	Leflunomide66,67	Rheumatoid Arthritis	Regulation of autoimmune lymphocytes
	Fostamatinib71–83	Chronic immune thrombocytopenia	Tyrosine kinase inhibitor
BCL-2	Navitoclax31,110	Solid tumors	Targets BCL-2 family proteins

Drugs targeting the genes specific for population were also searched. Only the population of Ukraine and Brazil showed drugs targeting the genes CRABP1 and MAPK4 respectively (Table 12).

Table 12

List of Drugs Targeting the Genes Highly Upregulated in Population Specific Set

Population	Gene	Drugs Known to Target the Gene	Conditions Associated	Mechanism
Ukraine	CRABP1	Alitretinoin, Tretinoin	Vit A deficiency, eczema	Activates retinoid receptors
Brazil	MAPK4	Fostamatinib	Chronic immune thrombocytopenia	Inhibitor of spleen tyrosine kinase
South Korea	LAMB3	–	–	–

In vivo Validation of DEGs Using NGS

The 3 genes identified from in silico analysis; EGFR, PTK2B and KCNN4 were validated on a cohort of 6 well characterized thyroid carcinoma tissue biopsies collected from patients in the UAE. RNA seq data analyzed for the expression of the genes identified from in silico analysis revealed a significantly higher transcript value for the genes KCNN4 (p < 0.001) and EGFR (p < 0.05) in late PTC samples in comparison to early thyroid cancer samples. PTK2B showed relatively higher expression trend in late PTC samples compared to early (Figure 11).

Figure 11

Differential gene expression in six tissue biopsies from thyroid cancer patients from UAE. *p < 0.05, ***p < 0.01

In vivo Validation of DEGs Using TNM Plot

The expression values for the three genes identified from in silico analysis; EGFR, KCNN4 and PTK2B was examined using an independent larger cohort of RNAseq data from 502 thyroid cancer. The data showed that in this cohort the three genes had higher expression in tumour compared to normal thyroid samples with EGFR (p < 0.01), KCNN4 (p < 0.0001) and PTK2B (p < 0.0001). Thus, the expression fold change confirmed the results from both the microarray and the tissue biopsy for the expression of the EGFR, PTK2B and KCNN4 genes (Figure 12).

Figure 12

TNM Plot output of the three differentially expressed genes identified from in silico analysis on large independent cohort of 58 normal and 502 non-aggressive and 8 metastatic thyroid cancer cases. (A) differential expression of EGFR, (B) differential expression of PTK2B and (C) differential expression of KCNN4

Discussion

This study identified cellular pathways unique to non-aggressive and metastatic PTC as well as common between the two different entities. Interestingly, many of the genes and pathways overlapped between the two clinical groups, these include calcium and potassium ion transport and tyrosine kinase and protein phosphatase pathways. The NAG group showed more unique association with regulation of hormone levels and cell signaling related to hormone whereas the study identified more impact of MAPK activation as well as activation of other cancer hallmark pathways such as regulation of apoptosis and growth in the metastatic pathways.

Validating the data using pathway analysis from the differential expressed genes across different populations showed that MAPK is active in diverse populations including that from Ukraine (Europe), Brazil (South America) and South Korea (Asia). However, each population had set of unique cellular pathways activated. Ukrainian patients had more pathways linked to response to hormones and hormone metabolic processes. Brazilian patients had pathways linked to environmental triggers including response to inorganic substances and vitamin A metabolism suggesting the perhaps pollution20 and poor diet21 may have partially contributed to PTC cases. The South Korean patients seem to have more cancer hallmark pathways activated such as apoptosis and growth in addition to possibly pro-inflammatory pathway activation as shown by the activation of viral entry into the cell pathway suggesting that there might be additional genetic components within that population leading to chronic inflammatory response which when not treated immediately might lead to PTC. Overall, pathway analysis indicated that PTC is highly complex disease with high level intra-tumoral heterogeneity as shown by the activation of different cellular pathways across different populations but with a commonly activated pathway such as the MAPK related pathways. This difference in activated pathways might be reflected by the diverse genetics between the different population. For example, Brazilian population have shown prevalence to inherited TP53 mutation which can lead to tumours in multiple tissue types as characterised by Li-Fraumeni patients.22 Also the Brazilian population show prevalence to mutations in BRCA1 and BRCA2 which are DNA repair genes.23 Ukrainian population show more prevalence to a different DNA repair gene; RAD51.24 However, the South Korean population seem to have prevalence of deletion mutations in Sialic Acid Binding Ig Like Lectin 14 (SIGLEC14)25 which is linked to inflammasome activation in macrophage.26 Thus, the different mutations in the different populations may have role in shaping the transcriptomics profile in those populations. In addition, GSEA of the pathways identified unique upregulated genes for each of the population; CRABP1 for Ukrainian population, MAPK4 for Brazil and LAMB3 for South Korea. DrugBank search using those genes identified retinoid receptors for CRABP1 such as Tretinoin and Alitretinoin which are used for eczema, Fostamatinib for MAPK4 which is used for chronic immune thrombocytopenia. Taken together, the transcriptomic analysis indicated the possibility of repurposing different drugs in different populations for thyroid cancer treatment.

The immune response analysis suggested an imbalance in the tumour immune microenvironment as it showed that NAG has more inflammatory component than metastatic thyroid. This is partly demonstrated by the fact that NAG has both resting and activated NK fraction and higher memory:naïve B-cell ratio whereas metastatic cancer did not show the NK fraction suggesting that the disease stage has passed the inflammatory stage to the cancer stage. This is supported by the fact that other studies showed that in PTC, NK cell infiltration is in early stages of PTC is higher compared to the metastatic stages.27 Additionally, the immune analysis showed an imbalance in the M1/M2 ratio in both the the NAG and MET types of PTC with slightly higher ratio in the metastatic stage indicating that in this cohort might have different mechanism of PTC pathogenesis warranting further studies of the genes involved in the M1 and M2 polarization in thyroid cancer as done in previous studies.28

The study identified the following targets linked to PTC initiation and progression: BCL2, CACNA1D, KCNQ1, KCNN4, EGFR and PTK2B.

B cell Lymphoma-2 (BCL2) is anti-apoptotic protein responsible for inhibiting programmed cell death or apoptosis.29 Aksoy et al found that lower BCL2 expression in thyroid cancer supports the formation of oncocytic neoplasms in early thyroid cancer stages by inhibiting apoptosis of tumor cells.30 This finding from this study supports the results obtained in our study where the frequency of BCL2 overexpression is present in both the NAG and metastatic groups. In addition, few studies have shown that BCL2 is likely to be involved in early PTC as few studies have shown that BCL2 expression decreases in microcarcinomas of PTC30 which indicates that it is probably not a reliable prognostic marker since it is probably involved in very early PTC and continues in the metastatic phase. However, its related drugs such as Navitoclax31 might be useful in treating some forms of PTC.

One of the recurrent activated pathways identified from this study is related to ion transport and more specifically calcium and potassium transport. Many genes related to calcium and potassium transport were identified. CACNA1D gene is responsible for regulating positively charged calcium channels (CaV1.3) across cell membranes and specifically adrenal gland to form alpha-1 subunit. These subunits act as pores to calcium ions to flow through. It is also involved in the regulation of adrenal hormones production such as aldosterone which maintains blood pressure and fluid balance in the body.32,33 Somatic mutations of CACNA1D is associated with tumorigenesis such as in adrenal aldosterone-producing adenomas.33 Interestingly, it has been shown that cancer cells can undergo oncogenic switch by transforming apoptosis inducing Ca influx pathway to proliferative calcium influx which in turn can promote growth and apoptosis resistance in cancerous cells.34 This was also confirmed by the fact that pathway analysis showed the activation of calcium ion transport pathways in both NAG and metastatic PTC. The results from this study, showed that CACNA1D is more frequently overexpressed in the NAG and metastatic PTC compared to healthy suggesting that it is probably involved in PTC progression. In addition, the results indicates that although drugs that targets CACNA1D such as Isradipine35 are used to treat hypertension by regulating the calcium transport, they may help in treating some of the thyroid cancer patients.

Another gene implicated is the KCNQ1 gene which belongs to family genes responsible for potassium channels formation. Channels formed by KCNQ1 genes are located in the inner ear, cardiac muscles, kidney, liver, intestine and stomach. Voltage gated K+ channels (Kv1.3) were identified as novel tumor markers36 Somatic mutations of KCNQ1 and specifically KCNA3 promoter’s methylation contributes to gene silencing37 and dysregulation of potassium ion transport which in turn causes several disease such as cardiovascular diseases, sudden infant death syndrome and cancers.38–41 The results showed that KCNQ1 was the top most frequently present gene across the significantly activated cellular pathways in NAG PTC. However, it remains to be seen whether the drugs that targets KCNQ1 such as Promethazine and Enflurane which are sedative drugs (Table 10) might be worth considered for thyroid cancer pending future studies.

Another gene identified from this study that is implicated in thyroid cancer is Potassium Calcium-Activated Channel Subfamily N Member 4 (KCNN4), a known oncogene, very recently was reported to be upregulated in PTC and was proposed as a diagnostic and prognostic marker for PTC.42 Apart from thyroid cancer, differential expression of KCNN4 in various cancers was indicated either in poor prognosis, drug resistance and/or poor survival.43–45 In the present study, KCNN4 was occurring in both the datasets with high frequency and showed approximately 2 fold change in expression. Hence, potassium calcium activated channels can be targeted to control the progression of PTC to metastatic phase. Interestingly, the results of the RNAseq from the clinical biopsies of both the NGS as well as the TNMplot carried out in this study showed that KCNN4 is significantly overexpressed in metastatic and non-aggressive compared to normal PTC (p < 0.001).

The results also found Epidermal Growth Factor Receptor (EGFR) to be associated with thyroid cancer. EGFR is known to mediate cell proliferation and survival signaling pathways. The transmembrane tyrosine kinase receptor is expressed in different subtypes of cancers such as thyroid carcinoma, glioblastoma and lung cancer.46 EGFR signaling pathways are altered in human cancers due to somatic mutation, gene amplification and protein overexpression which are associated with aggressiveness of the disease and poor survival.47 In this study, EGFR is the most frequently differentially expressed gene in metastatic PTC and also present halfway in the non-aggressive set suggesting that EGFR play a key role in PTC progression and metastasis. Interestingly, search for Thyroid cancer treating drugs identified many FDA approved drugs that targets EGFR including Cabozantinib-S-Malate and Vandetanib.48–60 The results of the RNAseq of the clinical biopsies carried out in this study showed that EGFR is significantly overexpressed in metastatic and non-aggressive compared to normal PTC (p < 0.05).

Another protein identified is Protein tyrosine kinase 2 beta (PTK2B). This has many functions including regulator of cell growth, survival, proliferation and invasion.61 It encodes a cytoplasmic protein tyrosine kinase that is involved in calcium-induced regulation of ion channels and activation of the MAP kinase signaling pathway. Methylated PTK2B favouring overexpression is linked to c-Src activation, development of Pyk2/c-Src complex and the activation of ERK/MAPK signaling pathway. Activation of ERK/MAPK signaling pathway is responsible for regulating the activation of more than 160 downstream signaling transcription factors affecting cancer progression62. The results of the RNAseq from the clinical biopsies of both the NGS as well as the TNMplot carried out in this study showed that PTK2B is significantly overexpressed in metastatic and non-aggressive compared to normal PTC (p < 0.001). Since it is involved in calcium ion regulation and MAPK activation, the drugs which target PTK2B include Genistein63,64 (used to treat calcium deficiency), Leflunomide65–70 (used to treat rheumatoid Arthritis) and Fostamatinib71–83 (used to treat chronic immune thrombocytopenia). Few studies have shown some links between EGFR and PTK2B. Notably, a recent report indicated that overexpression of EGFR and focal adhesion kinases (FAKs) correlated with PTC progression more specifically in aggressive clinicopathological condition and lymph node metastasis84. PTK2B is one of the FAKs, also known to be associated with lymph nodes, tumor size and pathologic state in thyroid cancer samples.85 Moreover, a combinatorial drug Crizotinib (receptor tyrosine kinase targeting drug) was proven effective in reducing tumor size of triple negative breast cancer graft when used along with erlotinib (EGFR targeting drug).86 PTK2B is also one of the targets for the drug Crizotinib.87 The current study also indicates that the metastatic samples were enriched with EGFR and PTK2B genes and the combination might be effective in treating aggressive or metastatic PTC. Therefore, since PTK2B is linked to EGFR, MAP kinase activation and calcium ion transport, it is probably an attractive therapeutic target and since it is linked with poor survival it can be a good prognostic target.

In summary, the absolute gene set enrichment and the pathway analysis indicated strongly that most of pathways overlapped among the non-aggressive and metastatic PTC. The key regulatory proteins among the pathways integrated in PTC pathophysiology are receptor protein tyrosine kinases, calcium channels, potassium channels, potassium activated calcium channels and MAP/ERK kinase family. The genes involved in these processes were seen occurring in high frequency and also seen upregulated in both the datasets which could be used as potential therapeutic targets to treat PTC.

Conclusions

In conclusion, the differentially activated cellular pathways and genes from this study showed the involvement of ion transport as well as other cancer related pathways including tyrosine kinase and protein phosphatase and the modulation of MAPK related pathways in the initiation of PTC during the non-aggressive phase and further progression to PTC metastatic phase. Transcriptomic analysis in different populations highlighted common and unique pathways involved in thyroid cancer pathogenesis further highlighting its heterogeneity. Understanding of genes mediated pathways during carcinogenesis, invasion and metastasis can have significant clinical outcome in developing better prognostic assays and molecular inhibitors that can replace classic generalized PTC treatments.

In addition, the transcriptomics analysis in this study identified interesting putative prognostic targets including EGFR, PTK2B, KCNQ1, KCNN4, BCL2 and CACNA1D which may be involved in key mechanisms of thyroid cancer. EGFR, PTK2B and KCNN4 showed significant higher expression in non-aggressive and metastatic compared to normal using PTC clinical biopsies. Search for corresponding drugs identified FDA approved drugs such as Vandetanib as well as other drugs that may prove useful treating the PTC.