Characterizing dedifferentiation of thyroid cancer by integrated analysis

Abstract

Understanding of dedifferentiation, an indicator of poo prognosis for patients with thyroid cancer, has been hampered by imprecise and incomplete characterization of its heterogeneity and its attributes. Using single-cell RNA sequencing, we explored the landscape of thyroid cancer at single-cell resolution with 46,205 cells and delineated its dedifferentiation process and suppressive immune microenvironment. The developmental trajectory indicated that anaplastic thyroid cancer (ATC) cells were derived from a small subset of papillary thyroid cancer (PTC) cells. Moreover, a potential functional role of CREB3L1 on ATC development was revealed by integrated analyses of copy number alteration and transcriptional regulatory network. Multiple genes in differentiation-related pathways (e.g., EMT) were involved as the downstream targets of CREB3L1, increased expression of which can thus predict higher relapse risk of PTC. Collectively, our study provided insights into the heterogeneity and molecular evolution of thyroid cancer and highlighted the potential driver role of CREB3L1 in its dedifferentiation process.

Fig. 1. Single-cell profile of thyroid cancers.

(A) Overview of the experiment procedure. Ten scRNA-seq data were generated from normal, PTC, and ATC tissues using 10× Genomics protocol. We analyzed the transcriptome of 46,205 individual cells. In parallel, WES was performed. FFPE, formalin-fixed, paraffin-embedded; H&E, hematoxylin and eosin. (B) UMAP demonstrates the 16 cell clusters. Assigned cell types were labeled on the figure. TAMC, tumor-associated myeloid cell. (C) Violin plots showing the normalized expression level of cell type-specific markers. (D) Bubble plot showing the genes that were uniquely expressed in normal follicular and each cancer cell type. The degree of color represents the average expression value, and the size of dot represents the expression percentage in each cluster. (E) Volcano plot showing the DEGs between two types of PTC clusters (PTC1 and PTC2). The x axis represents log₂ fold change, and the y axis represents −log₁₀ P value. (F) Representative IHC staining for PBK in normal, PTC, and ATC tissues. (G) Left: Cell type proportion of cancer samples compared to a normal sample. The x axis represents each cell type, and the y axis represents the relative cell type proportion fold change compared to the normal sample. Blue and red boxes indicate PTC and ATC, respectively. Right: Deconvolution of bulk thyroid RNA-seq data (GSE33630) consisting of 45 normal samples, 49 PTC samples, and 11 ATC samples by using scRNA-seq data. The color indicates relative cell fraction between the samples.

Fig. 2. Evolution trajectory and transcriptional fluctuation during dedifferentiation of thyroid cancer.

(A) Distribution pattern of the cells derived from ATC samples. Red color indicates the cells derived from the ATC, and the circle indicates PTC clusters. (B) The t-Stochastic Neighbor Embedding (T-SNE) plot showing the 10 subclusters of PTC cell clusters (PTC1 and PTC2). On the left side, clusters of PTC cells were divided into 10 subclusters and each color represents each subcluster. On the right side, the cells derived from ATC samples are marked. Red color indicates the cells derived from the ATC samples. (C) Cell trajectory analysis of tumor progression from normal to cancer cells. Cells are labeled with each cell cluster: Follicular cells derived from normal (NOM_follicular), PTC (PTC_ Follicular), and ATC (ATC_ Follicular) samples and PTC1, PTC2, aPTC, and ATC cells. (D) Levels of TDS and ERKS (ERK score) on the trajectory. TDS consists of 16 genes and ERKS consists of 51 genes. Each score is represented by the average expression of the genes corresponding to each pathway. Color scale represents expression level. (E) PBK and LGALS3 expression levels on the trajectory. Color scale represents expression level. (F) IF staining showed a few ATC cells double-stained by PBK (red fluorescence) and galectin-3 (green fluorescence). Yellow circle indicates double-stained ATC cells. (G) Heatmap showing the expression changes of the 800 highly variable genes along the ATC-PTC axis of the trajectory. Significantly enriched functional annotations are shown on the right side of the heatmap. Color scale indicates z score. (H) Example GSEA plots showing the gene enrichment patterns of EMT and mTORC1 signaling genes on the ATC-PTC axis of the trajectory. Normalized enrichment score (NES) and false discovery rate (FDR) values are shown on the plots.

