Identification of Key Functional Gene Signatures Indicative of Dedifferentiation in Papillary Thyroid Cancer

Background: Differentiated thyroid cancer (DTC) is the most common type of thyroid cancer. Many of them can relapse to dedifferentiated thyroid cancer (DDTC) and exhibit different gene expression profiles. The underlying mechanism of dedifferentiation and the involved genes or pathways remained to be investigated.

Methods: A discovery cohort obtained from patients who received surgical resection in the Fudan University Shanghai Cancer Center (FUSCC) and two validation cohorts derived from Gene Expression Omnibus (GEO) database were used to screen out differentially expressed genes in the dedifferentiation process. Weighted gene co-expression network analysis (WGCNA) was constructed to identify modules highly related to differentiation. Gene Set Enrichment Analysis (GSEA) was used to identify pathways related to differentiation, and all differentially expressed genes were grouped by function based on the GSEA and literature reviewing data. Least absolute shrinkage and selection operator (LASSO) regression analysis was used to control the number of variables in each group. Next, we used logistic regression to build a gene signature in each group to indicate differentiation status, and we computed receiver operating characteristic (ROC) curve to evaluate the indicative performance of each signature.

Results: A total of 307 upregulated and 313 downregulated genes in poorly differentiated thyroid cancer (PDTC) compared with papillary thyroid cancer (PTC) and normal thyroid (NT) were screened out in FUSCC cohort and validated in two GEO cohorts. WGCNA of 620 differential genes yielded the seven core genes with the highest correlation with thyroid differentiation score (TDS). Furthermore, 395 genes significantly correlated with TDS in univariate logistic regression analysis were divided into 11 groups. The areas under the ROC curve (AUCs) of the gene signature of group transcription and epigenetic modification, signal and substance transport, extracellular matrix (ECM), and metabolism in the training set [The Cancer Genome Atlas (TCGA) cohort] and validation set (combined GEO cohort) were both >0.75. The gene signature based on group transcription and epigenetic modification, cilia formation and movement, and proliferation can reflect the patient's disease recurrence state.

Conclusion: The dedifferentiation of DTC is affected by a variety of mechanisms including many genes. The gene signature of group transcription and epigenetic modification, signal and substance transport, ECM, and metabolism can be used as biomarkers for DDTC.

Introduction

Papillary thyroid cancer (PTC) is the most common type of thyroid cancer, and the majority of PTCs exhibit a relatively good prognosis (1, 2). However, it has been observed recently that some PTCs may dedifferentiate in some situations. When PTC appears to be dedifferentiated, its prognosis becomes very poor, and conventional surgical treatment cannot achieve good therapeutic effects. Such patients often experience relapse or metastasis in a short period of time (3–5). At present, the treatment methods for poorly differentiated thyroid cancer (PDTC) and anaplastic thyroid cancer (ATC) are limited, and their progression mechanism is still unclear.

It has been shown that many ATC and PDTC result from dedifferentiation of DTC (6). In addition, some genetic abnormalities such as TERT and TP53 mutation may play an important role (7, 8). A considerable number of studies have shown that occurrence and development of PDTC and ATC are closely related to immune microenvironment and epigenetic changes (9–12). Our previous study also revealed that some genes may have a significant impact on the initiation and progression of dedifferentiated thyroid cancer (DDTC) through metabolism-related pathways (13). However, considering that dedifferentiation of DTC is accompanied by a great increase in the degree of malignancy, it is likely that dedifferentiation must involve more than one mechanism. The mechanisms of how genes affect DTC dedifferentiation remain to be studied.

This study was oriented toward mining out differentially expressed genes among PDTC, PTC, and normal thyroid (NT) at the level of transcriptome and then classifying them into different groups based on their biological functions to explore possible dedifferentiation-related processes. We expect that our findings could provide a plausible basis for further study of PDTC, thereby helping to indicate prognosis and development of PTC and exploring the possibility of reversing dedifferentiation or re-differentiation.

Methods

Sample Collection

Six NT, five PTC, and five PDTC specimens were obtained from eight patients who underwent surgical management in the Fudan University Shanghai Cancer Center (FUSCC) (Supplementary Table 1). The information of the eight patients and 16 samples was described in our previous study (13). These 16 samples were included in a discovery cohort and used for high-throughput RNA sequencing (RNA-seq) to identify differentially expressed genes. Written informed consent was obtained from each patient before his/her specimens were used in this study, and the study was approved by the Medical Ethics Committee of the FUSCC. All procedures performed in this study were in accordance with the Declaration of Helsinki.

RNA-Seq Analysis

