Identification and Validation of Novel Genes in Anaplastic Thyroid Carcinoma via Bioinformatics Analysis

Introduction

Thyroid carcinoma is the most common endocrine malignancy,¹ which has become the second most common malignant tumor of female endocrine systems.² Anaplastic thyroid carcinoma (ATC) is rare (2% of thyroid carcinoma), but it is associated with rapid progression, high invasiveness and high mortality.³ The treatments of thyroid carcinoma are mainly compose of surgery, drugs, and radioactive iodine ¹³¹I (RAI). These are effective for differentiated thyroid carcinoma (DTC), however, less effective, even no effective for ATC.^4,5 Because ATC does not have the ability to uptake iodine, nor regulated by thyroid stimulating hormone (TSH), RAI and TSH inhibiting therapy are also not effective for ATC.^6,7 Up to now, there is no available established treatment that prolongs survival. In recent decades, with the understanding and research of the individual diversity of tumors, genetic analysis has been used as a more conducive method for accurate diagnosis and effective treatment. Therefore, exploring novel molecular biomarkers of ATC will help to find new therapeutic methods.

Microarray technology has been widely used to profile gene expression and explore the potential pathogenesis of cancers and is suitable for the screening of the differentially expressed genes (DEGs) in recent years.^8–10 Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) was used for the bioinformatics data mining of gene expression profiles.¹¹ Microarray technology and bioinformatics analysis methods provide new ideas for the study on molecular mechanisms of diseases and are the best ways to explore gene regulation during disease development currently.¹² Combining microarray technology with bioinformatics analysis allows for comprehensive analysis of the DEGs, provides a new and effective way to explore pivotal pathways and molecular mechanisms of ATC.

In the present study, we identified overlapping DEGs of ATC via analyzing five microarray datasets downloaded from the GEO database. The aim of the study is to find DEGs and explore new potential therapeutic targets for ATC. We screened the novel genes in ATC via a bioinformatics analysis and preliminarily validated these genes through in vitro experiments. Our study results indicated that CENPF, KIF2C, and TPX2 might be novel biomarkers for the treatment of ATC.

Materials and Methods

ATC Dataset Preprocessing

The microarray data GSE9115, GSE29265, GSE33630, GSE53072, and GSE65144 of ATC were downloaded from GEO database in the National Center for Biotechnology Information (NCBI) by searching the following key terms: ((anaplastic thyroid cancer) OR anaplastic thyroid carcinoma) OR anaplastic thyroid tumor) OR anaplastic thyroid neoplasm) OR ATC) AND “Homo sapiens”. All included chip datasets have the same characteristics, that is, the comparison of ATC tissues and normal specimens, and contained at least 6 corresponding samples and all microarray data available to the public.

The GSE9115 dataset included five ATC tissues and four normal thyroid tissues (noncancerous samples, normal thyroid specimens). The GSE29265 consisted of nine ATC tissues and ten adjacent normal thyroid tissues. The GSE33630 contained eleven ATC tissues and forty-five non-tumoral thyroid tissues. The GSE53072 included five ATC tissues and four normal thyroid tissues. The GSE65144 consisted of twelve ATC tissues and thirteen normal thyroid tissues.

