Identification of biomarkers associated with cervical lymph node metastasis in papillary thyroid carcinoma: Evidence from an integrated bioinformatic analysis

Abstract

Integrated analysis of accumulated data is an effective way to obtain reliable potential diagnostic molecular of cervical lymph node metastases (LNM) in papillary thyroid carcinoma (PTC). The benefits of prophylactic lymph node dissection (PLND) for these clinically node-negative (cN0) patients remained considerable controversies. Hence, elucidation of the mechanisms of LNM and exploration of potential biomarkers and prognostic indicators are essential for accurate diagnosis of LNM in PTC patients. Up to date, advanced microarray and bioinformatics analysis have advanced an understanding of the molecular mechanisms of disease occurrence and development, which are necessary to explore genetic changes and identify potential diagnostic biomarkers. In present study, we performed a comprehensive analysis of the differential expression, biological functions, and interactions of LNM-related genes. Two publicly available microarray datasets GSE60542 and GSE129562 were available from Gene Expression Omnibus (GEO) database. Differentially expressed genes between clinically node-positive (cN1) and cN0 PTC samples were screened by an integrated analysis of multiple gene expression profile after gene reannotation and batch normalization. Our results identified 48 differentially expressed genes (DEGs) genetically associated with LNM in PTC patients. Gene ontology (GO) analyses revealed the changes in the modules were mostly enriched in the regulation of MHC class II receptor activity, the immune receptor activity, and the peptide antigen binding. Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis of DEGs displayed the intestinal immune network for IgA production, staphylococcus aureus infection, and cell adhesion molecules (CAMs). To screen core genes related to LNM of PTC from the protein-protein interaction network, top 10 hub genes were identified with highest scores. Our results help us understand the exact mechanisms underlying the metastasis of cervical LNM in PTC tissues and pave an avenue for the progress of precise medicine for individual patients.

Keywords

Lymph node metastasis (LNM)papillary thyroid carcinoma (PTC)hub genes protein-protein interaction network (PPI)

1 Introduction

Papillary thyroid carcinoma (PTC) is exceeding malignancies with a rise in novel cases worldwide, which is largely due to the widespread utilization of high-resolution ultrasound and fine-needle aspiration (FNA) biopsy [1 –3]. Thyroidectomy couple with therapeutic neck dissection are preferred for patients with known lymph node metastasis in preoperative examination such as Ultrasound (US), Computed Tomography (CT) and Magnetic Resonance Imaging (MRI). However, the capability of these examination to recognize lymph node metastasis is common restricted for the overlying thyroid gland [4, 5]. Up to date, unfortunately, there remains considerable controversy about the benefits of prophylactic lymph node dissection (PLND) for these clinically node-negative (cN0) patients. Some insist PLND will be helpful to reduce the recurrence of PTC and avoid the possibility of secondary surgery, while others argue PLND increased risk of recurrent laryngeal nerve damage and permanent hypoparathyroidism [6, 7]. Contrast-enhanced ultrasonography (CEUS) is a potential tool to investigate the microvascular flow of different tissues/organs, which displays the microvascularity of cancer tissue and achieves better temporal and spatial resolution than the convention Doppler techniques [8]. Ki-67 is a crucial biomarker that affects cell proliferation and a fundamental indicator for evaluating the recurrence and metastasis of various malignant tumors, including breast cancer, colorectal cancer, lung cancer, and thyroid cancer [9]. Moreover, p53 protein is a tumor suppressor, which can promote the transformation and excessive proliferation of malignant tumor cells [10]. Although extensive efforts have been made, studies providing LNM gene expression profiles of PTC are relatively rare. Thus, it is a pressing task to identify molecular biomarkers and explore underlying mechanisms involved in the LNM in PTC patients, which could help develop novel diagnostic and individual thera-peutic strategies.