Total RNA was extracted from all samples using TRIzol reagent (Life Technologies, Carlsbad, CA). We used RiboMinus eukaryote kit (Qiagen, Valencia, CA) to remove ribosomal RNA of total RNA (~3 mg) before RNA-seq libraries construction. Strand-specific RNA-Seq libraries were prepared using the Illumina workflow (New England BioLabs, Beverly, MA). Next, the samples were fragmented, reverse-transcribed, and ligated to Illumina adaptors. We purified the ligated cDNA products to remove second-strand cDNA. After 13–15 cycles of amplification, libraries were controlled for quality and quantified using with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and sequenced by a HiSeq 2000 sequencing system (Illumina, San Diego, CA) on a 100-bp paired-end run. Clustal Omega was used for sequence alignment. Human genome version GRCh38.100 was used throughout. RNA-seq data were normalized by fragments per kilobase per million (FPKM). Significant differences were determined by Limma package (version 3.11; https://bioconductor.org/packages/limma/).

The Cancer Genome Atlas and Gene Expression Omnibus Cohorts

The Cancer Genome Atlas (TCGA) mRNA expression data (FPKM) was sourced from the University of California Santa Cruz (UCSC) Cancer Browser (14). The raw data of Gene Expression Omnibus (GEO) combined GEO cohort including the GSE29265 cohort (20 NTs, 20 PTCs, and 9 ATCs) (15), GSE33630 (16), GSE53157 (17), GSE65144 (18), and GSE76039 (19) were obtained from the GEO database (20, 21) and were background adjusted and normalized. The ComBat method was performed to remove batch effects by the R package “sva.” According to the annotation file, probes were matched with gene symbol, and probes that were not matched to gene symbol were deleted. When more than one probe matched the same gene symbol, average value was calculated as the final expression value. We performed a correlation analysis on expression matrix of PDTC and ATC samples in the combined GEO cohort, and we found that PDTC and ATC are highly correlated (Supplementary Figure 1). Therefore, we used the GSE29265 cohort and GSE33630 cohort as verification cohorts to verify the effectiveness of gene signatures. All patients were staged by the 8th edition of the TNM staging system published by the American Joint Committee on Cancer. Only samples for which all clinical data and thyroid differentiation score (TDS) could be obtained were included in the analysis.

Gene Function Annotation and Interaction Analysis

Gene Set Enrichment Analysis (GSEA) was performed using GSEA software. C5 collection [Gene Ontology (GO)] was utilized to identify GO terms that were differentially regulated in different comparisons. GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed by DAVID online function annotation tool (22, 23) to classify differentially expressed genes into functional categories. The enriched group was ranked by p-value (p < 0.05) and false discovery rate (FDR < 0.25). The GeneMANIA prediction server was used for interaction analysis between genes (24).

Weighted Gene Co-expression Network Analysis

Weighted gene co-expression network analysis (WGCNA) was performed by the R WGCNA package (25) (v1.66) to identify TDS-related modules and their gene members. The modules were identified by Dynamic Hybrid Tree Cut algorithm. We chose the yellow module with the highest correlation coefficient to screen out hub genes.

Statistical Analysis

A comparison of categorical variables was performed using chi-square test and Fisher's exact test. Descriptive statistics are presented in the tables. TDS was proposed by TCGA project. It summarizes the expression of 16 genes (DIO1, DIO2, DUOX1, DUOX2, FOXE1, GLIS3, NKXX2-1, PAX8, SLC26A4, SLC5A5, SLC5A8, TG, THRA, THRB, TPO, and TSHR) related to thyroid metabolism and function, which can be an index to determine the degree of thyroid-specific differentiation. We have eliminated these genes in subsequent calculations to avoid bias (26). According to median of TDS, TCGA dataset was divided into two groups: highly differentiated and poorly differentiated. Least absolute shrinkage and selection operator (LASSO) regression was performed on each function group to reduce the number of variables. Multivariate logistic regression was used to determine genes to be finally included in the signature; expression levels of the differential genes as signatures in each patient were integrated into a risk score fitted by logistic regression (Table 1). Non-parametric receiver operating characteristic (ROC) analysis was performed for each signature, and we calculated area under the ROC curve (AUC) to test its indicative power for dedifferentiation. The Kaplan–Meier method was used to construct the disease-free survival (DFS) curve. Two-side p < 0.05 was considered statistically significant. Statistical analysis and visualization were carried out through SPSS (v25.0; IBM Corporation, Armonk, NY) and R (v3.6.3; R Foundation for Statistical Computing, Vienna, Austria) software.

TABLE 1

Table 1. Function based differential gene grouping.

Results

Identification of Genes Related to Differentiation Status in Dedifferentiated Thyroid Cancer

To identify genes that may relate to dedifferentiation of DTC, we first analyzed changes in high-throughput transcriptome expression profile of five PDTCs, five PTCs, and six NTs from the FUSCC cohort. Then, 1,465 deregulated genes (690 upregulated and 775 downregulated, fold change ≥ 2, p < 0.05) were found among three groups (Figure 1A). Then we found out the corresponding expression values of these genes in two GEO datasets (GSE29265 and GSE33630) and confirmed their differential expression (307 upregulated genes and 313 downregulated genes, fold change ≥ 1, p < 0.05), as shown in Venn diagrams and heat maps, which finally validated our target gene set for an in-depth study (Figures 1B–F).

FIGURE 1

Figure 1. Identification of deregulation of genes in DDTC. (A) Volcano maps: left map, 1,728 downregulated genes (green dots) and 988 upregulated genes (red dots) in PDTCs compared with NT tissues; right map, 2,630 downregulated genes (green dots) and 1,246 upregulated genes (red dots) in PDTCs compared with PTCs; fold change (FC) ≥ 2, p < 0.05. (B) Compared with PTC and NTs, overlap of abnormal expressed genes in PDTC. (C) The Venn diagram showed that 1,465 genes abnormally expressed in the FUSCC cohort were validated into GSE29265 and GSE33630, and 313 downregulated genes and 307 upregulated genes were screened out. (D–F) The heat map shows the expression data of 620 dysregulated genes in the FUSCC cohort, GSE29265 cohort, and GSE33630 cohort. DDTC, dedifferentiated thyroid cancer; PDTC, poorly differentiated thyroid cancer; NT, normal thyroid; PTC, papillary thyroid cancer; FUSCC, Fudan University Shanghai Cancer Center.

Weighted Gene Co-expression Network Analysis Obtained a Co-expression Module Containing Seven Hub Genes

The WGCNA of 620 valid genes was carried out. Seven gene co-expression modules were detected. WGCNA assigned colors to name each module, and the yellow module showed the highest correlation with TDS (Figures 2A–D). Then we screened genes by direct correlation between genes and specified traits, module identity, and weighted correlation; METTL7A, KCNQ1, ALDH9A1, C16orf46, PLAUR, BCL2, and TPMT were defined as the hub genes. The expression levels of these seven hub genes in each patient were fitted to a dedifferentiation risk score through logistic regression. The AUC value of this risk score in TCGA cohort reached 0.91(Figure 2F) and 0.73 in GEO combined cohort (Supplementary Figure 2). However, after the patients were divided into two groups based on the risk score, the difference in DFS between two groups was not significant (Figure 2E).

FIGURE 2

Figure 2. Gene modules identified by weighted gene co-expression network analysis (WGCNA) and indicative ability of hub gene signature constructed by logistic regression for PTC differentiation and disease free survival rate. (A) Cluster the gene dendrograms according to topological overlap and assign different colors to each module. Finally, seven co-expression modules were constructed. The size of these modules varies according to the number of genes they contain. (B) Module–TDS associations. Each row corresponds to a module eigengene; each column corresponds to TDS. Each cell contains the corresponding correlation and p-value. (C) Using heat map plots to visualize gene networks. The heat map depicts the topological overlap matrix (TOM) between all genes in the analysis. (D) Heat map of module adjacency. Red indicates high adjacency (positive correlation), and blue indicates low adjacency (negative correlation). (E) The receiver operating characteristic (ROC) shows indicative ability of hub gene signature for TDS. (F) Hub gene signature fails to distinguish DFS well. PTC, papillary thyroid cancer; TDS, thyroid differentiation score; DFS, disease-free survival.

Functional Classification of Differentially Expressed Genes

The seven hub genes obtained by WGCNA showed no good results in functional enrichment and thus cannot fully reflect the impact of abnormally expressed genes on differentiation. Therefore, we screened 620 genes through univariate logistic regression analysis and found 396 genes that significantly correlated with TDS. Next, we tried to group these genes into functional groups based on the results of GSEA. We found that the results of GSEA were mostly concentrated in proliferation, immunity, metabolism, and other related pathways (Supplementary Tables 2, 3). The most significant phenotype of tumor cells after dedifferentiation was the proliferation ability (27, 28); as a result, some phenotypes of interest may have been masked by proliferation during GSEA.

Therefore, we combined the results of GSEA and data from the literature to divide all deregulated genes into 11 groups including cell differentiation and rhythm, substance and signal transduction, transcription and epigenetic modification, protein modification, cilia formation and movement, extracellular matrix (ECM), apoptosis, proliferation, invasion, metabolism, and immunity (Table 2). In addition, there were 44 genes that cannot be clearly functionally classified (Table 1; Supplementary Table 4).

TABLE 2

Table 2. Risk score for each functional group.

Gene Signatures of Signal and Substance Transport, Transcription and Epigenetic Modification, Extracellular Matrix, and Metabolism Group Can Indicate Dedifferentiation of Papillary Thyroid Cancer

Among the 396 genes obtained in multivariate regression analysis, many genes have been extensively studied. For example, EZH2, a methyltransferase, is closely related to tumor metastasis and proliferation (29, 30). TNC, an extracellular matrix protein, plays an important role in tumor cell invasion and metastasis (31, 32). PRDM16, which also plays an important role in protein modification, is also closely related to fat metabolism and tumor growth (33, 34). In addition, there are many important tumor-related genes. We have further screened them in each group to construct a functionally related gene signature.

We first used LASSO regression to screen each group of genes (Supplementary Figure 3) and then performed multivariate logistic regression analysis adjusted by T stage, lymph node metastasis (LNM), and BRAF^V600E on each group of genes selected to determine genes to be incorporated into the signature, and we evaluated the indicative ability of each signature by calculating the AUC value in the training set (TCGA cohort) and validation set (GEO combined cohort). We found that the AUC value in the group signal and substance transport, transcription and epigenetic modification, ECM, and metabolism was >0.75 in both the training set and validation set (Figure 3). It was worth noting that in two datasets, the AUC values of the invasion group reached 0.919297 and 0.746622, and the cilia formation and movement group also reached 0.805777 and 0.709767 (Supplementary Figure 4).

FIGURE 3

Figure 3. The AUC of four signatures of metabolism, signal and substance transport, transcription and epigenetic modification, and extracellular matrix. (A–D) The AUC of the group metabolism, signal and substance transport, transcription and epigenetic modification, and extracellular matrix was >0.75 in both the training set and validation set. AUC, area under the receiver operating characteristic curve.

Correlations of the Gene Signatures With Clinical Features of Papillary Thyroid Cancer

In further analyses, we investigated the clinicopathological significance of the gene signatures associated with dedifferentiation to reveal their further research priority and translational potential.

Based on the gene signature of each group, we divided all TCGA patients into high-risk and low-risk groups, and we analyzed their survival differences. We found that although some risk score of gene signatures can reflect degree of differentiation well, they cannot accurately reflect the patient's tumor recurrence status. Among all gene signatures, there were significant differences in DFS between the high-risk and low-risk groups including transcription and epigenetic modification, cilia formation and movement, and proliferation (Figure 4; Supplementary Figure 5). Next, we further analyzed the relationships between high or low risk and clinical parameters in TCGA cohort (Table 3), and we found significant differences in BRAF^V600E mutation, TNM stage, T stage, and TDS. In addition, in the cilia formation group, the risk score and LNM also had a significant correlation (Figure 5).

FIGURE 4

Figure 4. Gene signatures were tested for DFS analyses in patients with thyroid cancer. (A–C) There were significant differences in DFS between the high-risk group and low-risk group in cilia formation and movement, transcription and epigenetic modification, and proliferation groups. DFS, disease-free survival.

TABLE 3

Table 3. Clinicopathologics of PTC patients in TCGA cohort.

FIGURE 5

Figure 5. Correlation of risk score with clinical parameters in TCGA cohort. (A) The risk score in cilia formation and movement group was significantly correlated with TDS, BRAF^V600E mutation, LNM, T stage, and extrathyroidal extension (ETE). (B) The risk score in the transcription and epigenetic modification group was significantly correlated with TDS, BRAF^V600E mutation, LNM, T stage, and ETE. (C) The risk score in the proliferation group was significantly correlated with TDS, BRAF^V600E mutation, T stage, and ETE. *p < 0.05, **p < 0.01, and ***p < 0.001. TCGA, The Cancer Genome Atlas; TDS, thyroid differentiation score; LNM, lymph node metastasis.

Discussion

Undifferentiated thyroid cancer is relatively rare, but it has a very high degree of malignancy, which brings great difficulties in exploring its pathogenesis. In recent years, many studies have explored undifferentiated thyroid cancer at pre-transcription, transcription, and translation levels, revealing that a considerable number of undifferentiated cancers develop from dedifferentiation of DTC. Studies have shown that dedifferentiation of colorectal cancer is closely related to TGF-β, Wnt, and Hedgehog signaling pathways (35), while dedifferentiation of glioblastoma is closely related to hypoxia (36). Our previous studies have also shown that metabolic pathways play an important role in dedifferentiation of thyroid cancer (13). Most of these studies focused on a certain pathway or mechanism; however, according to previous experience, the occurrence and development of malignant tumors, including dedifferentiation, involve multiple mechanisms and a large number of abnormal gene expressions. Moreover, research on the mechanism of dedifferentiation of DTC is still not comprehensive. In order to explore the potential mechanism of dedifferentiation of DTC from a broader perspective, we retrospectively obtained 16 samples of eight patients with DDTC who had undergone surgical treatment in our institution (FUSCC cohort). We used intermittent sampling to obtain tissue samples with a gradient of differentiation. Through the joint analysis of 16 samples of high-throughput sequencing data, TCGA cohort, and GEO chip data, we obtained 620 differentially expressed genes related to dedifferentiation of PTC, and all of them were subjected to further analysis.

First, we performed WGCNA on all 620 genes and divided them into seven modules according to their expression. Among them, the yellow module demonstrated the strongest correlation with TDS. After screening genes by direct correlation between genes and specified traits, module identity, and weighted correlation, seven genes including METTL7A, KCNQ1, ALDH9A1, C16orf46, PLAUR, BCL2, and TMPT have been certified as hub genes. Through analysis of GeneMANIA prediction server, we found that these seven hub genes and their interacting genes directly or indirectly interact and co-localize, but they were not functionally enriched together (Supplementary Figure 6). Thus, we believe that the results of WGCNA cannot fully explain the dedifferentiation process of DTC.

In order to reflect the differentiation status more comprehensively from multiple angles, we first performed GSEA on all genes and found that the pathways related to TDS were mostly enriched on proliferation and apoptosis, which is consistent with the biological performance of PDTC. Many studies also focused on pathways related to cell proliferation, invasion, and immunity (26, 37–39); but also many new tumor-related phenotypes, such as ECM, cilia formation and movement, and epigenetic modification, have attracted our interest (40–43). However, these phenotypes have been rarely studied in thyroid cancer, so there is a certain value for further research (9, 44, 45).

We evaluated the gene signature of each group by calculating AUC, and we verified it in the combined GEO cohort. Finally, we identified gene signatures in group signal and substance transport, transcription and epigenetic modification, ECM, invasion, and metabolism groups, which showed very good indicative performance. In the gene signature of signal and substance transport group, PEBP1, FZD5, KCNQ1, and TOM1L2 play a role in protein kinase binding together (Supplementary Figure 7A), while PEBP1 is also associated with autophagy and ferroptotic death (46). In the transcription and epigenetic modification group, SMAD6 and PRDM16 showed a high degree of consistency in function (Supplementary Figure 7B), as functions of these two genes are concentrated on TGF-β/SMAD (47, 48), a pathway that plays an important role in cell differentiation. SMAD6 is a component of this pathway; and its importance is self-evident, as PRDM16, an important methyltransferase and transcription factor, has also been reported as a repressor of the TGF-β/SMAD pathway (49). In the ECM group, TNC was in a central position, and its main interaction target was Integrin Alpha V (ITGAV), a protein expressed in thyroid tissues higher than in other tissues in the human body, suggesting that it specifically influences the differentiation of thyroid cancer through the phenotype of extracellular matrix (Supplementary Figure 7C). As a traditional malignant tumor phenotype, the gene signature of the invasion group consisted of NHSL1, OCLN, EPB41L4B, PLAUR, PIK3CB, EPB41L5, ID4, TPM4, and MARVELD2. TPM4 showed the opposite effects in colon cancer and lung cancer. In our analysis, its high expression corresponded to better differentiation, which is consistent with the research results in colon cancer (50, 51). In the metabolism group, CHPT1, PDE8B, HSD17B8, ADAMTSL1, ST3GAL1, SLC7A5, ACOT7, and PYGL were combined to make up a signature. PDE8B, HSD17B8, ST3GAL1, and ACOT7 also appeared in a signature constructed in our previous study (13).

Through the STRING database, we found that some genes have an interaction relationship, including intra-group and inter-group interactions (52) (Supplementary Figure 8). In the transcription and epigenetic modification group, there are gene fusion and co-expression between TBX3 and SMAD6 (53). OCLN and MARVELD2 are both invasion-related genes, and they share a protein homology (54). MARVELD2 also has a co-expression relationship with the signal and substance transport gene ILDR2 (54). KLHDC1 in the signal and substance transport group, and IQGAP3 in the proliferation group also has co-expression (55). These possible interactions indicate that the joint effect of these genes may have a more important impact on differentiation of DTC. Since they have not been studied in thyroid cancer, it is worthy of further investigation.

At the same time, we also paid due attention to the relationship between risk scores and other clinical parameters. We found that the cilia formation group could better reflect the condition of LNM than other groups. Other studies have shown that cilia formation is closely related to autophagy (56); therefore, high expression of cilia formation-related proteins (DNAH7, TAPT1, and BBS1) in the low-risk group may be due to exuberant autophagy activity (Supplementary Figure 7D). Additionally, the formation and movement of cilia are also very important to the movement of cells, which in turn affect the tumor cell migration and invasion (57). In recent years, studies have found that the sensitivity of ATC to chemotherapy drugs is closely related to autophagy (58–61); therefore, we believe that genes involved in cilia formation and movement are promising as targets for drug therapy and should be further researched.

Our research also had certain limitations. First, we selected traditional malignant tumor-related phenotypes and some phenotypes that we were interested in as groupings, so this method is not very objective. There may be some important phenotypes that we did not include in the groups. This has also led to 44 genes that could not be included in specific group. At the same time, it should be noted that a more detailed and comprehensive grouping may divide the genes too finely, resulting in some interactions that cannot be reflected in intergroup. Second, limited by the number of cases and difficulty of obtaining material, the distinction between PDTC and PTC in the sequenced samples is not fine enough. We plan to collect more samples and use more precise microdissection methods for the next step of research.

In conclusion, we analyzed the genes that may affect the differentiation of thyroid cancer from 11 different perspectives and constructed gene signatures, revealing the possibility of multiple mechanisms that lead to dedifferentiation of DTC.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: National Center for Biotechnology Information (NCBI) BioProject, https://www.ncbi.nlm.nih.gov/bioproject/, PRJNA702648.

Ethics Statement

This study was approved by the Medical Ethics Committee of the FUSCC. All procedures performed in our study were in accordance with the ethical standards of our institutional research committee and the Declaration of Helsinki.

Author Contributions

WX, CL, BM, and WW contributed to conception and design of the study. HJ, XW, and YY organized the data. WX, CL, and YucW performed the statistical analysis. All authors contributed to writing and review of the manuscript and approved the final submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (Grants 81572622 and 81772854 to QJ and 82072951 to YuW) and the Natural Science Foundation of Shanghai (Grant 19ZR1410900 to WW).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The UCSC Xena Platform for cancer genomics datasets, the GEO platform of the National Center for Biotechnology Information, and the Science and Technology Commission of Shanghai Municipality are gratefully acknowledged. We were also very grateful to Haihan Chen for his help in the current research.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.641851/full#supplementary-material

References

1. Aschebrook-Kilfoy B, Ward MH, Sabra MM, Devesa SS. Thyroid cancer incidence patterns in the United States by histologic type, 1992-2006. Thyroid. (2011) 21:125–34. doi: 10.1089/thy.2010.0021

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Fagin JA, Wells SA Jr. Biologic and clinical perspectives on thyroid cancer. N Engl J Med. (2016) 375:2307. doi: 10.1056/NEJMc1613118

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Zhang C, Fu X, Chen P, Bao X, Li F, Sun X, et al. Dedifferentiation derived cells exhibit phenotypic and functional characteristics of epidermal stem cells. J Cell Mol Med. (2010) 14:1135–45. doi: 10.1111/j.1582-4934.2009.00765.x

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ibrahimpasic T, Ghossein R, Shah JP, Ganly I. Poorly differentiated carcinoma of the thyroid gland: current status and future prospects. Thyroid. (2019) 29:311–21. doi: 10.1089/thy.2018.0509

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Molinaro E, Romei C, Biagini A, Sabini E, Agate L, Mazzeo S, et al. Anaplastic thyroid carcinoma: from clinicopathology to genetics and advanced therapies. Nat Rev Endocrinol. (2017) 13:644–60. doi: 10.1038/nrendo.2017.76

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Ma B, Xu W, Wei W, Wen D, Lu Z, Yang S, et al. Clinicopathological and survival outcomes of well-differentiated thyroid carcinoma undergoing dedifferentiation: a retrospective study from FUSCC. Int J Endocrinol. (2018) 2018:2383715. doi: 10.1155/2018/2383715

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Pozdeyev N, Gay LM, Sokol ES, Hartmaier R, Deaver KE, Davis S, et al. Genetic analysis of 779 advanced differentiated and anaplastic thyroid cancers. Clin Cancer Res. (2018) 24:3059–68. doi: 10.1158/1078-0432.CCR-18-0373

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Yoo SK, Song YS, Lee EK, Hwang J, Kim HH, Jung G, et al. Integrative analysis of genomic and transcriptomic characteristics associated with progression of aggressive thyroid cancer. Nat Commun. (2019) 10:2764. doi: 10.1038/s41467-019-10680-5

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Sasanakietkul T, Murtha TD, Javid M, Korah R, Carling T. Epigenetic modifications in poorly differentiated and anaplastic thyroid cancer. Mol Cell Endocrinol. (2018) 469:23–37. doi: 10.1016/j.mce.2017.05.022

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Xing M. Genetic alterations in the phosphatidylinositol-3 kinase/Akt pathway in thyroid cancer. Thyroid. (2010) 20:697–706. doi: 10.1089/thy.2010.1646

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Bastman JJ, Serracino HS, Zhu Y, Koenig MR, Mateescu V, Sams SB, et al. Tumor-infiltrating T cells and the PD-1 checkpoint pathway in advanced differentiated and anaplastic thyroid cancer. J Clin Endocrinol Metab. (2016) 101:2863–73. doi: 10.1210/jc.2015-4227

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Giannini R, Moretti S, Ugolini C, Macerola E, Menicali E, Nucci N, et al. Immune profiling of thyroid carcinomas suggests the existence of two major phenotypes: an ATC-like and a PDTC-like. J Clin Endocrinol Metab. (2019) 104:3557–75. doi: 10.1210/jc.2018-01167

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ma B, Jiang H, Wen D, Hu J, Han L, Liu W, et al. Transcriptome analyses identify a metabolic gene signature indicative of dedifferentiation of papillary thyroid cancer. J Clin Endocrinol Metab. (2019) 104:3713–25. doi: 10.1210/jc.2018-02686

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Goldman MJ, Craft B, Hastie M, Repecka K, McDade F, Kamath A, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. (2020) 38:675–8. doi: 10.1038/s41587-020-0546-8

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Tomas G, Tarabichi M, Gacquer D, Hebrant A, Dom G, Dumont JE, et al. A general method to derive robust organ-specific gene expression-based differentiation indices: application to thyroid cancer diagnostic. Oncogene. (2012) 31:4490–8. doi: 10.1038/onc.2011.626

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Dom G, Tarabichi M, Unger K, Thomas G, Oczko-Wojciechowska M, Bogdanova T, et al. A gene expression signature distinguishes normal tissues of sporadic and radiation-induced papillary thyroid carcinomas. Br J Cancer. (2012) 107:994–1000. doi: 10.1038/bjc.2012.302

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Pita JM, Banito A, Cavaco BM, Leite V. Gene expression profiling associated with the progression to poorly differentiated thyroid carcinomas. Br J Cancer. (2009) 101:1782–91. doi: 10.1038/sj.bjc.6605340

PubMed Abstract | CrossRef Full Text | Google Scholar

18. von Roemeling CA, Marlow LA, Pinkerton AB, Crist A, Miller J, Tun HW, et al. Aberrant lipid metabolism in anaplastic thyroid carcinoma reveals stearoyl CoA desaturase 1 as a novel therapeutic target. J Clin Endocrinol Metab. (2015) 100:E697–709. doi: 10.1210/jc.2014-2764

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Landa I, Ibrahimpasic T, Boucai L, Sinha R, Knauf JA, Shah RH, et al. Genomic and transcriptomic hallmarks of poorly differentiated and anaplastic thyroid cancers. J Clin Invest. (2016) 126:1052–66. doi: 10.1172/JCI85271

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res. (2013) 41(Database issue):D991–D5. doi: 10.1093/nar/gks1193

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Edgar R, Domrachev M, Lash AE. Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. (2002) 30:207–10. doi: 10.1093/nar/30.1.207

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. (2009) 4:44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Galperin MY, Cochrane GR. Nucleic acids research annual database issue and the NAR online molecular biology database collection in 2009. Nucleic Acids Res. (2009) 37(Database issue):D1–D4. doi: 10.1093/nar/gkn942

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Warde-Farley D, Donaldson SL, Comes O, Zuberi K, Badrawi R, Chao P, et al. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. (2010) 38(Web Server issue):W214–W20. doi: 10.1093/nar/gkq537

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell. (2014) 159:676–90. doi: 10.1016/j.cell.2014.09.050

CrossRef Full Text | Google Scholar

27. Hu YB, Yan C, Mu L, Mi YL, Zhao H, Hu H, et al. Exosomal Wnt-induced dedifferentiation of colorectal cancer cells contributes to chemotherapy resistance. Oncogene. (2019) 38:1951–65. doi: 10.1038/s41388-018-0557-9

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Parker TM, Henriques V, Beltran A, Nakshatri H, Gogna R. Cell competition and tumor heterogeneity. Semin Cancer Biol. (2020) 63:1–10. doi: 10.1016/j.semcancer.2019.09.003

CrossRef Full Text | Google Scholar

29. Duan R, Du W, Guo W. EZH2: a novel target for cancer treatment. J Hematol Oncol. (2020) 13:104. doi: 10.1186/s13045-020-00937-8

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Yuan H, Han Y, Wang X, Li N, Liu Q, Yin Y, et al. SETD2 restricts prostate cancer metastasis by integrating EZH2 and AMPK signaling pathways. Cancer Cell. (2020) 38:350–65.e7. doi: 10.1016/j.ccell.2020.05.022

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Oskarsson T, Acharyya S, Zhang XH, Vanharanta S, Tavazoie SF, Morris PG, et al. Breast cancer cells produce tenascin C as a metastatic niche component to colonize the lungs. Nat Med. (2011) 17:867–74. doi: 10.1038/nm.2379

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hoshino A, Kim HS, Bojmar L, Gyan KE, Cioffi M, Hernandez J, et al. Extracellular vesicle and particle biomarkers define multiple human cancers. Cell. (2020) 182:1044–61.e18. doi: 10.1016/j.cell.2020.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Kundu A, Nam H, Shelar S, Chandrashekar DS, Brinkley G, Karki S, et al. PRDM16 suppresses HIF-targeted gene expression in kidney cancer. J Exp Med. (2020) 217:e20191005. doi: 10.1084/jem.20191005

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wang W, Ishibashi J, Trefely S, Shao M, Cowan AJ, Sakers A, et al. A PRDM16-Driven metabolic signal from adipocytes regulates precursor cell fate. Cell Metab. (2019) 30:174–89.e5. doi: 10.1016/j.cmet.2019.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Oku Y, Shimoji T, Takifuji K, Hotta T, Yokoyama S, Matsuda K, et al. Identification of the molecular mechanisms for dedifferentiation at the invasion front of colorectal cancer by a gene expression analysis. Clin Cancer Res. (2008) 14:7215–22. doi: 10.1158/1078-0432.CCR-08-0370

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wang P, Yan Q, Liao B, Zhao L, Xiong S, Wang J, et al. The HIF1α/HIF2α-miR210-3p network regulates glioblastoma cell proliferation, dedifferentiation and chemoresistance through EGF under hypoxic conditions. Cell Death Dis. (2020) 11:992. doi: 10.1038/s41419-020-03150-0

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Caillou B, Talbot M, Weyemi U, Pioche-Durieu C, Al Ghuzlan A, Bidart JM, et al. Tumor-associated macrophages (TAMs) form an interconnected cellular supportive network in anaplastic thyroid carcinoma. PLoS One. (2011) 6:e22567. doi: 10.1371/journal.pone.0022567

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Chang CJ, Chao CH, Xia W, Yang JY, Xiong Y, Li CW, et al. p53 regulates epithelial-mesenchymal transition and stem cell properties through modulating miRNAs. Nat Cell Biol. (2011) 13:317–23. doi: 10.1038/ncb2173

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Garcia-Rendueles ME, Ricarte-Filho JC, Untch BR, Landa I, Knauf JA, Voza F, et al. NF2 loss promotes oncogenic RAS-induced thyroid cancers via YAP-dependent transactivation of RAS proteins and sensitizes them to MEK inhibition. Cancer Discov. (2015) 5:1178–93. doi: 10.1158/2159-8290.CD-15-0330

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Xie Q, Wu TP, Gimple RC, Li Z, Prager BC, Wu Q, et al. N(6)-methyladenine DNA modification in glioblastoma. Cell. (2018) 175:1228–43.e20. doi: 10.1016/j.cell.2018.10.006

CrossRef Full Text | Google Scholar

41. Liu H, Kiseleva AA, Golemis EA. Ciliary signalling in cancer. Nat Rev Cancer. (2018) 18:511–24. doi: 10.1038/s41568-018-0023-6

CrossRef Full Text | Google Scholar

42. Insua-Rodriguez J, Oskarsson T. The extracellular matrix in breast cancer. Adv Drug Deliv Rev. (2016) 97:41–55. doi: 10.1016/j.addr.2015.12.017

CrossRef Full Text | Google Scholar

43. Kalluri R. The biology and function of fibroblasts in cancer. Nat Rev Cancer. (2016) 16:582–98. doi: 10.1038/nrc.2016.73

CrossRef Full Text | Google Scholar

44. Lee J, Yi S, Won M, Song YS, Yi HS, Park YJ, et al. Loss-of-function of IFT88 determines metabolic phenotypes in thyroid cancer. Oncogene. (2018) 37:4455–74. doi: 10.1038/s41388-018-0211-6

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Nucera C, Lawler J, Parangi S. BRAF(V600E) and microenvironment in thyroid cancer: a functional link to drive cancer progression. Cancer Res. (2011) 71:2417–22. doi: 10.1158/0008-5472.CAN-10-3844

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Zhao J, Dar HH, Deng Y, St Croix CM, Li Z, Minami Y, et al. PEBP1 acts as a rheostat between prosurvival autophagy and ferroptotic death in asthmatic epithelial cells. Proc Natl Acad Sci U S A. (2020) 117:14376–85. doi: 10.1073/pnas.1921618117

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Jiao J, Zhang R, Li Z, Yin Y, Fang X, Ding X, et al. Nuclear Smad6 promotes gliomagenesis by negatively regulating PIAS3-mediated STAT3 inhibition. Nat Commun. (2018) 9:2504. doi: 10.1038/s41467-018-04936-9

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Su L, Lei X, Ma H, Feng C, Jiang J, Jiao J. PRDM16 orchestrates angiogenesis via neural differentiation in the developing brain. Cell Death Differ. (2020) 27:2313–29. doi: 10.1038/s41418-020-0504-5

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Yoshida M, Nosaka K, Yasunaga J, Nishikata I, Morishita K, Matsuoka M. Aberrant expression of the MEL1S gene identified in association with hypomethylation in adult T-cell leukemia cells. Blood. (2004) 103:2753–60. doi: 10.1182/blood-2003-07-2482

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Yang R, Zheng G, Ren D, Chen C, Zeng C, Lu W, et al. The clinical significance and biological function of tropomyosin 4 in colon cancer. Biomed Pharmacother. (2018) 101:1–7. doi: 10.1016/j.biopha.2018.01.166

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Zhao X, Jiang M, Wang Z. TPM4 promotes cell migration by modulating F-actin formation in lung cancer. Onco Targets Ther. (2019) 12:4055–63. doi: 10.2147/OTT.S198542

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. (2019) 47:D607–D13. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Prados B, Gomez-Apinaniz P, Papoutsi T, Luxan G, Zaffran S, Perez-Pomares JM, et al. Myocardial Bmp2 gain causes ectopic EMT and promotes cardiomyocyte proliferation and immaturity. Cell Death Dis. (2018) 9:399. doi: 10.1038/s41419-018-0442-z

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Suzuki H, Nagatake T, Nasu A, Lan H, Ikegami K, Setou M, et al. Impaired airway mucociliary function reduces antigen-specific IgA immune response to immunization with a claudin-4-targeting nasal vaccine in mice. Sci Rep. (2018) 8:2904. doi: 10.1038/s41598-018-21120-7

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Ghosh S, Zang S, Mitra PS, Ghimbovschi S, Hoffman EP, Dutta SK. Global gene expression and ingenuity biological functions analysis on PCBs 153 and 138 induced human PBMC in vitro reveals differential mode(s) of action in developing toxicities. Environ Int. (2011) 37:838–57. doi: 10.1016/j.envint.2011.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Pampliega O, Orhon I, Patel B, Sridhar S, Diaz-Carretero A, Beau I, et al. Functional interaction between autophagy and ciliogenesis. Nature. (2013) 502:194–200. doi: 10.1038/nature12639

CrossRef Full Text | Google Scholar

57. Mansini AP, Peixoto E, Jin S, Richard S, Gradilone SA. The chemosensory function of primary cilia regulates cholangiocyte migration, invasion, and tumor growth. Hepatology. (2019) 69:1582–98. doi: 10.1002/hep.30308

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Liu F, Zhang J, Qin L, Yang Z, Xiong J, Zhang Y, et al. Circular RNA EIF6 (Hsa_circ_0060060) sponges miR-144-3p to promote the cisplatin-resistance of human thyroid carcinoma cells by autophagy regulation. Aging (Albany NY). (2018) 10:3806–20. doi: 10.18632/aging.101674

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Feng H, Cheng X, Kuang J, Chen L, Yuen S, Shi M, et al. Apatinib-induced protective autophagy and apoptosis through the AKT-mTOR pathway in anaplastic thyroid cancer. Cell Death Dis. (2018) 9:1030. doi: 10.1038/s41419-018-1054-3

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Tesselaar MH, Crezee T, Schuurmans I, Gerrits D, Nagarajah J, Boerman OC, et al. Digitalislike compounds restore hNIS expression and iodide uptake capacity in anaplastic thyroid cancer. J Nucl Med. (2018) 59:780–6. doi: 10.2967/jnumed.117.200675

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Krug SM, Bojarski C, Fromm A, Lee IM, Dames P, Richter JF, et al. Tricellulin is regulated via interleukin-13-receptor α2, affects macromolecule uptake, and is decreased in ulcerative colitis. Mucosal Immunol. (2018) 11:345–56. doi: 10.1038/mi.2017.52

PubMed Abstract | CrossRef Full Text | Google Scholar