Data Processing

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r) is an interactive web tool that allows users to compare different groups of Samples in a GEO Series.¹¹ Literature reviewing showed that most similar articles^13–15 used GEO2R to analyze the data from different platforms. In our study, GEO2R was used to screen the differentially expressed genes (DEGs) between ATC and normal thyroid samples in five datasets, respectively. The selection criteria of DEGs were based on published literatures.^16–18 Cao’s study used |log FC|≥1 and the adj. P < 0.05 as the selection criteria of DEGs and DE miRNAs selection.¹⁶ Chen’s study used log FC (fold change) >1 and adj. P < 0.01 as the selection criteria of DEGs.¹⁷ The screening conditions for the DEGs were absolute value of p < 0.05 and |log FC|>1 in Yan’s study.¹⁸ Therefore, compared with the normal thyroid specimens, the GEO2R was conducted to screen DEGs in each chip under the condition of |log FC| > l, adjusted P‐values (adj. P) <0.05. VENNY online tool (http://bioinfogp.cnb.csic.es/tools/venny) was performed to obtain the overlapping DEGs in the five datasets.

Functional Enrichment Analysis of DEGs

The overlapping DEGs were uploaded to the Database for Annotation Visualization and Integrated Discovery (DAVID, version 6.8, https://david.ncifcrf.gov) for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses.^19,20 In addition, the cut-off value for GO and KEGG analyses were set to P < 0.05.

Protein–Protein Interaction Network Construction

The Retrieval of Interacting Genes (STRING; http://stringdb.org) database was used to calculate protein–protein interaction (PPI) networks of DEGs with a combined score >0.4 as the cut-off criteria.^21,22 The PPI network was visualized by Cytoscape software (version 3.2.1).²³ Subsequently, the Cytoscape plugin Molecular Complex Detection (MCODE) was utilized to identify the hub genes and critical modules in the PPI network with the default parameters (Degree cut-off 2, Node score cut-off 2, K-core 2 and Max depth = 100).²⁴ The Cytoscape plugin CytoHubba provided a simple interface to analyze a network with eleven scoring methods.²⁵ Among the 11 methods, the newly proposed method Maximal Clique Centrality (MCC) has a better performance in the PPI network. In addition, the crucial nodes were identified based on Degree Centrality (DC), nodes with a high degree were identified as hub genes.²⁶

Therefore, the hub genes in the PPI network were selected according to the screening conditions: i) The PPI network model was analyzed by Molecular Complex Detection plugin (MCODE). In the result of default condition analysis, target points corresponding to the model with score >3 and number of nodes >5 were selected; ii) The top 50 targets of MCC algorithm in CytoHubba; iii) The target points larger than average degree value two times. Finally, the hierarchical clustering analysis was performed for the hub genes.

Survival Analysis of the Hub Genes

The publically available database Kaplan–Meier plotter (KM plotter, http://kmplot.com/analysis/) was utilized to analyze the prognostic value of the hub genes.²⁷ In order to analyze the overall survival (OS) of patients with ATC, the patients were divided into two groups, the high expression group and the low gene expression group (control group) based on the median expression of the selected gene. In addition, the hazard ratio (HR) with 95% confidence intervals (CIs) and log-rank P value were calculated and presented on the plot, and the log-rank P value <0.05 was considered statistically significant.

Cell Culture

Human normal thyroid follicular epithelial cell line (Nthy-ori 3–1 cell) and human anaplastic differentiated thyroid cancer cell lines (ARO cell, CAL-62 cell) were employed in our study. Nthy-ori 3–1 cells and ARO cells were obtained from the First Affiliated Hospital of Shandong First Medical University. CAL-62 cells were purchased from Procell (Wuhan, China). As shown in Supplementary Figure S1, thyroglobulin (TG) expression was detected by immunofluorescence in ARO cells. TG is a specific marker protein of thyroid tumor cells. Therefore, ARO cell line is a thyroid-derived tumor cell line. As shown in the Report of Cell Line Identification, we have done STR profile identification of ARO cells. Therefore, ARO cell line is the cells of human anaplastic thyroid carcinoma. The use of the Nthy-ori 3–1 cells, ARO cells, and CAL-62 cells were approved by the institutional research ethics committee of the First Affiliated Hospital of Shandong First Medical University. Cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco Life Technology, Grand Island, NY, USA) containing 10% fetal bovine serum (FBS) and 1% streptomycin at 37°C with 5% CO₂. In addition, cells with 80–90% adherence were used for subculturing and were respectively harvested with trypsin (Gibco Life Technology, Grand Island, NY, USA).

Quantitative RT-PCR (RT‐qPCR)