Fig. 3. CNAs during dedifferentiation of thyroid cancer.

(A) Large-scale CNAs of ATC and normal follicular cells identified using scRNA-seq data. Normal follicular cells were used as a reference. The CNA patterns for the normal follicular cells, ATC cells, and PTC cells from three ATC patients are shown. X axis, chromosome position; y axis, individual cells. Red color represents amplification, while blue color represents deletion. (B) CNA plots showing the relative CNA patterns in chromosome 11 that were identified by both scRNA-seq and WES data. X axis, position; y axis, copy number. (C) Bubble plot showing expression of the amplified genes in the Chr11-CNA. The location of the gene set is marked by box in (B). Dot sizes and color indicate the expression percentage and average expression, respectively. (D) Top: Cell trajectory analysis of tumor progression from normal to ATC cells. Cells are labeled with each cell cluster. Bottom: Average expression level of the genes in Chr11 CNA on the trajectory. Color scale represents expression level.

Fig. 4. Key transcriptional regulator of CREB3L1 in thyroid cancer evolution.

(A) Coexpression modules in each cell cluster and their master regulators were identified using SCENIC. Cell clusters and their origin samples are shown at the top of the figure. Color scale indicates regulon activity levels. (B) Ratio of samples (NOM, PTC, and ATC) between CREB3L1-positive and CREB3L1-negative in the ATC cell cluster. Color indicates each sample. (C) CREB3L1 expression level on the trajectory. Color scale represents expression level. (D) Staining of CREB3L1 in PTC sample and ATC and PTC parts from one concurrent ATC sample. Lowercase represents the ATC part, which was cytoplasm-positive in CREB3L1 with high density; uppercase represents the PTC part, which was also positive in CREB3L1. (E) Distribution pattern of the CREB3L1-binding peaks in different genomic elements. Genomic elements and percent of peaks are shown on the right side. (F) GSEA plot showing the gene enrichment patterns of CREB3L1 target genes on the trajectory from normal to ATC cells. Normalized enrichment score and FDR values are shown on the plots. (G) Bottom: CREB3L1 target gene enrichment scores for top significant hallmarks of GSEA analysis (FDR < 0.25). The x axis represents normalized enrichment scores. Top: Enrichment of overlapped genes between EMT- or mTORC1-related genes and CREB3L1 targets along normal thyroid tissue to ATC. (H) CREB3L1-binding and expression patterns of CREB3L1 target genes. Left: Genome browser screenshots showing the CREB3L1-binding patterns near COL1A1 and MMP11 gene loci. Right: Cell trajectory plots showing the expression levels of COL1A1 and MMP11.

Fig. 5. Involvement and functional status of CREB3L1 in thyroid cancer progression.

(A) GSEA plots showing the gene enrichment patterns of EMT and mTORC1 signaling genes between normal and ATC samples from four different publicly available bulk RNA-seq datasets. Normalized enrichment score and FDR values are shown on the right side of each figure. 1: GSE27155, 2: GSE33630, 3: GSE65144, and 4: GSE126698. (B) Volcano plot showing the DEGs between CREB3L1 high and low thyroid cancer groups. Red and blue indicate up-regulated genes in CREB3L1 high and low, respectively. The x axis represents log₂ fold change, and the y axis represents −log₁₀ P value. (C) GSEA plots showing the gene enrichment patterns of EMT, oxidative phosphorylation, KRAS, and mTORC1 signaling genes between CREB3L1 high and low thyroid cancer samples from TCGA. (D) Heatmap showing the expression patterns of EMT and mTORC1 signaling genes in CREB3L1 high and low groups. (E) Kaplan-Meier survival curves showing the overall survival and disease-specific survival in thyroid cancer patients from TCGA separated by extremely high CREB3L1 expression. (F) Comparison of positive ratio of CREB3L1 staining between recurrent group (2-year clinical recurrent PTC samples) and nonrecurrent group (10-year disease-free PTC samples). Representative images of CREB3L1 staining in both groups are shown below the x axis. (G) CREB3L1 overexpression significantly increased invasion ability of the TPC-1 cell line compared with control (P = 0.004). (H) Expression changes of EMT-related markers in TPC-1 cells with CREB3L1 overexpression versus control. (I) Schematic summary of the thyroid cancer progression revealed by this study.