Analysis of thyroid carcinoma composition and spatial architecture in the progression of dedifferentiation, lymphatic metastasis, and gastric metastasis

Abstract

Background

Gastrointestinal metastases are rare in patients with thyroid carcinoma (TC), and their underlying mechanisms remain unclear. Thus, in this study, we aimed to explore the spatial distribution characteristics of TCs and associated gastrointestinal metastatic cells.

Methods

We used spatial transcriptomics to generate an atlas that captures spatial gene expression patterns in papillary thyroid cancer (PTC), anaplastic thyroid carcinoma (ATC), ATC-associated lymphatic metastasis (ATC-LM), and rare ATC-associated gastric metastasis (ATC-GM).

Results

We demonstrated that tumor-specific myeloid cells with high SFRP4 expression were correlated with TC dedifferentiation and poor prognosis. Moreover, we validated their close localization to CD44⁺ tissue stem cells using immunofluorescence staining and spatial transcriptomics. We also demonstrated that ATC-LM and ATC-GM tissues exhibited high levels of CD44⁺PKHD1L1⁺ cells, which could serve as markers for these two pathological types.

Conclusions

These findings highlight the dynamic changes in cell composition, intercellular communication, and potential markers associated with TC dedifferentiation and distant metastasis. Further research based on our findings may contribute to improving diagnostic and therapeutic strategies for patients with TC.

Background

Thyroid cancer (TC) mainly includes aggressive papillary thyroid cancer (PTC) and anaplastic thyroid carcinoma (ATC). PTC, a well-differentiated type, has been observed at a higher rate than that of other pathological types. ATC is an undifferentiated type of cancer that exhibits stem cell-like properties, high proliferative potential, and resistance to current therapies. Histologically, PTC progresses to ATC via dedifferentiation, a biological process in pervasive cancers that induces the transition from a highly differentiated to a poorly differentiated status [1].

Dedifferentiated TC presents a risk of distant metastasis. The most common sites of ATC metastasis are the lungs, bones, and brain; however, gastrointestinal metastases (GM) are rare in patients with TC [1, 2]. The early diagnosis of GM in ATC is of great clinical significance, as it is crucial for guiding treatment and predicting prognosis. However, it faces numerous challenges, including nonspecific symptoms, variability in metastatic sites, imaging and molecular limitations, pathological complexity, individual variability, and comorbidities. At present, there is no well-established and effective method to accurately GM in ATC, largely due to the unclear understanding of its development mechanisms, and research progress is slow as samples are scarce [2].

The dedifferentiation and distant metastasis of TC constitute an interconnected process in which cancer cells gradually lose their differentiated characteristics, gain invasiveness, breach the basement membrane, enter the bloodstream or lymphatic system, and migrate to distant organs. To elucidate the underlying mechanisms, it is essential to obtain pathological samples from different stages (such as PTC, ATC, ATC associated metastasis samples) within the same individual for study, as they will comprehensively reflect the dynamic chain of disease progression.