Total RNA was respectively extracted from Nthy-ori 3–1 cells, ARO cells and CAL-62 cells using RNA TriZol reagent (TaKaRa Bio, Otsu, Japan) according to the manufacturer’s protocol. The cDNA was synthesized using the First Strand cDNA Synthesis Kit (TaKaRa Bio, Otsu, Japan). Quantitative real‐time polymerase chain reaction (RT‐qPCR) was performed with the Ultra SYBR One Step RT-qPCR kit (CWBIO, Beijing, China) in a total volume of 20μL. The primers for ANLN, CENPF, KIF2C, TPX2, NDC80, and GAPDH were designed (Table 1) and synthesized by TaKaRa Bio (Otsu, Japan). The expression level of each target gene was calculated using the 2^−ΔΔCt method and GAPDH was used as the internal control.

Table 1 Primer Sequences Used to Amplify Target Genes by Real-Time RT-PCR

Statistical Analysis

Statistical analysis was performed using SPSS 21.0 statistical software (SPSS, Chicago, IL, USA). Each experiment was done at least three independent measurements. All the data were presented as means ± standard deviation (SD). Student’s t-test was used to analyze statistical significance. P<0.05 was considered as significant (*: P<0.05, **: P<0.01, *** P<0.001).

Results

Identification of DEGs

A total of 1388, 2548, 4796, 3004, and 4644 DEGs from GSE9115, GSE29265, GSE33630, GSE53072, and GSE65144 datasets were identified after standardization of the microarray results. As shown in the volcano diagram (Figure 1), red represents upregulated genes and blue downregulated genes. There were 293 overlapped DEGs among the five datasets (Figure 1F). After removing nine genes whose log FC values were inconsistent in the five datasets, 284 overlapped DEGs remained, which included 121 upregulated genes and 163 downregulated genes (Supplementary Table S1). The flowchart of the work was shown in Supplementary Figure S2.

Figure 1 Volcano plot of gene expression profile data in ATC and normal samples. (A) Volcano plot of GSE9115, (B) Volcano plot of GSE29265, (C) Volcano plot of GSE33630, (D) Volcano plot of GSE53072, and (E) Volcano plot of GSE65144. The red color represents the upregulation, and the green color represents the downregulation. (F) Venn diagram of the overlapping DEGs.

Function and Pathway Enrichment Analyses of DEGs

GO and KEGG analyses were performed using the DAVID to investigate the biological function of 284 overlapped DEGs between ATC and normal thyroid tissues. GO analysis showed that the biological processes (BP) terms of the DEGs were significantly enriched in cell division, mitotic nuclear division, metaphase plate congression, mitotic sister chromatid segregation, etc. The molecular function (MF) terms also indicated that the DEGs were significantly enriched in ATP binding, signal transducer activity, protein binding, etc. In the cell component (CC) ontology, the DEGs were mainly enriched in kinetochore, extracellular exosome, proteinaceous extracellular matrix, spindle midzone, etc. Figure 2 shows the top 10 enriched functions of BP, CC, and MF, respectively. In addition, the DEGs were mainly enriched in seven KEGG pathways (p<0.05), which included Cell cycle, Fanconi anemia pathway, Oocyte meiosis, DNA replication, Central carbon metabolism in cancer and HIF-1 signaling pathway (Figure 2D).

Figure 2 Functional and pathway enrichment analyses of the overlapping DEGs in ATC. (A) Biological Process (BP), (B) Cellular Component (CC), (C) Molecular Function (MF), (D) and KEGG pathway analyses. The x-axis represents the Q value, and the y-axis represents the GO term. The GO terms were measured by the rich factor, Q value, and number of genes enriched.

Protein–Protein Interaction (PPI) Network Construction

