PMMB 2022, 5, 1; a0000269. doi: 10.36877/pmmb.a0000269 http://journals.hh-publisher.com/index.php/pmmb Original Research Article Differential gene expression analysis of papillary thyroid carcinoma reveals important genes for lymph node metastasis Wan Fahmi Wan Mohamad Nazarie1, Azliana Mohamad Yusof 2, Francis Yew Fu Tieng1, Nur Fadhlina Mohamad Pakarulrazy1, Rohaizak Muhammad3, Shahrun Niza Abdullah Suhaimi3, Nani Harlina Md Latar3, Sazuita Saidin1, Isa Mohamed Rose4, Learn-Han Lee5, Nurul-Syakima Ab Mutalib1, 5, 6* Article History 1UKM Medical Molecular Biology Institute, Universiti Kebangsaan Malaysia, Bandar Tun Razak, 56000 Cheras, Kuala Lumpur, Malaysia; wanfahmi5785@gmail.com (WFWMN); francistieng@yahoo.com.my (FYFT); nurfadhlinamohamadpakarulrazy@gmail.com (NFMP); sazuita@ukm.edu.my (SS) 2Cytogenetics and Molecular Diagnostics Laboratory, Pantai Premier Pathology Sdn. Bhd., 55100 Kuala Lumpur, Kuala Lumpur, Malaysia; azlianayusof@gmail.com (AMY) 3Department of Surgery, Faculty of Medicine, Universiti Kebangsaan Malaysia, Bandar Tun Razak, 56000 Cheras, Kuala Lumpur, Malaysia; rohaizak@ppukm.ukm.edu.my (RM); shahrunniza@ppukm.ukm.edu.my (SNAS), naniharlinalatar@ukm.edu.my (NHML) 4Department of Pathology, Faculty of Medicine, Universiti Kebangsaan Malaysia, Bandar Tun Razak, 56000 Cheras, Kuala Lumpur, Malaysia; isa@ppukm.ukm.edu.my (IMR) 5Novel Bacteria and Drug Discovery Research Group, Microbiome and Bioresource Research Strength, Jeffrey Cheah School of Medicine and Health Sciences, Monash University Malaysia, 47500 Subang Jaya, Selangor, Malaysia; lee.learn.han@monash.edu (L-HL) 6Faculty of Health Sciences, Universiti Kebangsaan Malaysia, 50300 Kuala Lumpur, Kuala Lumpur, Malaysia *Corresponding author: Nurul Syakima Ab Mutalib, UKM Medical Molecular Biology Institute, Universiti Kebangsaan Malaysia, Bandar Tun Razak, 56000 Cheras, Kuala Lumpur, Malaysia; syakima@ppukm.ukm.edu.my (N-SAM) Received: 05 June 2022; Received in Revised Form: 11 July 2022; Accepted: 13 July 2022; Available Online: 17 July 2022 Abstract: Background: Papillary thyroid carcinoma (PTC) is the most prevalent thyroid cancer. We explored the differential gene expression (DEG) of the transcriptional regulation in thyroid cancer with the main aim of determining the genes involved in lymph node metastasis (LNM) in PTC. Methods: We employed a bioinformatics pipeline for RNA-Seq analysis in PTC with and without LNM at the gene expression level. We performed read mapping, read quantitation, and DEG analysis using STAR, Cufflinks, and Cuffdiff, respectively. Subsequently, functional annotation and pathway enrichment were carried out using FunRich (functional enrichment analysis tool). Results: Expression profiling revealed changes in the PTC with LNM (33 genes at p-value < 0.05 and log2 fold change |1.0|) mailto:wanfahmi5785@gmail.com mailto:francistieng@yahoo.com.my mailto:sazuita@ukm.edu.my mailto:azlianayusof@gmail.com mailto:rohaizak@ppukm.ukm.edu.my mailto:shahrunniza@ppukm.ukm.edu.my mailto:isa@ppukm.ukm.edu.my PMMB 2022, 5, 1; a0000269 2 of 19 compared to the adjacent normal thyroid, whereas 69 genes showed differential expressions in the PTC with lymph node negative (LNN) versus adjacent normal thyroid. We identified 31 significant DEGs in PTC LNM versus PTC LNN; and 44 significant DEGs between adjacent normal thyroid tissues from PTC LNM and PTC LNN. The biological processes of genes expressed at higher levels in PTC LNM compared to PTC LNN were involved in cell communication, energy pathway, and metabolism, whereas ion transport, energy pathways, and regulation of nucleobase, nucleoside, nucleotide, and nucleic acid metabolisms, were downregulated. Conclusion: These findings provide further evidence for the role of cellular transport regulatory processes in metastatic disease. Keywords: Papillary thyroid carcinoma; lymph node metastasis; lymph node negative; differentially expressed genes 1. Introduction The number of new cases of thyroid cancer including papillary thyroid carcinoma (PTC) is steadily increasing over the past decade especially in Asian countries compared to European countries[1]. Thyroid cancer has a high survival rate because of early detection through screening and advances in medical treatment[2,3]. There are four types of thyroid cancer namely follicular (FTC), medullary (MTC), anaplastic (ATC) and papillary (PTC), and the latter is the most common, contributing 80% to 85% of all thyroid cancers[4,5]. The survival rate for PTC is high (98%) in Europe and North America due to improved diagnosis, management, and treatment of the disease[6]. Yet, a subset of PTC patients displayed aggressive features associated with poor prognosis, and one of the features is the presence of lymph node metastasis (LNM)[7]. Several studies have shown the different molecular landscapes of LNM in various types of cancer, including oral, pancreatic, colorectal, and breast cancer[8]. However, the gene expression profile in LNM of PTC has not been extensively done, with only a handful of published findings[9]. Understanding the molecular and cellular mechanisms of LNM in PTC is required to improve prognostication. In this study, we investigated the gene expression profiles of PTC with and without LNM and conducted a comprehensive analysis of gene expression of thyroid tissue samples via RNA-Seq. The differentially expressed genes (DEGs) identified in PTC with and without LNM might be the potential biomarkers and provide the basis for further exploration of the mechanisms of metastatic disease in PTC. 2. Materials and methods 2.1. Ethics statement and sample collection This study was conducted according to the Universiti Kebangsaan Malaysia Research Ethics Committee (UKMREC) (reference: UKM 1.5.3.5/244/UMBI-2015-002). Ten fresh- PMMB 2022, 5, 1; a0000269 3 of 19 frozen tumour PTC tissues specimens from the cancer Biobank were subjected to cryosectioning and stained with haematoxylin and eosin (H&E). The pathologist reviewed the percentage of tumour cells and normal cells in the slides. Tumour tissues containing more than 80% of cancer cells and normal tissues with less than 20% necrosis were subjected to RNA extraction. 2.2. Total RNA isolation According to the manufacturer’s protocol, total RNA was isolated from the samples using AllPrep DNA/RNA/miRNA Isolation Kit (Qiagen, Germany). The total RNA quality and quantity were carried out using an ND-1000 spectrophotometer (Thermo Scientific, USA) and Qubit™ fluorometer (Invitrogen, USA). Before library preparation, the RNA quality was measured using an Eukaryote Total RNA Nano chip on Bioanalyzer 2100 (Agilent Technologies, USA). 2.3. Library preparation and RNA sequencing Libraries were constructed using the Ion Total RNA-Seq v2 kit (Life Technologies, USA) according to the manufacturer’s protocol. Briefly, 2 µg of total RNA was subjected to ribosomal removal using RiboMinus Eukaryote kit (Life Technologies, USA). ERCC RNA Spike-In Control Mixes (Thermo Scientific, USA) was added into rRNA depleted-RNA before RNA fragmentation using RNase III. The fragmented RNA was quantified using Qubit™ fluorometer (Invitrogen, USA) and Agilent RNA 6000 Nano Kit (Agilent Technologies, USA), followed by hybridization and ligation to adaptors, reverse-transcribed, purified, size-selected, and amplified. The final library concentration was assessed using Agilent High Sensitivity DNA chip (Agilent Technologies, USA) and normalized from 100 pM to 130 pM for clonal amplification. Clonal amplification was performed on Ion Chef System using Ion PI IC 200 kit followed by sequencing using PI BC v2 chip on Ion Proton system (all from Life Technologies, USA). 2.4. Read mapping and gene quantification Raw reads (FASTQ files) were processed using AfterQC[10] to clean the data and obtain high-quality reads by removing reads that contain adapter or poly-N and low quality reads. GC content of the clean data was then calculated. All downstream analyses were based on clean data with high quality. Sequencing reads were mapped with STAR (v2.0.9[11,12] to the human genome sequence assembly (GRCh38.89). The mapped reads of each sample were assembled using Cufflinks[13]. 2.5. Differential gene expression Differential gene expression (DEG) analysis was performed using Cuffdiff[13]. Four comparisons were performed: PTC LNM (lymph node metastasis) versus PTC LNN (lymph node negative), PTC LNM versus normal thyroid, PTC LNN versus normal thyroid, and normal thyroid from PTC LNM patients versus normal thyroid from PTC LNN. Genes with PMMB 2022, 5, 1; a0000269 4 of 19 a p-value of less than 0.05 and showing a log2 fold change |1| were considered as significantly differentially expressed. 2.6. Gene Ontology and KEGG pathway enrichment analysis Next, we performed a Gene Ontology (GO) analysis of these DEGs. Functional enrichment analysis focusing on biological processes was performed on datasets by using FunRich software[14]. Settings were set at a default parameter, and fold enrichment with unadjusted p-values were computed using a modification of the Fisher's exact test. We used KEGG, a database resource for understanding high-level functions and utilities of the biological system, to test the statistical enrichment of differential expression genes in KEGG pathways[15]. 3. Results 3.1. Data characterization of sequencing and mapping Our RNA-Seq yielded 22.2 to 80.3 million reads. Using STAR aligner, the proportion of reads that were successfully mapped to the Ensembl reference genes ranged from 77% to 85%, indicating that the data produced are of high quality. Gene expression levels were estimated using Cufflinks plotted in Figure 1 with Pearson correlation coefficients of the estimated abundances between the two tissues. It showed that the gene quantification with Cufflinks was consistent across different comparisons. Figure 1. Differential gene expression analysis of PTC and normal thyroid tissue. (A) Density plots showing the expression level distribution for all genes in the two tissues. FPKM = fragments per kb of transcript per million fragments mapped. (1B) The scatter plot of global expression for each comparison between samples where the Pearson correlation coefficient (PCC) is shown. Each dot represents one gene that has detectable expression in either tissue. PMMB 2022, 5, 1; a0000269 5 of 19 3.2. Differential gene expression We identified 33, 69, 31, and 44 significant DEGs between the PTC LNM versus normal, PTC LNN versus normal, PTC LNM versus PTC LNN, and normal thyroid from PTC LNM versus normal thyroid from PTC LNN, respectively (Table 1). The complete lists of the DEGs were provided in Supplementary Table S1. Up- and down-regulated genes comparing PTC LNM versus PTC LNN were shown in Figure 2. In the group comparison of PTC LNM versus normal thyroid, up-regulated transport-related genes in the PTC LNM included GABRE and NUP160, while down-regulated genes were RYR3 and HBA2. In the PTC LNN versus normal thyroid, the upregulated genes included LCN2, RYR3, KCNQ3, GRIA3, SLC13A5, and SLC34A2 but downregulated transport-related genes were not found. While in the group comparison of PTC LNM versus PTC LNN, there were two up-regulated transport-related genes, SLC13A5 and SYBU, and three down-regulated transport-related genes, which were GRID1, SCN8A, and TRPM2. In the normal thyroid from both PTC LNM and PTC LNN, 17 up-regulated and 27 down-regulated with only two transport-related genes were being differentially expressed. The down-regulated transport-related genes were RYR3 and TMEM199. The overlap of the DEGs among the four different groups was shown in a Venn diagram in Figure 3. Table 1. Significant differently expressed genes (p-value < 0.05 and log2 fold change |1|) in PTC LNM versus PTC LNN, normal thyroid (PTC LNM) versus normal thyroid (PTC LNN), PTC LNM versus normal, and PTC LNN versus normal group. Group PTC LNM versus normal thyroid PTC LNN versus normal thyroid PTC LNM versus PTC LNN Normal thyroid (PTC LNM) versus normal thyroid (PTC LNN) Up-regulated gene 14 50 14 17 Down-regulated gene 19 19 17 27 TOTAL 33 69 31 44 LNM: lymph node metastasis; LNN: lymph node negative; PTC: papillary thyroid cancer PMMB 2022, 5, 1; a0000269 6 of 19 Figure 2. Expression level of five up- and down-regulated in PTC LNM versus PTC LNN. (2A) Boxplots illustrated FPKM expression value in PTC LNM versus PTC LNN (SLC13A5 and SYBU) and (2B) boxplots demonstrated FPKM expression value in PTC LNM versus PTC LNN (GRID1, SCN8A, and TRPM2). Figure 3. Venn diagram representing differentially expressed transport-related genes overlapping between the samples. The Gene list from each group comparison showed the frequency of transport-related genes. PMMB 2022, 5, 1; a0000269 7 of 19 3.3. Enriched functions of differentially expressed gene The functional analysis showed that the primary biological processes altered in PTC LNM versus normal thyroid were protein metabolism, metabolism, energy pathways, cell growth and/or maintenance, and transport. These processes are frequently involved in cancer. To improve our understanding of the function of the DEGs, we performed an enrichment analysis of Gene Ontology (GO) for the up-and down-regulated genes, as illustrated in Figure 4. First, we performed enrichment tests to detect the functional categories for significantly up- and down-regulated genes detected by Cuffdiff in the PTC LNM, PTC LNN, and adjacent normal tissue using the FunRich software. The GO categories that were greatly enriched in the regulated genes from PTC LNM versus normal thyroid, PTC LNN versus normal thyroid and PTC LNM versus PTC LNN were selected. The up- and down-regulated genes in all comparisons in PTC were categorized into 177 functional categories. We identified 14 up- and 19 down-regulated genes in PTC LNM versus normal thyroid and 50 up-regulated and 19 down-regulated genes in PTC LNN versus normal thyroid. For instance, the biological processes for up-regulated genes, which include “protein targeting”, “inflammatory response”, “regulation of cell growth”, “cell differentiation” and “transport” are essential to cancer progression[16,17]. Whilst the upregulated genes in PTC LNN versus normal included “cell growth and maintenance”, “energy pathway”, “pyrimidine salvage”, “transport”, “ion transport” and “protein metabolism”. On the contrary, genes that were downregulated in PTC LNN versus normal thyroid included those involved in “cell communication”, “cell differentiation”, “cell proliferation”, “steroid hormone receptor signaling pathway”, and “vesicle-mediated transport” and “apoptosis”. Interestingly, the biological process of “transport” appeared in most of the group comparisons suggesting that transport genes could potentially govern the metastasis cancer progression in thyroid cancer. Detailed analysis of functional annotation was performed using FunRich to analyze which KEGG pathways were enriched with PTC-specific down-regulated genes. The pathways enriched with DEGs are listed in Table 2. RHO GTPases activate CIT interaction pathway was affected in two out of four pairwise comparisons (PTC LNM versus normal thyroid and PTC LNM versus PTC LNN, and other gene regulation alterations in the extracellular matrix (ECM) organization pathway were involved in cancer tissue. Furthermore, the ChREBP activates metabolic gene expression pathway was enriched in the DEGs identified from the PTC tissue. PMMB 2022, 5, 1; a0000269 8 of 19 Figure 4. Gene enrichment biological processes across the different groups. PMMB 2022, 5, 1; a0000269 9 of 19 Table 2. KEGG pathway of enriched differentially expressed genes. Comparison Pathway ID Pathway Name p-value FDR PTC LNM versus normal thyroid 1270056 Acyl chain remodeling of DAG and TAG 1.272E- 4 2.480E- 2 1269514 RHO GTPases activate CIT 4.759E- 4 4.640E- 2 PTC LNN versus normal thyroid 1270114 ChREBP activates metabolic gene expression 6.05E- 04 2.62E- 01 1270254 Non-integrin membrane-ECM interactions 9.19E- 04 3.99E- 01 1470923 Interleukin-4 and 13 signaling 1.29E- 03 5.59E- 01 1270244 Extracellular matrix organization 1.44E- 03 6.26E- 01 PTC LNM versus PTC LNN 1269514 RHO GTPases activate CIT 5.60E- 04 5.93E- 02 1457800 MET activates PTK2 signaling 7.79E- 04 8.25E- 02 1108786 mucin core 1 and core 2 O- glycosylation 1.22E- 03 1.30E- 01 1270253 Laminin interactions 1.22E- 03 1.30E- 01 1383022 FGFR2 alternative splicing 1.53E- 03 1.63E- 01 PTC LNM (normal) versus PTC LNN (normal) 852705 MicroRNAs in cancer 2.38E- 05 7.58E- 03 142173 Serine biosynthesis (phosphorylated route) 2.98E- 03 9.48E- 01 CIT: citron rho-interacting serine/threonine kinase; ChREBP: carbohydrate response element binding protein; DAG: diacylglycerol; ECM: extracellular matrix; LNM: lymph node metastasis; LNN: lymph node negative; MET: MNNG HOS transforming gene; PTC: papillary thyroid cancer; PTK2: protein tyrosine kinase; TAG: triacylglycerol PMMB 2022, 5, 1; a0000269 10 of 19 GO analysis was further refined to identify which DEGs were involved in transport, which was an important component of metastasis in PTC and illustrated in Table 3. For in silico validation, the expression signatures of selected thyroid-specific genes (GRID1, KCNQ3, and SLC34A) were retrieved from the RNA-Seq data of primary tumors in The Cancer Genome Atlas (TCGA) portal (Figure 5). Table 3. Transport-related genes expression in RNA-Seq data of thyroid cancer. Group PTC LNM versus normal thyroid PTC LNN versus normal thyroid PTC LNM versus PTC LNN Normal thyroid (PTC LNM) versus normal thyroid (PTC LNN) Up- regulated gene GABRE, NUP160 LCN2, RYR3, SLC34A2, KCNQ3, SLC13A5, GRIA3 SLC13A5, SYBU - Down- regulated gene RYR3, HBA2 - GRID1, SCN8A, TRPM2 RYR3, TMEM199 GABRE: gamma-aminobutyric acid type A receptor subunit Epsilon; GRIA3: glutamate ionotropic receptor AMPA type subunit 3; GRID1: glutamate receptor delta-1 subunit; HBA2: haemoglobin alpha 2; KCNQ3: potassium voltage-gated channel subfamily Q member 3; LCN2: lipocalin 2; LNM: lymph node metastasis; LNN: lymph node negative; NUP160: nucleoporin 160; PTC: papillary thyroid cancer; RYR3: ryanodine receptor 3; SCN8A: sodium voltage-gated channel alpha subunit 8; SLC13A5: solute carrier family 13 member 5; SLC34A2: solute carrier family 34 member 2; SYBU: syntabulin; TMEM199: transmembrane protein 199; TRPM2: transient receptor potential cation channel, subfamily M, member 2 PMMB 2022, 5, 1; a0000269 11 of 19 Figure 5. Gene expression overview of TCGA data; higher in thyroid cancer. (5A) GRID1 (5B) SLC34A2 and (5C) KCNQ3. PMMB 2022, 5, 1; a0000269 12 of 19 4. Discussion Our study aimed to identify genes responsible for lymph node metastasis in PTC. We profiled the whole transcriptome of PTC LNM, PTC LNN, and normal thyroid tissues using ribosomal-reduced RNA-Seq. Research on thyroid cancer progression has been well-studied using microarray and RNA-Seq. We identified several transporter-related genes associated with PTC LNM, PTC LNN, and adjacent normal thyroid. Ion transport process happens during the multi-stage cancer progression, from a normal cell to cancer cell and finally to metastasis. Our results suggested an important role of transporter genes during lymph node metastasis in thyroid cancer. We identified 177 DEGs within the various comparison groups. However, a KEGG pathway enrichment analysis showed that the regulated genes were significantly enriched in the RHO GTPases activated CIT pathway in two group comparisons: PTC LNM versus PTC LNN and PTC LNM versus normal thyroid. This pathway is closely related to cancer progression[18,19]. For instance, aberrant activity of Rho small G-proteins, particularly Rac1 and their regulators, is a hallmark of cancer and contributes to the tumorigenic and metastatic phenotypes of cancer cells[20]. Two up-regulated transporter genes (GABRE and NUP160) and two down-regulated genes (RYR3 and HBA2) were identified in PTC LNM versus normal thyroid. The GABRE gene encodes for the gamma-aminobutyric acid type A receptor epsilon subunit, which is a multisubunit chloride channel that facilitates the synaptic transmission in the central nervous system (CNS). The miR‐224/452 cluster is co‐expressed with its host gene GABRE located in the intron of GABRE gene coding for gamma-aminobutyric acid (GABA) receptor[21], and its expression is directly activated by transcription factor E2F1 through transactivation of the GABRE gene[22]. This gene is up-regulated in medulloblastomas, lung cancer, bladder cancer, prostate cancer, and colorectal cancer[23–26]. Many studies reveal that miRNAs are involved in carcinogenesis, and miRNA expression has been found in many types of cancer[27–29]. In contrast, the GABRE∼miR-452∼miR-224 locus has been reported to be down- regulated and hypermethylated in prostate cancer, suggesting its potential as a new epigenetic candidate biomarker for prostate cancer diagnosis and prognosis[24,30]. It has also been up- regulated in hepatocellular carcinoma[31]. MiR-224 is also up-regulated in PTC, MTC, FTC and ATC[32]. Many studies on different types of cancer suggest that miR-224 promotes tumorigenesis in colorectal cancer[33,34], non-small cell lung cancer[35], medullary thyroid carcinoma[36] and hepatocellular carcinoma[37]. In contrast, some studies show that it suppresses tumorigenesis in prostate cancer[24,30], and breast cancer[38,39]. GABRB2, which is part of the GABA receptors gene family, plays a crucial role in the lymph node metastasis of PTC[40]. Its function is closely related to the nervous system, but its role in cancer remains uncertain[41]. We identified two down-regulated genes in PTC LNM versus normal thyroid, namely, NUP160 and RYR3. The NUP160 gene encodes the nuclear pore complex protein Nup160, which mediates molecule transport across the nuclear envelope[42] and is thought to be PMMB 2022, 5, 1; a0000269 13 of 19 involved in cancer progression[43]. Up-regulation of this gene may indicate its role in lymph node metastasis. One study reported a novel frequent fusion gene, NUP160-SLC43A3, in angiosarcoma patients and deduced that this fusion protein might cause rapid tumour progression[44]. This novel fusion could also provide a possible therapeutic target for many human malignancies[45,46]. Our findings suggest an oncogenic role of GABRE and NUP160 in thyroid carcinogenesis. The RYR3 gene encodes a ryanodine receptor, which releases calcium from intracellular storage for use in many cellular processes[47,48]. It is imperative for the growth, morphology, and migration of breast cancer cells, and studies have shown that it is commonly expressed in breast cancer[49,50]. Functional ryanodine receptor is also involved in the apoptosis of in vitro human prostate cancer LNCaP cells[51]. This gene's Gene Ontology (GO) analysis included calcium ion binding and calmodulin binding. Additionally, this gene is also involved in transport, suggesting that it might be involved in the process of metastasis in PTC. In comparing tumour samples of PTC LNM and PTC LNN, we identified two up- regulated transport-related genes (SLC13A5 and SYBU) and three down-regulated transport- related genes (GRID1, SCN8A, and TRPM2). The SLC13A5 (solute carrier family 13 members 5) gene encodes for proteins that transport citrate from the cytoplasm into cells and is mostly expressed in liver cancer[52] via regulation of metabolic processes[53]. It is a sodium- coupled transporter that facilitates the cellular uptake of citrate and plays vital role in synthesizing fatty acids and cholesterol[52]. Syntabulin-kinesin-1 family member 5B (SYBU), a subset of the kinesin motor-adaptor complex is important for axonal transport and is involved in neuronal development [54]. It encodes syntabulin, a microtubule-related protein that facilitates the transport of vesicles to neuronal processes[55]. Functional analysis of SYBU revealed its involvement in various binding function such as protein binding, microtubule binding, syntaxin-1 binding, and kinesin binding[56]. It has been reported to be over-expressed in pancreatic adenocarcinoma, breast cancer and bladder cancer[57]. Thus far, no study has been conducted on this gene involvement in thyroid cancer. In our study, three glutamate receptor subunits: GRID1 (ionotropic glutamate receptor), SCN8A (sodium voltage-gated channel) and TRPM2 (transient receptor potential cation channel) were significantly up-regulated in PTC LNM tissues. Ion transport is believed to be involved in cancer progression, starting with transforming a normal cell into a cancer cell and until the metastasis stage[58]. Moreover, this derangement in the ion transport mechanism is one of the hallmarks of cancer and could be a key event in tumour metastasis [59]. The study by Stepulak et al. also suggested glutamate as a potential growth factor in tumorigenesis [60]. Glutamate receptor subunits: NR1-NR3B, GluR1-GluR7, KA1, KA2, and mGluR1-mGluR8, have been shown to be differentially expressed in multiple cancer cell lines (TE671, SK-NA- S, FTC 238, SK-LU-1, MOGGCCM, RPMI 8226, A549, HT 29, Jurkat E6.1, T47D, LS180, U87-MG, and U343)[61,62] and tumours including thyroid cancer, hepatocellular carcinoma PMMB 2022, 5, 1; a0000269 14 of 19 (GRIK1), colorectal cancer (mGluR4)[63], melanoma (mGluR1)[64], and serous ovarian adenocarcinoma (GRIA2)[65]. The GRID1 gene showed increased expression in colorectal cancer and breast cancer[66], and down-regulated expression in non-small cell lung cancer of node-positive cases[67]. Thereby, GRID1 might play an essential role in the progression of PTC. It could also be used as a prognostic biomarker and a therapeutic target[68]. Another gene up-regulated in PTC LNM is SCN8A which encodes a protein of the ion pore region of the voltage-gated sodium channel (VGSC) and is crucial for the rapid membrane depolarization[69]. This gene is over-expressed in cervical cancer and prostate cancer[70]. Additionally, the expression of VGSC is associated with cancer cell migration, invasion, and metastasis[71]. Inhibition of VGSC subunits has been suggested as a suitable treatment strategy to decrease the metastatic spread of prostate cancer[70]. A detailed review of VGSC by Wang et al. indicated that upregulation of SCN8A gene might increase cancer cell growth and promote metastasis. They also suggested that inhibiting sodium channels in both immune and cancer cells might reduce metastases[37]. Based on our results, the TRPM2 gene, which encodes for the non-selective calcium- permeable cation channel, was up-regulated in PTC LNM. Our study showed concordance with other studies where TRPM2 was over-expressed in many cancers, including but not limited to bladder cancer[72], gastric cancer[73(p2)], lung cancer, and breast cancer[74]. In 2013, Liu and the coauthors reported activation of TRPM2 gene by irradiation via PARP1 activation. They contributed to the irreversible loss of salivary gland function, suggesting that the disordered expression of TRPM2 might be associated with the occurrence of head squamous cell carcinoma (HSCC)[75(p2)]. In addition, a group of researchers from China proposed that TRPM2 was essential in the regulation of migration and survival of SCC cancer cells. At the same time, another study described that TRP channels were clinically valuable for cancer diagnosis and prognosis[76]. Our study identified a prominent association between GRID1, SCN8A, and TRPM2 with LNM in PTC. These findings suggest that the upregulation of GRID1 and SCN8A is related to PTC metastasis. Nevertheless, further studies are required to explore the potential of GRID1 and SCN8A as molecular biomarkers of PTC. Three thyroid-specific gene expression signatures (GRID1, KCNQ3, and SLC34A) were up-regulated in thyroid cancer based on the RNA-Seq data of primary tumours in The Cancer Genome Atlas (TCGA) portal, and this observation was concordant with our results. In our study, the KEGG pathway and GO enrichment analysis showed many transport pathways involved in lymph node metastasis. Future studies could focus on these transporter genes and its related pathways in cancer, such as the glucose and hexose pathways, which had been discussed previously by Adekola et al.[77]. The functional analyses of genes differentially enriched in PTC LNM and PTC LNN further supported the hypothesis that lymph node metastasis was associated with specific essential genes and pathways. The results from our study showed the similarities found in terms of the role of crucial transport pathways, especially from the perspective of metastasis in PTC. We screened the DEGs linked to the transport-related pathways to PTC by using the PMMB 2022, 5, 1; a0000269 15 of 19 TCGA database to confirm our findings. Among the 17 transport-related DEGs, only four genes (GABRE, HBA2, KCNQ3, and SLC34A2) showed similar expressions to our conclusions. Another 13 genes required further validation to confirm the results. 5. Conclusion In conclusion, we have identified key mRNA signatures of PTC LNM and PTC LNN. The DEGs might be promising biomarkers of PTC LNM and could provide a basis for further exploration of the tumorigenesis PTC. These genes may play critical roles in metastasis and survival in PTC patients. These results offer new insights on metastasis in PTC and eventually lead to better diagnosis and possibly new therapeutic targets for thyroid cancer. Author Contributions: N-SAM conceived the idea, supervised the project, and applied for the grant that funded the work. AMY and SS performed laboratory experiments. WFWMN performed the bioinformatics analyses and generated the tables and figures. RM, SNAS and NHML are the thyroid surgeons involved in specimen collection, and IMR is the pathologist who confirms the diagnosis. WFWMN, N-SAM, FYFT, NFMP and L-HL wrote and edited the manuscript. Funding: This study was funded by the Fundamental Research Grant Scheme (FRGS) from the Ministry of Higher Education Malaysia (FRGS/1/2014/SKK01/UKM/03/1). Conflicts of Interest: The authors declare no conflict of interest. References 1. Vaccarella S, Franceschi S, Bray F, et al. Worldwide thyroid-cancer epidemic? The increasing impact of overdiagnosis. N Engl J Med 2016; 375(7): 614-617. 2. Sung H, Ferlay J, Siegel RL, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2021; 71(3): 209-249. 3. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin 2022; 72(1): 7-33. 4. Schlumberger MJ. Papillary and follicular thyroid carcinoma. N Engl J Med 1998; 338(5): 297-306. 5. Limaiem F, Rehman A, Mazzoni T. Papillary thyroid carcinoma. In: StatPearls. StatPearls Publishing; 2022. Accessed April 15, 2022. http://www.ncbi.nlm.nih.gov/books/NBK536943/ 6. La Vecchia C, Malvezzi M, Bosetti C, et al. Thyroid cancer mortality and incidence: a global overview. Int J Cancer 2015; 136(9): 2187-2195. 7. Tan J, Tan J, Qian X, et al. Integrated bioinformatics analysis reveals that the expression of cathepsin S is associated with lymph node metastasis and poor prognosis in papillary thyroid cancer. Oncol Rep 2018; 40(1): 111-122. 8. Hoie J, Stenwig AE, Kullmann G, Lindegaard M. Distant metastases in papillary thyroid cancer. A review of 91 patients. Cancer 1988; 61(1): 1-6. 9. Ab Mutalib NS; Othman SN, Yusof AM, et al. Integrated microRNA, gene expression and transcription factors signature in papillary thyroid cancer with lymph node metastasis. PeerJ 2016; 4: e2119. 10. Chen S, Huang T, Zhou Y, et al. AfterQC: automatic filtering, trimming, error removing and quality control for fastq data. BMC Bioinformatics 2017; 18(Suppl 3): 80. PMMB 2022, 5, 1; a0000269 16 of 19 11. Dobin A, Gingeras TR. Mapping rna-seq reads with star. Curr Protoc Bioinforma Ed Board Andreas Baxevanis Al 2015; 51: 11.14.1-11.14.19. 12. Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinforma Oxf Engl 2013; 29(1): 15-21. 13. Trapnell C, Hendrickson DG, Sauvageau M, et al. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol 2013; 31(1): 46-53. 14. Pathan M, Keerthikumar S, Ang CS, et al. FunRich: An open access standalone functional enrichment and interaction network analysis tool. Proteomics 2015; 15(15): 2597-2601. 15. Kanehisa M, Furumichi M, Tanabe M, et al. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res 2017; 45(D1): D353-D361. 16. Prete A, Borges de Souza P, Censi S, et al. Update on fundamental mechanisms of thyroid cancer. Front Endocrinol 2020; 11: Article 102. 17. Xing M. Molecular pathogenesis and mechanisms of thyroid cancer. Nat Rev Cancer 2013; 13(3): 184- 199. 18. Bustelo XR. RHO GTPases in cancer: known facts, open questions, and therapeutic challenges. Biochem Soc Trans 2018: BST20170531. 19. Haga RB, Ridley AJ. Rho GTPases: Regulation and roles in cancer cell biology. Small GTPases 2016; 7(4): 207-221. 20. Kazanietz MG, Caloca MJ. The rac gtpase in cancer: from old concepts to new paradigms. Cancer Res 2017; 77(20): 5445-5451. 21. Gokhale A, Kunder R, Goel A, et al. Distinctive microRNA signature of medulloblastomas associated with the WNT signaling pathway. J Cancer Res Ther 2010; 6(4): 521-529. 22. Knoll S, Fürst K, Kowtharapu B, et al. E2F1 induces miR-224/452 expression to drive EMT through TXNIP downregulation. EMBO Rep 2014; 15(12): 1315-1329. 23. Ke TW, Hsu HL, Wu YH, et al. Microrna-224 suppresses colorectal cancer cell migration by targeting cdc42. Dis Markers 2014; 2014: e617150. 24. Kristensen H, Haldrup C, Strand S, et al. Hypermethylation of the GABRE\textasciitildemiR- 452\textasciitildemiR-224 promoter in prostate cancer predicts biochemical recurrence after radical prostatectomy. Clin Cancer Res Off J Am Assoc Cancer Res 2014; 20(8): 2169-2181. 25. Boguslawska J, Piekielko-Witkowska A, Wojcicka A, et al. Regulatory feedback loop between T3 and microRNAs in renal cancer. Mol Cell Endocrinol 2014; 384(1): 61-70. 26. Yang J, Bai Y, Cui Y, et al. Expression and biological significance of GABRE in colon cancer: An analysis based on data mining of Oncomine and TCGA databases. Chin J Cancer Biotherapy 2020; 6: 1399-1405. 27. Zhao L, Chen X, Cao Y. New role of microRNA: carcinogenesis and clinical application in cancer. Acta Biochim Biophys Sin 2011; 43(11): 831-839. PMMB 2022, 5, 1; a0000269 17 of 19 28. Osada H, Takahashi T. MicroRNAs in biological processes and carcinogenesis. Carcinogenesis 2007; 28(1): 2-12. 29. Ahmed FE. Role of miRNA in carcinogenesis and biomarker selection: a methodological view. Expert Rev Mol Diagn 2007; 7(5): 569-603. 30. Esfahani M, Ataei N, Panjehpour M. Biomarkers for evaluation of prostate cancer prognosis. Asian Pac J Cancer Prev APJCP 2015; 16(7): 2601-2611. 31. Ladeiro Y, Couchy G, Balabaud C, et al. MicroRNA profiling in hepatocellular tumors is associated with clinical features and oncogene/tumor suppressor gene mutations. Hepatol Baltim Md 2008; 47(6): 1955- 1963. 32. Nikiforova MN, Tseng GC, Steward D, et al. MicroRNA expression profiling of thyroid tumors: biological significance and diagnostic utility. J Clin Endocrinol Metab 2008; 93(5): 1600-1608. 33. Ling H, Pickard K, Ivan C, et al. The clinical and biological significance of MIR-224 expression in colorectal cancer metastasis. Gut 2016; 65(6): 977-989. 34. Fassan M, Cui R, Gasparini P, et al. Mir-224 is significantly upregulated and targets caspase-3 and caspase-7 during colorectal carcinogenesis. Transl Oncol 2018; 12(2): 282-291. 35. Cui R, Meng W, Sun HL, et al. MicroRNA-224 promotes tumor progression in nonsmall cell lung cancer. Proc Natl Acad Sci 2015; 112(31): E4288-E4297. 36. Cavedon E, Barollo S, Bertazza L, et al. Prognostic impact of miR-224 and RAS mutations in medullary thyroid carcinoma. Int J Endocrinol 2017; 2017: 4915736-4915736. 37. Wang Y, Toh HC, Chow P, et al. MicroRNA-224 is up-regulated in hepatocellular carcinoma through epigenetic mechanisms. FASEB J Off Publ Fed Am Soc Exp Biol 2012; 26(7): 3032-3041. 38. Liu F, Liu Y, Shen J, et al. MicroRNA-224 inhibits proliferation and migration of breast cancer cells by down-regulating fizzled 5 expression. Oncotarget 2016; 7(31): 49130-49142. 39. Shi Y, Ye P, Long X. Differential expression profiles of the transcriptome in breast cancer cell lines revealed by next generation sequencing. Cell Physiol Biochem 2017; 44(2): 804-816. 40. Jin Y, Jin W, Zheng Z, et al. GABRB2 plays an important role in the lymph node metastasis of papillary thyroid cancer. Biochem Biophys Res Commun 2017; 492(3): 323-330. 41. Barki M, Xue H. GABRB2, a key player in neuropsychiatric disorders and beyond. Gene 2022; 809: 146021. 42. Vasu S, Shah S, Orjalo A, et al. Novel vertebrate nucleoporins Nup133 and Nup160 play a role in mRNA export. J Cell Biol 2001; 155(3): 339-354. 43. Xu S, Powers MA. Nuclear pore proteins and cancer. Semin Cell Dev Biol 2009; 20(5): 620-630. 44. Shimozono N, Jinnin M, Masuzawa M, et al. Nup160–slc43a3 is a novel recurrent fusion oncogene in angiosarcoma. Cancer Res 2015; 75(21): 4458-4465. 45. Kudela E, Nachajova M, Biringer K, et al. Bilateral ovarian angiosarcoma arising from the mature cystic teratomas – A case report and review of the literature. Int J Surg Case Rep 2018; 42: 90-93. PMMB 2022, 5, 1; a0000269 18 of 19 46. Ishida Y, Otsuka A, Kabashima K. Cutaneous angiosarcoma: update on biology and latest treatment. Curr Opin Oncol 2018; 30(2): 107-112. 47. Mackrill JJ. Ryanodine receptor calcium channels and their partners as drug targets. Biochem Pharmacol 2010; 79(11): 1535-1543. 48. Marchi S, Pinton P. Alterations of calcium homeostasis in cancer cells. Curr Opin Pharmacol 2016; 29: 1-6. 49. Zhang L, Liu Y, Song F, et al. Functional SNP in the microRNA-367 binding site in the 3’UTR of the calcium channel ryanodine receptor gene 3 (RYR3) affects breast cancer risk and calcification. Proc Natl Acad Sci U S A 2011; 108(33): 13653-13658. 50. Abdul M, Ramlal S, Hoosein N. Ryanodine receptor expression correlates with tumor grade in breast cancer. Pathol Oncol Res 2008; 14(2): 157-160. 51. Mariot P, Prevarskaya N, Roudbaraki MM, et al. Evidence of functional ryanodine receptor involved in apoptosis of prostate cancer (LNCaP) cells. The Prostate 2000; 43(3): 205-214. 52. Uhlén M, Hallström BM, Lindskog C, et al. Transcriptomics resources of human tissues and organs. Mol Syst Biol 2016; 12(4): 862. 53. Li L, Li H, Garzel B, et al. SLC13A5 is a novel transcriptional target of the pregnane X receptor and sensitizes drug-induced steatosis in human liver. Mol Pharmacol 2015; 87(4): 674-682. 54. Cai Q, Pan PY, Sheng ZH. Syntabulin-kinesin-1 family member 5B-mediated axonal transport contributes to activity-dependent presynaptic assembly. J Neurosci 2007; 27(27): 7284-7296. 55. Ying Y, Li L, Cao W, et al. The microtubule associated protein syntabulin is required for glucose- stimulated and cAMP-potentiated insulin secretion. FEBS Lett 2012; 586(20): 3674-3680. 56. Su Q, Cai Q, Gerwin C, et al. Syntabulin is a microtubule-associated protein implicated in syntaxin transport in neurons. Nat Cell Biol 2004; 6(10): 941-953. 57. Sin MLY, Mach KE, Sinha R, et al. Deep sequencing of urinary RNAs for bladder cancer molecular diagnostics. Clin Cancer Res 2017; 23(14): 3700-3710. 58. Djamgoz MBA, Coombes RC, Schwab A. Ion transport and cancer: from initiation to metastasis. Philos Trans R Soc B Biol Sci 2014; 369(1638): 20130092. 59. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011; 144(5): 646-674. 60. Stepulak A, Rola R, Polberg K, Ikonomidou C. Glutamate and its receptors in cancer. J Neural Transm Vienna Austria 1996 2014; 121(8): 933-944. 61. Stepulak A, Luksch H, Gebhardt C, et al. Expression of glutamate receptor subunits in human cancers. Histochem Cell Biol 2009; 132(4): 435-445. 62. Luksch H, Uckermann O, Stepulak A, et al. Silencing of selected glutamate receptor subunits modulates cancer growth. Anticancer Res 2011; 31(10): 3181-3192. 63. Chang HJ, Yoo BC, Lim SB, et al. Metabotropic glutamate receptor 4 expression in colorectal carcinoma and its prognostic significance. Clin Cancer Res 2005; 11(9): 3288-3295. PMMB 2022, 5, 1; a0000269 19 of 19 64. Marín YE, Chen S. Involvement of metabotropic glutamate receptor 1, a G protein coupled receptor, in melanoma development. J Mol Med 2004; 82(11): 735-749. 65. Choi CH, Choi JJ, Park YA, et al. Identification of differentially expressed genes according to chemosensitivity in advanced ovarian serous adenocarcinomas: expression of GRIA2 predicts better survival. Br J Cancer 2012; 107(1): 91-99. 66. Fernández-Nogueira P, Bragado P, Almendro V, et al. Differential expression of neurogenes among breast cancer subtypes identifies high risk patients. Oncotarget 2015; 7(5): 5313-5326. 67. Takada M, Tada M, Tamoto E, et al. Prediction of lymph node metastasis by analysis of gene expression profiles in non-small cell lung cancer. J Surg Res 2004; 122(1): 61-69. 68. Ribeiro MPC, Custódio JBA, Santos AE. Ionotropic glutamate receptor antagonists and cancer therapy: time to think out of the box? Cancer Chemother Pharmacol 2017; 79(2): 219-225. 69. Catterall WA. The molecular basis of neuronal excitability. Science 1984; 223(4637): 653-661. 70. Shan B, Dong M, Tang H, et al. Voltage-gated sodium channels were differentially expressed in human normal prostate, benign prostatic hyperplasia and prostate cancer cells. Oncol Lett 2014; 8(1): 345-350. 71. Zhang TT, Han JD, Liu YF. Role of voltage-gated sodium channel in tumor metastasis. J Int Pharm Res 2018; 6: 569-574. 72. Tian Y, Guan Y, Su Y, et al. Trpm2-as promotes bladder cancer by targeting mir-22-3p and regulating gins2 mRNA expression. OncoTargets Ther 2021; 14: 1219-1237. 73. Almasi S, Sterea AM, Fernando W, et al. TRPM2 ion channel promotes gastric cancer migration, invasion and tumor growth through the AKT signaling pathway. Sci Rep 2019; 9: 4182. 74. Orfanelli U, Wenke AK, Doglioni C, et al. Identification of novel sense and antisense transcription at the TRPM2 locus in cancer. Cell Res 2008; 18(11): 1128-1140. 75. Liu X, Cotrim A, Teos L, et al. Loss of TRPM2 function protects against irradiation-induced salivary gland dysfunction. Nat Commun 2013; 4: 1515. 76. Park YR, Chun JN, So I, et al. Data-driven analysis of trp channels in cancer: linking variation in gene expression to clinical significance. Cancer Genomics Proteomics 2016; 13(1): 83-90. 77. Adekola K, Rosen ST, Shanmugam M. Glucose transporters in cancer metabolism. Curr Opin Oncol 2012; 24(6): 650-654. Author(s) shall retain the copyright of their work and grant the Journal/Publisher right for the first publication with the work simultaneously licensed under: Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0). This license allows for the copying, distribution and transmission of the work, provided the correct attribution of the original creator is stated. Adaptation and remixing are also permitted.