Microarray techniques and bioinformatic analysis have been widely adopted to recognize the differentially expressed genes (DEGs) and biological functional pathways [11, 12]. However, solid results may be difficult to be acquired because of small sample size and high false-positive rate in independent microarray analysis. In addition, different microarray platforms, various statistical methods and smaller sample sizes may lead to unconvincing results from these investigations. Accordingly, further comprehensive analysis is needed to identify more reliable biomarkers to overcome these inconsistencies. In our research, two mRNA microarray datasets from the Gene Expression Omnibus (GEO) were picked out for subsequent analysis. DEGs between positive-node (N1) status and negative-node (N0) status from patients with PTC were screened to determine key biomarkers. Subsequently, possible mechanisms underlying metastasis of lymph node were investigated through Ontology (GO) enrichment, Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment and protein-protein interaction (PPI) network analysis. In summary, a total of 48 DEGs and top 10 hub genes were identified, which may be prospective diagnosis and therapy biomarkers and therapeutic targets for LNM in PTC patients, specifically in deciding to perform PLND in the absence of clinical LNM.

2 Materials and methods

Our research was conducted according to Declaration of Helsinki and approved by approved by the Ethics Committee of the Affiliated Hospital of Jiangsu University.

2.1 Microarray data

GEO (http://www.ncbi.nlm.nih.gov/geo) is a public functional genomics data repository of high throughout gene expression data, chips and microarrays [13]. Two gene expression datasets GSE60542 and GSE129562 were downloaded from GEO. The probes were converted into the corresponding gene symbol according to the annotation information. The GSE60542 contained 19 positive-node PTC samples and 14 negative-node PTC samples. GSE129562 is involved in 5 present-node tissues and 3 absent-node tissues in PTC.

2.2 Identification of DEGs

GEO datasets were often based on different platforms, and samples were processed on different days, by different groups even by different people, we integrated all samples from the two datasets to improve the number of samples (21 positive-node PTC samples and 17 negative-node PTC samples). Batch normalization was applied to provide a balance between discovery of statistically significant genes and limitations of false-positives using sva and limma R package. Probe sets without corresponding gene symbols or genes with more than one probe set were removed or averaged, respectively. DEGs were screened by R package limma (http://bioconductor.org/packages/release/bioc/html/limma.html) [14]. Genes from each sample with the following criteria were retained: (1) a |log₂ (foldchange)|>0.5 and (2) P < 0.05. The integrated dysregulated gene lists were saved for following analysis.

2.3 Kyoto encyclopedia of genes and genomes and gene ontology enrichment analyses of DEGs

The identified DEGs were annotated by Gene Ontology and Kyoto encyclopedia of genes and genomes (KEGG) analysis using “clusterProfiler” Bioconductor package [15], which integrates biological data and analysis tools, and provides a comprehensive set of functional annotation of genes to extract biological information. GO is a major bioinformatics tool to annotate genes and analyze biological process of these genes [16]. KEGG is used to understanding high-level functions and biological systems from large-scale molecular datasets generated by high-throughput experimental technologies [17]. P value <0.05 was considered as statistical significance.

2.4 Protein-protein interaction network construction and module analysis

Protein-protein interaction (PPI) networks were predicted using the Search Tool for the Retrieval of Interacting Gene (STRING) online database (http://string-db.org; version 10.5) [18]. Each node in the PPI network represents a gene, protein or other molecule, and the connections between the nodes represent the interactions between these biomolecules. The most closely associated nodes may indicate hub genes with important physiological regulatory functions [19]. The analysis of functional interactions between proteins may provide insight into the mechanism of the occurrence or metastasis of the disease. In current investigation, the PPI network of DEGs was constructed via string database, and the interactions with comprehensive scores>0.4 was considered statistically significant.

2.5 Statistical analysis

Quantitative data are exhibited as mean±standard deviation (SD). Student’s t-test or one-way analysis of variance (ANOVA) were used for statistical analysis of all cell line experiments. Gene expression correlations were assessed using Spearman’s correlation. P-value <0.05 was considered statistically significant.

3 Results

3.1 Convergence of gene expression signatures across different experiments from GEO database

In present investigation, two datasets GSE60542 and GSE129562 were available from GEO database. The GSE60542 dataset included 19 positive-LNM and 14 negative-LNM in PTC samples. The GSE129562 dataset consisted of 5 present-LNM tissues and 3 absent-LNM in PTC tissues. The data from the two microarray datasets were batch normalized and then merged to increase sample number. Next, limma package was applied to screen differentially expressed genes according to the threshold of |log₂ (foldchange)| and P < 0.05. Total 48 DEGs were finally selected out, comprising 14 up-regulated genes and 34 down-regulated genes (Table 1). Volcano plot and heat maps for DEGs, illustrating the trends in DEG expression, were visualized in Fig. 1.

Table 1
All 48 differentially expressed genes (DEGs) were screen from two profile datasets

DEGs Genes name

Up-regulated RASD1, ELOVL6, CDO1, IGFBP3, PXDN, TM4SF1, SPOCK1,

LPL, ESM1, SLPI, IRX2, FAM155B, LRCH2, TMOD1.

Down-regulated IL15, HLA-DQA1, CBLN4, CX3CR1, TPST2, AQP1, NPFFR2,

CXCL17, COL15A1, SELE, ADRA2A, ADM, PON3, PLN,

PLXNC1, SELL, GEM, FAM46C, WISP2, FCGR2B, HLA-DPA1,

BHMT2, LAX1, WARS, HCP5, LRP8, CCL8, QPCT, CD38,

CD79A, HLA-DRA, HLA-DMA, SERPINA1, SHISA2.

DEGs	Genes name
Up-regulated	RASD1, ELOVL6, CDO1, IGFBP3, PXDN, TM4SF1, SPOCK1,
	LPL, ESM1, SLPI, IRX2, FAM155B, LRCH2, TMOD1.
Down-regulated	IL15, HLA-DQA1, CBLN4, CX3CR1, TPST2, AQP1, NPFFR2,
	CXCL17, COL15A1, SELE, ADRA2A, ADM, PON3, PLN,
	PLXNC1, SELL, GEM, FAM46C, WISP2, FCGR2B, HLA-DPA1,
	BHMT2, LAX1, WARS, HCP5, LRP8, CCL8, QPCT, CD38,
	CD79A, HLA-DRA, HLA-DMA, SERPINA1, SHISA2.

Fig. 1

(a, b) Volcano of DEGs among GSE60542 and GSE129562 after batch normalization. (b) Heatmap of 48 DEGs available from GEO database. The red color represents the upregulated genes, the green color represents the downregulated genes, and the black color represents genes without any change.

3.2 Gene Ontology (GO) functional and signaling pathway enrichment of DEGs

We then performed GO enrichment analysis of DEGs in LNM of PTC samples using “clusterProfiler” Bioconductor package. The GO enriched function of DEGs were mainly involved in immune receptor activity (5 genes), G protein-coupled receptor binding (4 genes), and MHC class II receptor activity (3 genes) (Table 2; Fig. 2a). Furthermore, 25 KEGG pathways were overrepresented in the DEGs. Most of the DEGs were predominantly enriched in intestinal immune network for IgA production (6 genes), staphylococcus aureus infection (5 genes), and cell adhesion molecules (CAMs) (5 genes) (Table 3; Fig. 2b).

Table 2
Gene ontology analysis of differentially expressed genes of cervical LNM in PTC tissues

Term Description Gene count, n P-value

GO:0140375 Immune receptor activity 5 1.11E-05

GO:0001664 G protein-coupled receptor binding 4 3.92E-03

GO:0032395 MHC class II receptor activity 3 1.37E-06

GO:0042605 Peptide antigen binding 3 4.96E-05

O:0004867 Serine-type endopeptidase inhibitor activity 3 1.34E-03

GO:0019838 Growth factor binding 3 3.90E-03

GO:0070492 Oligosaccharide binding 2 5.40E-04

GO:0023026 MHC class II protein complex binding 2 6.16E-04

GO:0034185 Apolipoprotein binding 2 6.97E-04

GO:0023023 MHC protein complex binding 2 1.52E-03

GO:0005520 Insulin-like growth factor binding 2 1.91E-03

GO:0071813 Lipoprotein particle binding 2 2.64E-03

GO:0071814 Protein-lipid complex binding 2 2.64E-03

Term	Description	Gene count, n	P-value
GO:0140375	Immune receptor activity	5	1.11E-05
GO:0001664	G protein-coupled receptor binding	4	3.92E-03
GO:0032395	MHC class II receptor activity	3	1.37E-06
GO:0042605	Peptide antigen binding	3	4.96E-05
O:0004867	Serine-type endopeptidase inhibitor activity	3	1.34E-03
GO:0019838	Growth factor binding	3	3.90E-03
GO:0070492	Oligosaccharide binding	2	5.40E-04
GO:0023026	MHC class II protein complex binding	2	6.16E-04
GO:0034185	Apolipoprotein binding	2	6.97E-04
GO:0023023	MHC protein complex binding	2	1.52E-03
GO:0005520	Insulin-like growth factor binding	2	1.91E-03
GO:0071813	Lipoprotein particle binding	2	2.64E-03
GO:0071814	Protein-lipid complex binding	2	2.64E-03

Fig. 2

(a) Gene Ontology analysis was used to obtain significant functional analysis. (b) KEGG analysis were used to obtain significant enriched KEGG terms. The color represents the P value.

Table 3

Kyoto encyclopedia of genes and genomes pathways of differentially expressed genes of cervical LNM in PTC tissues

ID	Term	Gene count, n	P-value
hsa04514	Cell adhesion molecules (CAMs)	6	7.61E-06
hsa04672	Intestinal immune network for IgA production	5	5.01E-07
hsa05323	Rheumatoid arthritis	5	1.23E-05
hsa05150	Staphylococcus aureus infection	5	1.44E-05
hsa04640	Hematopoietic cell lineage	5	1.67E-05
hsa05310	Asthma	4	3.00E-06
hsa05330	Allograft rejection	4	6.92E-06
hsa05332	Graft-versus-host disease	4	9.43E-06
hsa04940	Type I diabetes mellitus	4	1.14E-05
hsa05320	Autoimmune thyroid disease	4	2.65E-05
hsa05416	Viral myocarditis	4	4.35E-05
hsa05321	Inflammatory bowel disease (IBD)	4	5.96E-05

3.3 Protein-protein interaction (PPI) network creation and hub gene identification

All DEGs were screened using the STRING database to further investigate their properties and the interactions among them. The PPI network was constructed based on STRING database to provide deep insights into the mechanisms underlying metastasis of cervical LNM in PTC tissues. The PPI network consists of 20 nodes and 31 edges, in which the nodes correspond to genes and the edges represent the links between genes, and top 10 genes with the most connections in PPI relationship pairs were graphed in barplot diagram, including SELL, CX3CR1, FCGR2B, CD38, IL15, HLA-DRA, CCL8, CD79A, HLA-DMA and HLA-DPA1 (Table 4; Fig. 3).

Table 4
Functional roles of top 10 hub genes

No. Gene symbol Full name Function

1 SELL Selectin L Selectin L facilitates the migration into secondary lymphoid organs and inflammation sites.

2 CX3CR1 C-X3-C Motif Chemokine Receptor 1 The protein encoded by CX3CR1 is a receptor for fractalkine.

3 FCGR2B Fc Fragment of IgG Receptor IIb The encoded protein by FCGR2B is involved in the phagocytosis of immune complexes and in the regulation of antibody production by B-cells.

4 CD38 CD38 Molecule High expression of CD38 is associated with Richter’s Syndrome and Prolymphocytic Leukemia.

5 IL15 Interleukin 15 IL15 regulates T and natural killer cell activation and proliferation

6 HLA-DRA HLA Class II Histocompatibility Syndrome Antigen HLA-DRA is related to Graham-Little-Piccardi-Lassueur and Heart Lymphoma.

7 CCL8 C-C Motif Chemokine Ligand 8 CCL8 may contribute to tumor-associated leukocyte infiltration by recruiting leukocytes to sites of inflammation

8 CD79A CD79a Molecule CD79A is overexpressed in Agammaglobulinemia 3, Autosomal Recessive and Agammaglobulinemia, Non-Bruton Type.

9 HLA-DMA HLA Class II Histocompatibility Sarcoidosis Antigen HLA-DMA has been implicated as a regulator of Cardiac and Measles.

10 HLA-DPA1 HLA Class II Histocompatibility Antigen HLA-DPA1 plays a regulatory role in cytokinesis, cell growth and differentiation.

No.	Gene symbol	Full name	Function
1	SELL	Selectin L	Selectin L facilitates the migration into secondary lymphoid organs and inflammation sites.
2	CX3CR1	C-X3-C Motif Chemokine Receptor 1	The protein encoded by CX3CR1 is a receptor for fractalkine.
3	FCGR2B	Fc Fragment of IgG Receptor IIb	The encoded protein by FCGR2B is involved in the phagocytosis of immune complexes and in the regulation of antibody production by B-cells.
4	CD38	CD38 Molecule	High expression of CD38 is associated with Richter’s Syndrome and Prolymphocytic Leukemia.
5	IL15	Interleukin 15	IL15 regulates T and natural killer cell activation and proliferation
6	HLA-DRA	HLA Class II Histocompatibility Syndrome Antigen	HLA-DRA is related to Graham-Little-Piccardi-Lassueur and Heart Lymphoma.
7	CCL8	C-C Motif Chemokine Ligand 8	CCL8 may contribute to tumor-associated leukocyte infiltration by recruiting leukocytes to sites of inflammation
8	CD79A	CD79a Molecule	CD79A is overexpressed in Agammaglobulinemia 3, Autosomal Recessive and Agammaglobulinemia, Non-Bruton Type.
9	HLA-DMA	HLA Class II Histocompatibility Sarcoidosis Antigen	HLA-DMA has been implicated as a regulator of Cardiac and Measles.
10	HLA-DPA1	HLA Class II Histocompatibility Antigen	HLA-DPA1 plays a regulatory role in cytokinesis, cell growth and differentiation.

Fig. 3

(a) Protein-protein interaction (PPI) network constructed among DEGs with the criterion of interaction score >0.4. (b) Top 10 genes in PPI relationship pairs were graphed in barplot diagram. The x-axis represents the number of connections.

4 Discussion

PTC is associated with an excellent clinical outcome, as a result, the 5-year survival rate for PTC is approximately 98% and is one of the indolent diseases among cancers [20, 21]. Although the mortality rate is very low, recurrent and metastatic diseases can occur in approximately 20–30% of cases, which are mainly related to metastatic lymph nodes [22]. The role of routine lymph node dissection in the treatment for PTC remains controversial. In the preset, contrast-enhanced ultrasound (CEUS) is a common method to detect LNM in clinical practice, which usually presents centripetal or asynchronous perfusion, hyperenhancement, heterogeneous enhancement, perfusion defects and ring-enhancing margins in these metastasis LNM in PTC patients [23]. The quantitative CEUS has the higher sensitivity to detect the subtle hemodynamic changes of the lesions because metastatic LNMs often secret vascular growth factors to form lots of twisted, irregular diameter and incomplete vascular walls of immature micro-vessels [24]. However, their value in clinical application is still controversial attributed to inconsistent and unsatisfactory results [25, 26]. It is generally accepted therapeutic neck dissection is recommended to remove known metastatic lymph nodes identified either pre- or intraoperatively because this surgical management can reduce the incidence of PTC persistence/recurrence [27]. However, there are no randomized controlled trails to support the benefit of routine prophylactic LND for PTC patients with no pre-or intraoperative evidence of nodal metastasis [28, 29]. Therefore, due to the debates of current surgical treatment of PTC and prophylactic LND as well as the fuzziness in the literature, molecular biomarkers that allow risk stratification of the presence of LNM will be of great clinical value. US-guided fine-needle aspiration biopsy (FNAB) is widespread utilized to improve the diagnostic accuracy through cytologic and genomic detection [30]. More importantly, the capability to genomic detection in cytologic specimens is not subordinate to that in pathologic specimens [31].

Microarray is a promising and popular method for large-scale gene expression profiling, greatly facilitating the analysis of thousands of dysregulated key genes responsible for carcinogenesis [12 , 32–34]. In the present study, we identified biomarkers involved in metastatic lymph nodes of PTC patients by analyzing gene expression data downloaded from GEO databases. Since the datasets collected for two-dataset analysis were based on different platforms, and the samples were processed on different dates or in different groups, we used the R package to remove or minimize batch effects [35]. A total of 48 DEGs were identified to associate with LNM between clinically node-positive (cN1) and node-negative (cN0) PTC samples. GO annotation result of DEGs demonstrated immune receptor activity, G protein-coupled receptor binding, and MHC class II receptor activity were chiefly enriched. As the largest family of membrane protein receptors, G protein coupled receptors can be activated by a variety of ligands and play corresponding signal transduction functions, thus participating in crucial physiological processes in organisms, such as the regulation of the immune system, the regulation of the autonomic nervous system and the maintenance of homeostasis [36]. Upon G protein coupled receptor is stimulated, G protein related signal pathways including cAMP pathway, phosphatidylinositol signal pathway, and MAPK/Erk signal pathway [37 –39], which will largely affect cell survival and cause adverse stress reactions. Taken together, these observations imply that the activation of G protein coupled receptor may be a major contributor to metastatic lymph nodes. KEGG enrichment analysis of DEGs showed many of these genes were mapped to intestinal immune network for IgA production, staphylococcus aureus infection, and cell adhesion molecules (CAMs), which suggested a critical role of immune and inflammatory responses in metastatic of lymph nodes. The directional migration of lymphocytes is caused by the interaction between the lymphocyte homing receptor and vascular address on endothelial cells [40, 41]. In the PPI network identified in the present study, top 10 hub genes were highlighted as the most significant hub genes, with multiple interactions in the network.

5 Limitations

It is undeniable that there are still several limitations in our research. All DEGs retrieved from online databases, so this result has not been verified by clinical research. Also, although we integrated two GEO datasets to enlarge sample number, the total sample number was still not big enough to restrict the discovery of more DEGs in metastasis of lymph node. Therefore, further validation studies could lead to additional insights into the development of LNM as well as the identification of additional functional biomarker candidates for improved clinical diagnostic of LNM in PTC patients.

6 Conclusion

In summary, our comprehensive bioinformatic analysis has identified numerous useful molecular targets for the future investigation of the mechanisms and selection of biomarkers for LNM in PTC patients. Several important biological processes and pathways as well as the hub genes functioning in these processes may provide novel insights into the development and progression of metastasis of lymph node. Furthermore, further molecular biological exploration will be performed to confirm the function of the identified genes associated with LNM in PTC patients. Last but not the least, pre-operative hub genes analysis via minimally invasive US-FNA biopsy can guide the surgical management of PTC patients.

Footnotes

Acknowledgments

We greatly acknowledged the Zhenjiang Social Development Fund (SH2018035), Zhenjiang Key Research and Development Plan Fund (SH2016036) and the Doctoral Start-up Fund of Affiliated Hospital of Jiangsu University (jdfyRC2017013).

Disclosure

All authors report no conflicts of interest in this work.

References

, et al. Association of BRAF V600E Mutation and MicroRNA Expression with Central Lymph Node Metastases in Papillary Thyroid Cancer: A Prospective Study from Four Endocrine Surgery Centers. Thyroid: Official Journal of the American Thyroid Association. 2016;26(4):532–42.

Chen

, et al. A multivariable model of BRAFV600E and ultrasonographic features for predicting the risk of central lymph node metastasis in cN0 papillary thyroid microcarcinoma. Cancer Management & Research. 2019;11:7211–7.

Chen

, et al. Association of BRAFV600E mutation with ultrasonographic features and clinicopathologic characteristics of papillary thyroid microcarcinoma: A retrospective study of 116 cases. Clinical Hemorheology and Microcirculation. 2019;73(1):1–8.

Lee

, et al. Clinical characteristics of papillary thyroid microcarcinoma less than or equal to 5 mm on ultrasonography. European Archives of Oto Rhino Laryngology. 2013;270(11).

Liu

, et al. Conventional ultrasound characteristics, TI-RADS category and shear wave speed measurement between follicular adenoma and follicular thyroid carcinoma. Clinical Hemorheology and Microcirculation. 2020(6):1–11.

David , et al. Revised American Thyroid Association management guidelines for patients with thyroid nodules and differentiated thyroid cancer. Thyroid : Official Journal of the American Thyroid Association, 2009.

, GM

, DL

. The pros and cons of prophylactic central compartment lymph node dissection for papillary thyroid carcinoma. Thyroid: Official Journal of the American Thyroid Association. 2009;19(7):683–9.

Luo

, et al. Utility of quantitative contrast-enhanced ultrasound for the prediction of lymph node metastasis in patients with papillary thyroid carcinoma. Clinical Hemorheology and Microcirculation, 2020. doi: 10.3233/CH-200909

Yao

, et al. , Value of Ultrasound Combined with Immunohistochemistry Evaluation of Central Lymph Node Metastasis for the Prognosis of Papillary Thyroid Carcinoma. Cancer Management and Research. 2020;12:8787–99.

10.

Khalesi

, et al. Comparison of PTCH1, COX-2, p53, and Ki-67 protein expression in basal cell carcinomas of nodular and superficial subtypes arising on the head and trunk. International Journal of Dermatology, 2016.

11.

L L, et al. Screening and identification of key biomarkers in hepatocellular carcinoma: Evidence from bioinformatic analysis. Oncology Reports. 2017;38(5):2607–18.

12.

Comprehensive and Integrative Genomic Characterization of Hepatocellular Carcinoma. Cell. 2017;169(7):págs. 1327–1341.

13.

T B, et al. NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Research. 2013;41:D991–5.

14.

, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 2015;43(7):e47.

15.

G Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics: A Journal of Integrative Biology. 2012;16(5):284–7.

16.

Ashburner

, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nature Genetics. 2000;25(1):25–9.

17.

Kanehisa

, The KEGG database. Novartis Foundation Symposium. 2002;247(247):91.

18.

D S, et al. STRING v protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Research. 2015;43:D447–52.

19.

Zhang

, et al. Investigation of hypoxia networks in ovarian cancer via bioinformatics analysis. Journal of Ovarian Research. 2018;11(1):16.

20.

Chen

, et al. The diagnostic performances of conventional strain elastography (SE), acoustic radiation force impulse (ARFI) imaging and point shear-wave speed (pSWS) measurement for non-calcified thyroid nodules. Clinical Hemorheology and Microcirculation. 2017;65(3):259–73.

21.

Zhao

, et al. Factors associated with initial incomplete ablation for benign thyroid nodules after radiofrequency ablation: First results of CEUS evaluation. Clin Hemorheol Microcirc. 2017:1-13.

22.

Brown

, de

Souza JA

, Cohen

. Thyroid Cancer: Burden of Illness and Management of Disease. Journal of Cancer. 2011;2:193–9.

23.

Hong

, et al. Role of Contrast-Enhanced Ultrasound in the Pre-operative Diagnosis of Cervical Lymph Node Metastasis in Patients with Papillary Thyroid Carcinoma. Ultrasound in Medicine & Biology. 2017;43(11):2567–75.

24.

Nakopoulou

, et al. Expression of the vascular endothelial growth factor receptor-2/Flk-1 in breast carcinomas: Correlation with proliferation. Human Pathology. 2002;33(9):863–70.

25.

Dudau

, et al.

Can contrast-enhanced ultrasound distinguish malignant from reactive lymph nodes in patients with head and neck cancers?

Ultrasound in Medicine & Biology. 2014;40(4):747–54.

26.

Poanta

, et al. The place of CEUS in distinguishing benign from malignant cervical lymph nodes: a prospective study. Medical Ultrasonography. 2014;16(1):7–14.

27.

S D, et al. Revised American Thyroid Association management guidelines for patients with thyroid nodules and differentiated thyroid cancer. Thyroid : Official Journal of the American Thyroid Association, 2009.

28.

Calò

, et al.

Total thyroidectomy without prophylactic central neck dissection in clinically node-negative papillary thyroid cancer: is it an adequate treatment? World Journal of Surgical Oncology

12(1):152–152.

29.

Mazzaferri

, Doherty

, Steward

. The Pros and Cons of Prophylactic Central Compartment Lymph Node Dissection for Papillary Thyroid Carcinoma. Thyroid. 2009;19(7):683–9.

30.

Izquierdo

, et al. Ultrasound-guided fine-needle aspiration in the management of thyroid nodules in children and adolescents. Thyroid Official Journal of the American Thyroid Association. 2009;19(7):703.

31.

Kwak

, et al. Association of BRAFV600E mutation with poor clinical prognostic factors and US features in Korean patients with papillary thyroid microcarcinoma. Radiology. 2009;253(3):854–60.

32.

Wang

, et al. Identification of diagnostic markers for major depressive disorder by cross-validation of data from whole blood samples. PeerJ. 2019;7.

33.

, et al. Integrative Analysis of Siglec-15 mRNA in Human Cancers Based on Data Mining. Journal of Cancer. 2020;11(9):2453–64.

34.

Sang

, et al. Identification and validation of hub microRNAs dysregulated in esophageal squamous cell carcinoma. Aging. 2020;12(10).

35.

Holger

, et al. Batch-normalization of cerebellar and medulloblastoma gene expression datasets utilizing empirically defined negative control genes. Bioinformatics. 2019(18):18.

36.

, PH

. GPCRs: Emerging anti-cancer drug targets. Cellular Signalling. 2018;41:65–74.

37.

Rodríguez

, Setaluri

. Cyclic AMP (cAMP) signaling in melanocytes and melanoma. Archives of Biochemistry & Biophysics. 2014;563:22–7.

38.

Archbold

, Martin

, Sweet

. Towards selective lysophospholipid GPCR modulators. Trends in Pharmacological Sciences. 2014;35(5):219–26.

39.

Ishii

, Shimizu

. Platelet-activating factor (PAF) receptor and genetically engineered PAF receptor mutant mice. Progress in Lipid Research. 2000;39(1):41–82.

40.

Pandolfi

, et al. Integrins: Integrating the Biology and Therapy of Cell–Cell Interactions. Clinical Therapeutics. 2017:S0149291817310640.

41.

Mashud

, et al. Impaired lymphocyte trafficking in mice deficient in the kinase activity of PKN1. Sci Rep. 2017;7(1):7663.