The PPI networks of 284 DEGs were analyzed by using the String database, and 232 DEGs with a combined score >0.4 were imported into the Cytoscape software for further analysis. As shown in Figure 3A, the PPI network consisted of 232 nodes interacting via 1674 edges. The hub genes were obtained according to the screening criteria: i) The PPI network model was analyzed by MCODE. Among the results analyzed by default, there was one model with model score >3 and node number >5 (model score =47.16 and node number =52), including 52 genes; ii) Select the top 50 targets of MCC algorithm in CytoHubba; iii) The average degree value was 14.43103448, and 53 genes were more than twice the average degree. Finally, 50 genes that satisfied all three conditions were regarded as hub genes. These hub genes were CCNA2, TYMS, CCNB2, AURKA, TPX2, KIF11, KIF2C, MAD2L1, BUB1B, CDC6, CDCA8, UBE2C, PBK, CENPF, RAD51, NUSAP1, MKI67, CENPE, TRIP13, BUB1, MCM4, TTK, NCAPH, KIF23, CEP55, NDC80, FEN1, DTL, RAD54L, KIAA0101, MCM6, ESPL1, HMMR, FOXM1, PRC1, SMC4, NEK2, EZH2, HELLS, WDHD1, ASF1B, CKS2, UHRF1, POLE2, ANLN, KIF14, UBE2T, BLM, ZWILCH, FAM83D. The 50 hub genes were showed as a heat map in Figure 3B, and the details about the 50 hub genes were showed in Supplementary Table S2.

Figure 3 PPI network analysis of the DEGs. (A) The protein–protein network construction of the 233 DEGs. The yellow nodes denote the 50 hub genes, the pink nodes denote the remaining 183 DEGs. Triangles represent upregulated genes, and V-shape represents downregulated genes. (B) The heatmap of the 50 hub genes. Red represents upregulated genes.

Survival Analysis

Kaplan–Meier plotter was employed to analyze the prognostic values of the 50 hub genes in thyroid carcinoma. Consistent with the survival analysis, the five genes (ANLN, CENPF, KIF2C, TPX2, and NDC80) were uploaded to the database and Kaplan–Meier curves were plotted. As shown in Figure 4, the red line denoted thyroid carcinoma patients with high expression of the hub gene, and the black line denoted thyroid carcinoma patients with low expression of the hub gene. The low gene expression group was the control group. Compared with the low gene expression group, the high expressions of ANLN (HR=2.99; 95% CI: 1.04–8.61), CENPF (HR=3.33; 95% CI: 1.16–9.85), KIF2C (HR=2.96; 95% CI: 1.1–7.97), TPX2 (HR=3.27; 95% CI: 1.14–9.42), and NDC80 (HR=2.59; 95% CI: 0.97–6.92) were associated with the poor survival of ATC patients (P<0.05). While the high expressions of other hub genes were not related to the overall survival of thyroid carcinoma patients. So the five genes were chosen from the 50 hub genes.

Figure 4 The prognostic curve of five genes. The prognostic significance of the hub genes in patients with TC, according to the Kaplan–Meier plotter database. (A) ANLN; (B) CENPF; (C) KIF2C; (D) TPX2; (E) NDC80. The red lines denote patients with high gene expression, and black lines denote patients with low gene expression.

In order to strengthen the explanation of the possible connection between the five genes, we used the String database to analyze the five genes again (Supplementary Figure S3). As shown in the Supplementary Figure S3A, the protein–protein interaction (PPI) network of the five genes consisted of the 10 edges with an average degree of 4. The top 10 significantly enriched terms in biological processes (BP), all terms in cell component (CC) and KEGG correspond to the criterion (P<0.05) were showed in Figure S3 B, C and D.

The Transcriptional Expression Level of the Five Genes

The five genes (ANLN, CENPF, KIF2C, TPX2, and NDC80) were significantly upregulated in ATC compared to that of normal thyroid tissues in microarray datasets GSE9115, GSE29265, GSE33630, GSE53072, and GSE65144 (Figure 5).

Figure 5 The mRNA expression of the five genes in five microarray datasets. The mRNA expressions of ANLN, CENPF, KIF2C, TPX2 and NDC80 in five microarray datasets (GSE9115, GSE29265, GSE33630, GSE53072, and GSE65144) respectively. The red dots represent a high gene expression, and blue dots represent a low gene expression. ** P < 0.01, *** P < 0.001 vs nontumor.

