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

Introduction

Thyroid cancer was ranked as the most common endocrine malignancy.¹ 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.⁴ The mortality rate of thyroid cancer is considered to be low, whilst the reoccurrence and persistence of the disease is still considered high.⁵

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.⁶

Table 1 List of Subtypes of Thyroid Carcinoma and the Current Treatment Provided

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/BRAF⁷. However, although such studies suggested that thyroid cancer has high degree of intra-tumoral heterogeneity,⁸ 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 pathogenesis⁹. 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.¹⁰ 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

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.¹⁰ 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 (http://software.broadinstitute.org/gsea/downloads.jsp).¹¹ 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 (https://www.gsea-msigdb.org). 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.¹⁰ 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 (https://metascape.org/)¹³ 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 repository¹⁴ 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.¹⁰ 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

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.¹⁵ 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 DESeq2³⁰ 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 (https://tnmplot.com/analysis/)¹⁶ 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 Supplementary Material.

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

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

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

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

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

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 (Supplementary Material). The intersection of genes upregulated between both the NAG and MET samples were determined using InteractiVenn¹⁷ (http://www.interactivenn.net). 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.

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

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

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

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

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

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 pollution²⁰ and poor diet²¹ 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.²² Also the Brazilian population show prevalence to mutations in BRCA1 and BRCA2 which are DNA repair genes.²³ Ukrainian population show more prevalence to a different DNA repair gene; RAD51.²⁴ However, the South Korean population seem to have prevalence of deletion mutations in Sialic Acid Binding Ig Like Lectin 14 (SIGLEC14)²⁵ which is linked to inflammasome activation in macrophage.²⁶ 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.²⁷ 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.²⁸

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.²⁹ 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.³⁰ 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 PTC³⁰ 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 Navitoclax³¹ 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.³³ 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.³⁴ 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 Isradipine³⁵ 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 markers³⁶ Somatic mutations of KCNQ1 and specifically KCNA3 promoter’s methylation contributes to gene silencing³⁷ 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.⁴² 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.⁴⁶ 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.⁴⁷ 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.⁶¹ 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 progression⁶². 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 Genistein^63,64 (used to treat calcium deficiency), Leflunomide^65–70 (used to treat rheumatoid Arthritis) and Fostamatinib^71–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 metastasis⁸⁴. PTK2B is one of the FAKs, also known to be associated with lymph nodes, tumor size and pathologic state in thyroid cancer samples.⁸⁵ 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).⁸⁶ PTK2B is also one of the targets for the drug Crizotinib.⁸⁷ 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.