In this study, we report a patient (patient #6) with ATC-PTC mixing (ATC differentiated from PTC), accompanied by rare simultaneous GM and lymphatic metastases (LM). These lesions from the same patient represent different stages of TC progression, which is considerably beneficial for studying TC progression.

Single-cell RNA-sequencing (scRNA-seq) transcriptomes and recently developed spatial transcriptome (ST) technology have been used to identify the mechanisms of cancer development [3,4,5]. Therefore, we used formalin-fixed paraffin-embedded (FFPE) samples obtained from various lesions of patient #6, together with nodular goiter tissues and primary PTC and ATC tissues from other included patients, to study the spatial characteristics of specific tissue types based on ST analysis [5].

Methods

Human specimens

All study participants visited the Qingdao Municipal Hospital (Qingdao, China) during 2020–2022. The first cohort included eight FFPE specimens (Additional file 1), which were used for ST analysis. The second cohort included 11 FFPE specimens of primary single PTC (from 11 additional patients with PTC), eight FFPE specimens of primary single ATC (from eight additional patients with ATC), six FFPE specimens of primary gastric cancer (GC), and six paired lymph node metastasis (GC-LM) FFPE specimens (from six patients with GC-LM) (Additional files 2 and 3). Specimens from both cohorts were used for immunofluorescence (IF) validation. Two experienced pathologists independently confirmed the sample histology in a blinded manner and classified them according to the eighth edition of the AJCC TNM (American Joint Committee on Cancer Tumor Node Metastasis Classification) system. The procedure was approved by the Ethics Committee of Qingdao Municipal Hospital (approval number: 2023linshenzidi107) and performed in accordance with the Declaration of Helsinki. Written informed consent was obtained from all the study participants.

ST analysis

RNA quality was assessed by calculating the RNA integrity number (RIN) of freshly collected tissue sections (RIN should be ≥ 7). FFPE sections at 5 µm were prepared on a Visium Spatial Tissue Optimization Slide for methanol fixation and hematoxylin and eosin staining (10X Genomics). Polyadenylated mRNA released from the overlying cells was captured using primers on the spots. RT Master Mix containing reverse transcription reagents was added to the permeabilized tissue sections. The second-strand mix was added to the tissue sections on the slide to initiate second-strand synthesis. This was followed by the denaturation and transfer of cDNA from each capture area to a corresponding tube for amplification and library construction. After transferring cDNA from the slide, spatially barcoded full-length cDNA was amplified via PCR to generate sufficient mass for library construction. After library construction, 150PE mode sequencing was performed using the Illumina NovaSeq600 platform. For resolution, each spot is 55 µm in diameter with a 100 µm center to center distance between spots.

Data processing

The Space Ranger software provided by 10 × Genomics (version 1.3.1) was used to perform sample demultiplexing, alignment, tissue detection, fiducial detection, and UMI counting [6].

Gene ontology enrichment and kyoto encyclopedia of genes and genomes enrichment analysis

ClusterProfiler was used for the enrichment analysis [7], which shows only the enrichment graphs of the different clusters with graph-based clustering results.

Gene set variation analysis (GSVA)

GSVA was performed on 50 hallmark pathways annotated in the molecular signature database, as previously described [6].

Analysis of cell–cell communications

CellChat objects were created based on the UMI count matrix for each group using CellChat (https://github.com/sqjin/CellChat, R package, v.1) [8]. With “CellChatDB.human” set up as the ligand-receptor interaction database, cell–cell communication analysis was performed using the default settings. A comparison of the total number of interactions and interaction strength was obtained by merging the CellChat objects of each group using the function mergeCellChat. The differential number of interactions or interaction strengths among different cell populations were visualized using netVisual_diffInteraction.

IF staining

IF staining was performed as described in a previous publication [9]. FFPE human tumor and control specimens were cut into 3–4 mm-thick sections. Endogenous peroxidase activity and nonspecific Fc receptor binding were blocked with H₂O₂ 3% (Gifrer, 10603051) and Protein Block (Dako, X0909), respectively. After mounting on EcoMount (Biocare Medical, EM897L), the slides were scanned using a NanoZoomer (Hamamatsu). CD14/CD44/SFRP4 3-plex staining and CD44/PKHD1L1 2-plex staining were performed manually on FFPE tumor sections using tyramide system amplification (TSA) and a conventional fluorescent dye-conjugated secondary antibody system (all secondary antibodies were diluted 1:100). Nuclei were stained with DAPI Solution (Thermo Fisher Scientific, Waltham, MA, USA, 62248) at 2 mg/ml for 10 min. After mounting with ProLongTM Gold Antifade Mountant (Thermo Fisher Scientific, P36934), the slides were scanned at 20 X magnification using a Zeiss Axio scan Z1 device. The images were processed using ImageJ to calculate the fluorescence intensity associated with the targets. The quantification of fluorescence intensity is expressed as the mean fluorescence intensity per unit area.

The university of ALabama at birmingham CANcer data analysis portal (UALCAN) for validation

The relative expression of target genes in specific cancer species was analyzed based on the UALCAN analysis of data from The Cancer Genome Atlas (TCGA) [10].

Statistical analysis

Data are expressed as mean ± standard deviation. All experiments were performed in triplicates. SPSS 16.0 software (SPSS, Inc., Chicago, IL, USA) was used to conduct two-tailed t-tests to determine the differences between groups, and statistical significance was set at p < 0.05.

Results

Case report of patient #6

To the best of our knowledge, only sporadic publications have reported cases of ATC-GM [2, 11]. The present study was conducted based on the discovery of a rare case, patient #6; from her tissue samples, we obtained a mixed cancer type of PTC and ATC, as well as LM and GM.

The tumor was staged as pT3N1aM0. Two months after surgery, a cervical mass, confirmed to be metastasis to the cervical lymph nodes, was identified (Additional file 4a). Subsequently, the patient underwent radiotherapy targeting the affected cervical lymph nodes.

Four months after the surgery, CT revealed a mass in the gastric antrum (Additional file 4b). Gastroscopic examination indicated a large deep ulcer-like lesion at the gastric angle and antrum (Additional file 4c). Biopsy and immunohistochemistry revealed metastatic undifferentiated carcinoma originating from the thyroid in the gastric antrum and angle.

ST profiling

Next, we conducted ST analysis to visualize molecular changes in the transition from PTC to ATC and its associated LM and GM. In the first cohort, we included one ATC-PTC mixed tissue, one ATC-LM tissue, and one ATC-GM tissue from patient #6, together with two nodular goiter samples (from patients #1 and #2), two primary single PTC samples (from patients #3 and #4), and one primary single ATC sample (from patient #5) to perform ST analysis (Additional files 1 and 5a). The ST results revealed decreased TTF1, TG, and PAX8 expression levels in TC compared to that in the nodular goiter, and these expression levels were even lower in ATC areas than in PTC areas (Additional file 5b–d), consistent with a previous report [12]. IGF2BP1, which is a specific marker for ATC [13], was almost negatively expressed in tissues other than ATC (Additional file 5e). These results confirm the accuracy of the pathological typing.

Subsequently, the cells in these samples were annotated into eight major cell types according to the Celldex database [7] (Fig. 1a–c, Additional file 6a–c), and each cell type is visually and quantitatively presented in Additional file 7. Nodular goiter tissues exhibited high levels of endothelial cells, epithelial cells, and fibroblasts, and extremely limited amounts of T, natural killer (NK), B, and myeloid cells (Fig. 1c, Additional file 7). In TC tissues, the proportions of endothelial cells, epithelial cells, and fibroblasts were all decreased, while increased myeloid cell infiltration was observed, compared to those in the nodular goiter tissues (Fig. 1c, Additional file 7). Notably, myeloid cells were recruited to the ATC area in ATC-PTC mixed tissues (Sect. 6), whereas B, T, and NK cells exhibited uniform distribution in the ATC and PTC areas (Fig. 1c, Additional file 7).

Fig. 1

Expression and spatial distribution of major cell types in TC revealed by spatial transcriptomics. a Dot plots showing the average expression of known markers in indicated cell clusters. The dot size represents the percentage of cells expressing the genes in each cluster. b tSNE (left) and UMAP (right) plots of cells from all the tested samples, showing eight clusters in each plot. Each cluster was shown in a different color. c The spatial distribution of different cell types in representative pathological areas (ST maps for all eight FFPE samples are shown in Additional files 7 and 8). Left: Hematoxylin and eosin (H&E) staining of representative tissue sections; Right: Unbiased clustering of ST spots and defined cell types of each cluster in the corresponding tissue sections. d Pseudotime trajectory of major cell types analyzed by Monocle. e Volcano plots displaying the upregulated and downregulated genes in specific cell types between PTC and ATC. f Dot plots of hallmarks for differentially expressed genes in the global cell type between PTC and ATC. The intensity represents the average fold change of gene expression in PTC versus ATC. The dot size shows FDR for each hallmark. ATC anaplastic thyroid cancer, GO Gene Ontology, KEGG Kyoto Encyclopedia of Genes and Genomes, PTC papillary thyroid cancer, ST spatial transcriptomics, TC thyroid cancer, UMAP uniform manifold approximation and projection, FDR false discovery rate

The ATC-LM tissue (Sect. 7) was mainly composed of T, NK, and B cells. Additionally, some myeloid cells were present within the tumor area of Sect. 7 (Fig. 1c, Additional file 7). However, few myeloid cells were observed in the ATC-GM tissue (Sect. 8, Fig. 1c Additional file 7). The number of tissue stem cells in TC tissues was significantly higher than in nodular goiter tissues with a substantial variability observed in metastatic cancers (Fig. 1c, Additional file 7). A trajectory analysis suggested that during TC development, tissue stem cells were gradually activated and transformed into immune cells, such as myeloid cells, T, NK, and B cells (Fig. 1d).

CellChat predicts WNT signaling as the mediator of tissue stem cells and myeloid cell interaction changes in TC dedifferentiation

Integrative Gene Ontology term results and Kyoto Encyclopedia of Genes and Genome enrichment demonstrated that myeloid cells and tissue stem cells shared the most similar functions, which were mostly involved in cellular communication (Fig. 1e, f). Therefore, we performed CellChat analysis [8] to investigate intercellular communications and the associated pathways, especially between myeloid cells and tissue stem cells (we define it as “M-S-inter”). The signal strengths emitted by tissue stem cells, myeloid cells, and fibroblasts were stronger than those in other cells (Fig. 2a, Additional file 8). Moreover, we identified that the main outgoing signaling in tissue stem cells, i.e., myeloid cells, was characterized by pattern #1, which represents multiple pathways, such as COLLAGEN, LAMININ, and WNT. The main incoming signaling in tissue stem cells is characterized by pattern #1, which also includes the WNT pathway (Fig. 2b).

Fig. 2

Analysis of cell–cell communication identifies major signaling changes in PTC and ATC patients. a Circle plots depict the global interaction strength in the cell–cell communication network from all the tested samples. b The inferred outgoing communication patterns of secreting cells (left), and the inferred incoming communication patterns of target cells (right), show the correspondence between the inferred latent patterns and cell groups, as well as signaling pathways. The thickness of the flow indicates the contribution of the cell group or signaling pathway to each latent pattern. c Bar plots showing the number of inferred interactions (left) or interaction strength (right) in the cell–cell communication network analyzed via CellChat across the PTC and ATC. d Heatmaps of the differential number of interactions between PTC and ATC, showing the outgoing and incoming signaling change of each cell group in greater detail (the top-colored bar plot represents the sum of each column of values displayed in the heatmap (incoming signaling). The right-colored bar plot represents the sum of each row of values (outgoing signaling). e Comparison of the signaling pathway based on the relative information flow between pairwise datasets. f The heatmap shows the relative importance of each cell group based on the computed four network centrality measures of the WNT signaling network. g Circle plots showing differential WNT signaling networks between PTC and ATC. ATC anaplastic thyroid cancer, PTC papillary thyroid cancer

The total interaction numbers and strengths weakened from PTC to ATC (Fig. 2c). The number of “M-S-inter” exhibited the most significant difference (Fig. 2d). Notably, among the differential signaling pathways between PTC and ATC, the WNT pathway ranked high (Fig. 2e). These predicted results suggest that the WNT pathway may play an important role in the transition from PTC to ATC by mediating “M-S-inter” (Fig. 2f, g).

Identification of specific subtypes of myeloid cells and tissue stem cells in TC dedifferentiation

We subsequently identified five novel and distinguishable subpopulations of myeloid cells by clustering the sequencing results and identifying the markers of each subpopulation (Fig. 3a, b). Gene set variation analysis (GSVA) revealed that these five clusters had different levels of signaling pathways, indicating that they have different functions in tumor development (Fig. 3c). We also identified eight novel distinguishable clusters of tissue stem cells with specific markers (Fig. 3d, e) and different signaling pathway phenotypes (Fig. 3f). We mapped these clusters of myeloid cells and tissue stem cells onto the hematoxylin and eosin plane and found that the myeloid cells recruited to the ATC were generally in cluster 1, whereas tissue stem cells in the ATC area were mainly in cluster 4 (Fig. 3g, Sect. 6 as a visual representation), which are both considerably different from the PTC area (Fig. 3g). A subset of upregulated genes was shared between cluster 1 of myeloid cells and cluster 4 of tissue stem cells, including PLTP, SFRP4 (WNT antagonist), and COL3A1, suggesting a commonality (Fig. 3h, Additional files 9a–c and 10a–c).

Fig. 3

Identification of specific subtypes of myeloid cells and tissue stem cells in TC dedifferentiation. Myeloid cells are clustered into five subpopulations based on the specific marker expressions. a Heatmap showing the expression of marker genes in the indicated subpopulations; b tSNE (left) and UMAP (right) plots of subpopulations from all the tested samples. Each subpopulation was shown in a different color; c differences in pathway activities scored per cell via GSVA between the different subpopulations. Tissue stem cells are clustered into eight subpopulations based on the specific marker expressions. d Heatmap showing the expression of marker genes in the indicated subpopulations; e tSNE (left) and UMAP (right) plots of subpopulations from all the tested samples. Each subpopulation is shown in a different color; f differences in pathway activities scored per cell via GSVA between the different subpopulations. g The spatial distribution of different subpopulations in indicated pathological areas. Left: Hematoxylin and eosin (H&E) staining of representative tissue sections; Right: Unbiased clustering of ST spots and defining cell types of each subpopulation in the corresponding tissue sections. h Spatial feature plots of gene expression of PLTP, COL3A1, and SFRP4 in representative tissue sections. 3-plex immunofluorescence staining of human PTC and ATC tissue. i Representative immunofluorescence staining of CD44 (orange), DAPI (blue), SFRP4 (green), and CD14 (red), in individual and merged channels are shown; j Statistical differences in different fluorescence intensities or their merge rates. TC thyroid carcinoma, ATC anaplastic thyroid cancer, PTC papillary thyroid cancer

Analysis of the TCGA database indicated that SFRP4 expression was significantly lower in TC tissues than in normal controls, and high SFRP4 levels predicted poor prognosis (Additional file 11a, b). In contrast, PLTP and COL3A1 were not associated with the prognosis of patients with TC (Additional file 11c–f), suggesting the importance of the WNT antagonist SFRP4 in the development of TC. Therefore, we subsequently focused on SFRP4 and performed 3-plex IF imaging on sections from PTC and ATC tissues in cohort 2 (Additional file 2). The myeloid cell marker CD14 and tissue stem cell marker CD44 demonstrated a noticeably wide distribution pattern in both PTC and ATC areas. CD14 and CD44 exhibited limited spatial overlap (“myeloid stem cells”), with the majority primarily clustered near each other, providing a microenvironment conducive to their interactions. SFRP4 was more highly expressed in ATC areas than in PTC areas, with most SFRP4 overlapping with CD14, indicating that SFRP4 was mostly expressed in myeloid cells (CD14⁺SFRP4^high subtype, Fig. 3i, j).

Identification of specific subpopulations of tissue stem cells in ATC-GM

Tissue stem cell subtypes were mainly present in cluster 6 in the ATC-GM area in Sect. 8, which was different from that in primary ATC (Sect. 6) but similar to the phenotype in ATC-LM (Sect. 7) (Fig. 4a, b). The differentially expressed gene RNA-sequencing results demonstrated significant differences in gene expression between the primary ATC area and both the ATC-LM and ATC-GM areas (Fig. 4c), with a considerable portion displaying a consistently differential trend (Fig. 4d).

Fig. 4

Identification of specific subpopulations of tissue stem cells in ATC-GM. a The spatial distribution of different subpopulations in indicated pathological areas. Left: Hematoxylin and eosin (H&E) staining of representative tissue sections; right: Unbiased clustering of ST spots and defined cell types of each subpopulation in the corresponding tissue sections. b tSNE (up) and UMAP (down) plots of subpopulations from primary ATC, ATC-LM and ATC-GM samples. Each subpopulation is shown in a different color. Differentially expressed gene RNA-sequencing analysis between ATC and ATC-LM and between ATC and ATC-GM. c Heatmap showing the differentially expressed genes between ATC and ATC-LM (left), and ATC and ATC-GM (right); d Venn diagram showing the overlapped genes described in (c). e Spatial feature plots of gene expression of PKHD1L1 and TACSTD2 in tissue sections. f 2-plex immunofluorescence staining of human primary ATC, ATC-LM, and ATC-GM tissues: Top: Representative immunofluorescence staining of CD44 (red), DAPI (blue), and PKHD1L1 (green), in individual and merged channels are shown; Bottom: Statistical differences in different fluorescence intensities or their merge rates. g 2-plex immunofluorescence staining of human primary GC, and GC-LM tissues: Top: Representative immunofluorescence staining of CD44 (red), DAPI (blue), and PKHD1L1 (green), in individual and merged channels are shown; Bottom: Statistical differences in different fluorescence intensities or their merge rates. ATC anaplastic thyroid cancer, GC gastric cancer, GM gastric metastasis, LM lymphatic metastasis

Cluster 6 tissue stem cells, distributed in both ATC-LM and ATC-GM, specifically overexpressed PKHD1L1 and TACSTD2 (Additional file 12a–c). Notably, PKHD1L1 expression was low in both PTC and ATC and was specifically overexpressed in ATC-LM and ATC-GM (Fig. 4e). Analysis of the TCGA database revealed that PKHD1L1 expression was significantly lower in TC tissues than in normal controls, and high PKHD1L1 levels predicted a poor prognosis in patients with TC (Additional file 13a, b). PKHD1L1 was also weakly expressed in primary GC, although this expression did not significantly correlate with prognosis (Additional file 13c, d). Additionally, PKHD1L1 expression remained unchanged in most cancer types, whereas it was expressed at significantly lower levels in breast cancer (BRCA), colon adenocarcinoma (COAD), lung adenocarcinoma (LUAD), and lung squamous cell carcinoma (LUSC) (Additional file 13e). These results suggest that PKHD1L1 may specifically serve as a marker for metastatic cancer.

Subsequently, 2-plex IF imaging was performed to explore the distribution of cluster 6 tissue stem cells (CD44⁺PKHD1L1^high cells) in sections of different pathological sites. By comparing Sects. 6, 7, and 8, we confirmed that the number of CD44⁺PKHD1L1^high cells was significantly lower in the primary ATC area and increased in the ATC-LM and ATC-GM (Fig. 4f). Further, CD44⁺PKHD1L1^high cells were at present at extremely low levels in primary GC tissues and were significantly increased in GC-LM tissues (Fig. 4g). Collectively, these results suggest that the proportion of CD44⁺PKHD1L1^high cells may specifically increase in metastatic cancer and serve as a novel biomarker for ATC-LM and ATC-GM.

Discussion

The dedifferentiation of TC (From PTC to ATC) is closely linked to distant metastasis [14]. ScRNA-seq is a useful tool for exploring the molecular mechanisms of TC dedifferentiation and its associated metastasis. Luo et al. [3] and Lu et al. [15] reported that during the transition from PTC to ATC, cancer-associated fibroblasts exhibited strong interactions among mesenchymal cell types, macrophages shifted from M1 to M2 states, and T cells were reprogrammed from cytotoxic to exhausted states, highlighting new therapeutic opportunities for ATC treatment. Luo et al. [3] identified a potential functional role of CREB3L1 in TC dedifferentiation by integrated analyses of copy number alterations and transcriptional regulatory networks. However, scRNA-seq loses spatial information, which may be particularly important for ATC transformation. In this study, we performed ST analysis to visualize molecular changes in the transition from PTC to ATC and its associated LM and GM.

We observed that the specific macrophage type CD14⁺SFRP4^high increased during the PTC-to-ATC conversion. Macrophages secrete cytokines such as IL-6 and TGF-β, which activate signaling pathways like STAT3 and Smad, leading to tumor dedifferentiation by supporting the stemness of cancer cells [16, 17]. SFRP4 acts as a secreted antagonist of WNT signaling [18], which was closely related to macrophages functions [16]. Based on the CellChat prediction and IF analysis, we identified that the CD14⁺SFRP4^high subtype could communicate with tissue stem cells through the WNT pathway, and this communication was significantly enhanced in ATC, compared to that in PTC. Combined with literature reports above and our analysis, this study strongly suggests that CD14⁺SFRP4^high may regulate the dedifferentiation and its associated distant metastasis process through communication with tissue stem cells. Considering that CD14⁺SFRP4^high is not unique to TC, this might be a common molecular mechanism for pan-cancer metastasis, which could also help understand metastasis in other cancer types.

And, we demonstrated that high level of SFRP4 predicted poor prognosis in TC patients based on the TCGA database, which suggests that SFRP4 plays a promoting role in TC progression. To the best of our knowledge, SFRP4 plays a complex role in tumorigenesis. Generally, it is considered to inhibit cell proliferation, and promote apoptosis [18]. However, there is emerging evidence that suggests the SFRP4 may have alternate functions related to the promotion of carcinogensis in a context and tumor type-dependent manner. For instance, in uterine leiomyomas, SFRP4 overexpression, regulated by progesterone, enhances the proliferation of smooth muscle cells, thereby driving tumor progression [19]. In gastric cancer, hypomethylation of the SFRP4 promoter leads to its upregulation, which has been associated with accelerated tumor growth [20]. These results partially explain the biological regulatory role of SFRP4 in tumor development and poor prognosis.

Although PKHD1L1 is primarily associated with polycystic kidney and liver disease [21], its role in cancer has not been well established. Some studies have suggested a potential involvement in tumor progression and metastasis through its association with cellular adhesion, migration, and invasion [22]. Through analysis of TCGA database, we found that PKHD1L1 generally has low or no abnormal expression in various tumor cells, and its high expression predicts poor prognosis. These results suggest that PKHD1L1 might promote tumor progression and metastasis. The ST data and IF detection results showed that the CD44⁺PKH1L1^high subtype is mainly present in GM, which is substantially different from primary ATC or PTC and is similar to the phenotype of ATC-LM tissues. It is widely accepted that there is a sequential progression in tumors, in which the primary tumor first seeds LM, and LM further seeds distant metastases [23]. And, CD44 mediates the entry of CD44 + tumor cells into the lymphatic system to form lymphatic metastasis [24]. Therefore, we speculate that CD44⁺PKH1L1^high cells may have stronger lymphatic metastasis efficacy, leading to the preferential transfer of CD44⁺PKH1L1^high cells from ATC to lymph nodes, which explains the observed phenomenon that CD44⁺PKH1L1^high cells are scarce in ATC but abundant in ATC-LM and ATC-GM. Further studies are required to elucidate this specific role and its regulatory mechanisms.

In this study, we analyzed several existed TC markers and compared them with our findings. TTF1 and TG are typically expressed in tumors originating from thyroid follicular epithelial cells, while PAX8 is also a marker of thyroid follicular cells. However, the expression of these markers may vary in certain types of TC. For instance, in ATC, TTF1 and TG are usually negative, while PAX8 is positive in approximately half of the cases [25]. This study also found that these three markers are lowly expressed in ATC, which is not conducive to distinguishing ATC tissues from normal tissue cells. Haase et al. identified IGF2BP1 as the first positive marker for ATC diagnosis [13], and our study here confirms its reliability. Here, as the new markers we identified, CD14⁺SFRP4^high is specific over-expressed in ATC, while CD44⁺PKH1L1^high is specific over-expressed in ATC and ATC-LM and ATC-GM, and decreased in ATC. Further, these two markers were both highly associated with TC prognosis, which may serve as potential indicators for treatment outcomes. Overall, our findings, along with existing biomarkers, each have their own advantages. In clinical practice, disease detection sensitivity and specificity may be improved through appropriate combined testing methods.

There are several limitations to this study. (1) Limitations of ST technology: the key challenge is spatial resolution, as current methods may not capture fine-scale details of tissue architecture, especially in complex or heterogeneous samples. Additionally, batch effects can introduce variability between samples due to differences in processing, reagents, or platform-specific issues, potentially skewing results. Variability in tissue sampling, such as differences in collection or handling, can further contribute to inconsistencies, affecting the representativeness and accuracy of data. In addition, ST in a clinical setting is challenged by high equipment and resource expenses, limited availability of specialized facilities and expertise, and variability in tissue processing that can impact result consistency. These all need to be overcome and resolved in the subsequent technological development. (2) Limitations of sample size: This study is based on the analysis of a rare case, which might not reflect the overall characteristics of the population with this type of disease. Therefore, we need to validate the key findings from this study in larger cohorts in future work and explore the mechanisms underlying different TC subtypes. (3) Limitations of molecular mechanism elucidation: we did not elucidate how the specific cell type (including CD14 + SFRP4high and CD44 + PKH1L1high) contribute to the tumor microenvironment clearly. The further molecular mechanism exploration should be performed in the subsequent study.

Conclusions

The significance of this study is that it provides a marker for the differential diagnosis between primary ATC and its distant metastasis based on the differential expression of PKHD1L1.