To further verify the results of bioinformatics analysis, the mRNA levels of the five genes were determined in normal thyroid follicular epithelial Nthy-ori 3–1 cells, anaplastic differentiated thyroid cancer ARO cells and CAL-62 cells with RT-qPCR. Each RT-qPCR was done at least three independent measurements. As shown in Figure 6, compared with the Nthy-ori 3–1 cells group, the expression of ANLN was decreased by 23.1% (P<0.01) and 36.9% (P<0.01) in ARO cells and CAL-62 cells, respectively. The expression of CENPF was upregulated by 39.6% (P<0.001) and 42.4% (P<0.01) in ARO cells and CAL-62 cells, respectively. The expression of KIF2C was upregulated by 1.3 times (P<0.05) and 1.5 times (P<0.05) in ARO cells and CAL-62 cells, respectively. The expression of TPX2 was decreased by 0.5% (P>0.05) and upregulated by 81.9% (P<0.05) in ARO cells and CAL-62 cells. The expression of NDC80 was decreased by 23.18.1% (P<0.001) and 39.5% (P<0.001) in ARO cells and CAL-62 cells, respectively. The results showed that KIF2C and CENPF were significantly upregulated in ARO cells and CAL-62 cells compared to that of Nthy-ori 3–1 cells, while TPX2 was upregulated only in CAL-62 cells. Inconsistent with the results in the GEO database, ANLN and NDC80 were obviously decreased in ARO cells and CAL-62 cells (Figure 6).

Figure 6 The mRNA expression of the five genes in cell lines. The mRNA expressions of the five genes were detected in normal thyroid follicular epithelial Nthy-ori 3–1 cells, anaplastic differentiated thyroid cancer ARO cells and CAL-62 cells by quantitative real-time PCR. (A) ANLN; (B) CENPF; (C) KIF2C; (D) TPX2; (E) NDC80. * P < 0.05, ** P < 0.01, *** P < 0.001 vs Nthy ori 3–1.

Discussion

Anaplastic thyroid cancer (ATC) is an undifferentiated, aggressive malignancy with a tendency of early distant metastases. It is challenging to diagnose ATC in the early stage, and once diagnosed, it is mostly in the middle or late stages, with short survival periods and high mortality.²⁸ At present, the treatment method of ATC is mainly through traditional surgery, but it has not achieved sound therapeutic effects.^29,30 With the rapid development of molecular biology, especially in the technical breakthrough of tumor molecular medicine, we have a better understanding of the occurrence and development of tumors. It is generally believed that ATC is caused by multiple genes, multiple pathways, various phases, and numerous factors. In recent years, the development of high-throughput chip technology has dramatically helped the screening and research of differential genes, which are expected to become one of the development directions of ATC therapy.³¹ Some potential therapeutic targets, such as PBF, PTTG, BCL2, and hsa-miR-181a-5p, have been found in thyroid carcinoma by using chip technology.^32,33 The aim of this study is to identify DEGs with an analysis of comprehensive bioinformatics and establish a foundation for future research of ATC. In this study, the expression profile data of ATC in the GEO database were excavated. The results showed that high levels of ANLN, CENPF, KIF2C, TPX2, and NDC80 indicated poor prognosis in ATC, suggesting that intervention of these genes might be a promising strategy for ATC treatment. Liu’s meta-analysis of microarray datasets identified TRIP13, DLGAP5, HJURP, CDKN3, NEK2, KIF15, TTK, KIF2C, AURKA and TPX2 as cell cycle-related key genes potentially contributing to anaplastic thyroid carcinoma,³⁴ which were in consistence with some of our results (KIF2C, TPX2).

Anillin (encoded by ANLN), an actin‐binding protein, is mainly located in the cytoskeleton and nucleus. It was first identified in Drosophila and plays a critical role in cytokinesis. In recent years, ANLN has been involved in the occurrence and development of various malignant tumors such as colon cancer, breast cancer,^35–37 pancreatic cancer,³⁸ which can promote the growth, migration, and infiltration of tumor cells.

Centromere Protein F (CENPF), a cell cycle-associated nuclear antigen, involves in the regulation of DNA synthesis and cell cycle progression. It is identified as a marker of cell proliferation and metastasis in several human carcinomas.^39–41 A large number of studies have shown that the expression of CENPF is abnormally elevated in various malignant tumors such as hepatocellular carcinomas,⁴⁰ prostate cancer,^41,42 breast cancer³⁹ and so on, which is related to the prognosis of patients.

Kinesin family member 2C (KIF2C), a mitotic centromere associated kinesin (MCAK), mainly localizes to the nucleoplasm, regulating mitosis and cell cycle. Some studies have shown that KIF2C is highly expressed in breast cancer, colon cancer, and other tumors, and is not beneficial to the prognosis of patients.⁴³ Zaganjor et al showed that KIF2C was highly expressed in non-small cell lung cancer and enhanced the migration ability of tumor cells.⁴⁴

Nuclear division cycle 80 (NDC80), a nuclear division cycle protein, mainly localizes to the nucleoplasm and centrosome, which is highly expressed in the cells with vigorous division, such as tumor cells. NDC80 plays an essential role in normal cell mitosis, maintenance of chromosomal stability, and tumor development. It has been found that overexpression of NDC80 could cause the mitotic checkpoint to be continuously hyperactivated and induce tumor formation in vivo. In osteosarcoma, the expression level of NDC80 was related to the TNM stage and distant metastasis, and NDC80 is an independent prognostic indicator.⁴⁵ Ju et al showed that NDC80 induced hepatocarcinogenesis mainly through its value-added and anti-apoptosis effects.

Targeting protein for Xenopus kinesin-like protein 2 (TPX2), mainly localizes to the nucleoplasm, is a newly discovered protein regulates cell cycle and associates with nuclear proliferation.⁴⁶ Recent studies have shown that TPX2 is related to tumorigenesis, its expression level in malignant tumors increases, and it can induce tumor differentiation and growth. The expression level of TPX2 is associated with metastasis and recurrence of multiple tumors, such as breast cancer,⁴⁷ bladder cancer,⁴⁸ hepatocellular carcinoma,⁴⁹ and so on. Blocking the expression of TPX2 may be a potential target for tumor therapy.

As shown in Supplementary Figure S3, the results indicated that the five genes play important roles in regulations of cell division, mitotic cell cycle process, metaphase plate congression, mitotic cell cycle phase transition, chromosome segregation, etc. In the CC ontology, the five genes were mainly enriched in condensed chromosome kinetochore, condensed chromosome outer kinetochore, cytoskeletal part, microtubule cytoskeleton, microtubule organizing center. In addition, the five genes were mainly enriched in seven KEGG pathways (p<0.05), which included Amplification of signal from unattached kinetochores via a MAD2 inhibitory signal, Separation of Sister Chromatids, Resolution of Sister Chromatid Cohesion, RHO GTPases Activate Formins, Cell Cycle, and G2/M Transition. Among them, Cell Cycle pathways showed the largest number of genes count: CENPF, KIF2C, NDC80, and TPX2.

The above research indicated that the five genes were mainly involved in the regulation of the cell cycle, so they may affect the reoccurrence and metastasis of ATC by regulating the cell cycle. In this study, the RT-qPCR was used to further verify the expression of these five genes in ATC cells. Consistent with the survival analysis, CENPF and KIF2C have a high mRNA level in ARO cells and CAL-62 cells, which may affect the progress of ATC by regulating the cell cycle. The mRNA levels of TPX2 were upregulated significantly only in CAL-62 cells and not in ARO cells compared to that of Nthy-ori 3–1 cells, suggesting that TPX2 overexpression may only affect the progress of CAL-62 cells by regulating cell cycle. However, ANLN and NDC80 were obviously decreased in ARO cells and CAL-62 cells. The mRNA expressions of ANLN and NDC80 are inconsistent with the results in the GEO database. Gene locates in the nucleus of the cell. So, the genes in cells and their tissues are the same. As for the inconsistent results, the possible reasons are as follows: (1) The five genes were screened from 42 ATC tissues of different patients at different stages of the disease (it covers more gene characteristics than two cell lines), while the two cell lines were immortalized from two ATC patients. The two cell lines could only show some gene characteristics of the disease. (2) The results in the GEO database used the RNA-seq method, but we used RT-qPCR method to verify the final five genes. RNA-seq is used for large-scale screening, reflecting the overall gene expression change trend of the sample, but it cannot be guaranteed that the changing trend of each gene is consistent with RT-qPCR.⁵⁰ RNA-seq and RT-qPCR are two different experimental platforms themselves, and the two cannot correspond to each other.^50,51 (3) GEO databases are tissue samples, our experimental samples are cells. Pan’s study indicated that the hub genes showed inconsistent trends in ATC tissues, which may be due to the heterogeneity between cell line data and tumor tissue data.⁵² So, only three of the five genes screened were validated in cell lines. The inconsistent results should be understandable. We will verify the mechanisms of CENPF, KIF2C, and TPX2 in the reoccurrence and metastasis of ATC in ARO cells and CAL-62 cells.

Researches in thyroid cancer have been hot in recent years. There are many articles published every year on thyroid cancer. In order to keep step with the new developments and findings, knowledge needs to be updated every two or three years, even in the same area of research. Although DEGs have been studied in ATC, different from previous studies, five microarray datasets were used for the analysis of DEGs in this study, whereas Hu’s only used three microarray datasets.⁵³ In addition, the analysis process was totally different from the previous studies: (1) Our study used five microarray datasets, whereas Hu’s only used three microarray datasets. Our five microarray datasets included the Hu’s three and two more datasets. So, our study analyzed seventy-six normal thyroid tissues and forty-two anaplastic thyroid carcinoma tissues while Hu’s study analyzed twenty-one normal thyroid tissues and twenty-two anaplastic thyroid carcinoma tissues. The more samples, the less bias.⁵³ (2) We had the volcano diagram (Figure 1) of gene expression profiles in ATC and normal samples. It was much easy for the readers to find the upregulated genes (presented in red) and downregulated genes (presented in blue), while Hu’s and Pan’s study did not have.^52,53 (3) We have used three generally recognized screening methods of the hub genes, including Molecular Complex Detection plugin (MCODE), Maximal Clique Centrality (MCC) and Degree Centrality (DC) in PPI network, and selected the overlapping genes of the three methods to make the final results most valuable. The overlapping genes were then conducted survival analysis using the Kaplan–Meier plotter to get the most important five genes that closely related to the survival of ATC patients. However, many studies only used MCODE to screen hub genes in PPI network.^16,24,54 Cao’s study used MCODE and MCC methods to screen hub genes in PPI network.¹⁶ Zhang’s study identified the height nodes as hub genes.⁵⁴ The top 10 DEGs with highest node degree were regarded as hub genes of ATC in Pan’s study.⁵² Hu’s study did not do the screening with any of those screening methods. They simply get the top 10 hub genes by PPI network.⁵³ In addition, preliminary verification in vitro was done and three of the five important genes were identified, while Hu’s study did not verify their results of Bioinformatics analysis with experiments.⁵³ In summary, our study analyzed more microarray datasets, used more reasonable and rigorous screening methods to screen the hub genes, and has verified the final five genes through cell experiments. The results were more convincing in our study compared with previous studies.

In summary, the DEGs related to the reoccurrence and metastasis of ATC were analyzed by gene chip combined with integrated bioinformatics analysis. Our results showed that CENPF, KIF2C, and TPX2 have high mRNA levels in ATC cells, which were significantly associated with thyroid carcinoma prognosis. And the results showed that KIF2C has the highest mRNA expression in ARO cells and CAL-62 cells among the three genes. In addition, the bioinformatics results showed that CENPF, KIF2C, and TPX2 might affect the progress of ATC by regulating cell cycle. Therefore, the druggability of the biomarkers CENPF, KIF2C, and TPX2 are expected to treat ATC patients by regulating cell cycle, and the druggability of the biomarkers KIF2C may have the best results, while their preclinical and clinical values require further validation. Our study results indicated that the three genes might be novel biomarkers for the treatment of ATC, and further experiments are required to verify the effects and specific mechanisms of these genes on ATC.