Spatial transcriptomics reveals prognosis-associated cellular heterogeneity in the papillary thyroid carcinoma microenvironment

Abstract

Background

Papillary thyroid carcinoma (PTC) is the most common malignant endocrine tumour, and its incidence and prevalence are increasing considerably. Cellular heterogeneity in the tumour microenvironment is important for PTC prognosis. Spatial transcriptomics is a powerful technique for cellular heterogeneity study.

Methods

In conjunction with a clinical pathologist identification method, spatial transcriptomics was employed to characterise the spatial location and RNA profiles of PTC-associated cells within the tissue sections. The spatial RNA-clinical signature genes for each cell type were extracted and applied to outlining the distribution regions of specific cells on the entire section. The cellular heterogeneity of each cell type was further revealed by ContourPlot analysis, monocle analysis, trajectory analysis, ligand–receptor analysis and Gene Ontology enrichment analysis.

Results

The spatial distribution region of tumour cells, typical and atypical follicular cells (FCs and AFCs) and immune cells were accurately and comprehensively identified in all five PTC tissue sections. AFCs were identified as a transitional state between FCs and tumour cells, exhibiting a higher resemblance to the latter. Three tumour foci were shared among all patients out of the 13 observed. Notably, tumour foci No. 2 displayed elevated expression levels of genes associated with lower relapse-free survival in PTC patients. We discovered key ligand–receptor interactions, including LAMB3–ITGA2, FN1–ITGA3 and FN1–SDC4, involved in the transition of PTC cells from FCs to AFCs and eventually to tumour cells. High expression of these patterns correlated with reduced relapse-free survival. In the tumour immune microenvironment, reduced interaction between myeloid-derived TGFB1 and TGFBR1 in tumour focus No. 2 contributed to tumourigenesis and increased heterogeneity. The spatial RNA-clinical analysis method developed here revealed prognosis-associated cellular heterogeneity in the PTC microenvironment.

Conclusions

The occurrence of tumour foci No. 2 and three enhanced ligand–receptor interactions in the AFC area/tumour foci reduced the relapse-free survival of PTC patients, potentially leading to improved prognostic strategies and targeted therapies for PTC patients.

INTRODUCTION

Papillary thyroid carcinoma (PTC) is the most common malignant tumour of the endocrine system, and its incidence has been rapidly increasing in China over the past decade.^1–3 While the prognosis for most PTC patients is favourable, a substantial proportion of patients still experience progressive invasive disease. The prognosis and therapeutics of PTC can be affected by several factors, including tumour differentiation and dedifferentiation,^4–7 atypical follicular cells (FCs) (AFCs) (atypical cells with focal nuclear changes, extensive but mild nuclear changes)^8–12 and immune cells^13,14 in the tumour microenvironment.¹⁵ Therefore, it is crucial to explore the cellular heterogeneity in the PTC tumour microenvironment. Generally, the cell type and tumour area in PTC tissues are defined by a pathologist using classic pathological examination.^16,17 However, due to the variations in pathologists' expertise and insufficient spatial resolution of visual methods, pathological examinations are insufficient to fully elucidate the cellular heterogeneity of PTC, especially the discovery of valuable atypical tumour cells. Hence, methods with high precision and throughput, such as single-cell RNA sequencing (scRNA-seq),^6,7 transcriptomics (RNA-seq) and proteomics analyses,^5,18–20 have been gradually applied to cell heterogeneity studies in cancers, including PTC. For example, researchers characterised the landscape underlying the initiation and progression of thyroid cancer,^4–6,21 the dedifferentiation of thyroid cancer⁷ and the immune landscape of PTC^6,22 using scRNA-seq. They provided new insights into the understanding of PTC, but the omission of tissue in situ information limited these research efforts.

Spatial transcriptomics (SPT) has emerged as a powerful tool for investigating cellular heterogeneity based on spatial gene expression.^23–26 SPT has been extensively applied in the study of organ development in humans,^27–32 hematopoietic stem cells^33–35 and various types of tumours.^36–40 However, SPT has not been used to clarify the spatial cellular heterogeneity of PTC, particularly the spatial distribution of tumour cells, immune cells and normal tissue cells within or surrounding tumour foci, which significantly influence cellular functions and intercellular dynamics within tissues.^41–43

Analytical approaches for SPT data include clustering of spatial spots based on gene expression patterns,²⁵ identification of tumour areas through pathological observations and extraction of clusters using t-distributed stochastic neighbour embedding (tSNE) for dimensionality reduction analysis.³⁹ Mapping cell subpopulations from scRNA-seq onto SPT enables the identification of accurate spatial cell populations,^29,37 and signature genes derived from scRNA-seq clusters can be used to determine spatial cell types.⁴⁰ However, these previous studies have primarily focused on RNA-based analyses without integrating spatial location information. Moreover, these existing analytical methods have certain limitations, such as focusing solely on RNA information and lacking integration with clinical pathology data provided by pathologists. This may result in false positive and false negative results to some extent. To address these limitations, our research aims to combine spatial transcriptome analysis with clinicopathological diagnosis, leveraging the strengths of RNA-based clustering and clinical pathology information.

In this study, to characterise the spatial distribution and RNA profiles of PTC-associated cells, we collected thyroid cancer tissues from four PTC patients and used SPT in conjunction with clinical pathology identification methods. By extracting spatial RNA-clinical signature genes for each cell type, we outlined the distribution regions of specific cells throughout the tissue sections. Various analytical approaches, including ContourPlot analysis, monocle analysis, trajectory analysis, ligand‒receptor analysis and gene set enrichment analysis, were employed to reveal the cellular heterogeneity within each cell type. Moreover, special tumour foci and three enhanced ligand‒receptor interactions within the AFC area/tumour foci were found to be associated with the prognosis of PTC patients. By integrating RNA-based clustering with clinical pathology information, our method offers improved accuracy and comprehensive insights into spatial transcriptome analysis and prognosis-associated cellular heterogeneity in the PTC microenvironment, which provides valuable implications for prognosis evaluation in PTC patients.

MATERIALS AND METHODS

Human specimens

All procedures and the use of tissue were performed according to recent national ethical guidelines. Five patients diagnosed with PTC who underwent surgery at the Second Affiliated Hospital of the South China University of Technology were included. Table S1 provides a summary of the clinical information of the patients.

Sample preparation and selection

Plastic embedding moulds (8 mm × 8 mm) were filled with OCT (TissueTek Sakura) at room temperature. Fresh samples were washed twice with DMEM (Gibco), and areas of interest with a base area of less than 6 mm × 6 mm were excised. The excised tissues were then embedded in OCT, frozen using dry ice or liquid nitrogen (without direct contact with liquid nitrogen) and stored at −80°C in 50 mL centrifuge tubes for subsequent spatial experiments.

Fresh frozen samples, stored for less than 3 months, were cryosectioned at a thickness of 10 μm at −15°C using a Leica cryostat. The sections were placed on highly adhesive slides (Citotest) for sample selection. The tissue sections were fixed in prechilled methanol (−20°C), reactivated with isopropanol and subjected to an improved haematoxylin and eosin (H&E) staining method. The staining procedure included haematoxylin staining for 1 min, bluing for 15 min, 80% ethanol for 2 min, eosin staining for 1 s, 80% ethanol for 2 min, 95% ethanol for 2 min (twice), absolute ethanol for 2 min (twice) and environmental protection transparent agent for 2 min (twice).

Histology images of the stained sections were captured using an Aperio Imaging and Analysis System (Leica) with a 20× objective. Valid samples were selected for further 10× Visium spatial experiments based on their pathological features.^25,26,31

SPT sequencing

The selected frozen samples were cryosectioned and carefully placed on Visium Tissue Optimization Slides or Visium Spatial Gene Expression Slides (10x Genomics). To enhance the clarity and accuracy of histological information, mirror slice slides synchronised with the Visium spatial slides were prepared, where the front sides of the two slides represent adjacent areas. The spatial tissue sections were fixed, stained and scanned following the 10x Genomics Visium Spatial protocols (10× Genomics Document CG000241), while the synchronised mirror slices were stained using the improved H&E staining method (as described in Section 2.2). Through tissue optimisation time-course experiments, it was determined that a permeabilisation time of 12 min was optimal for PTC samples with a thickness of 10 μm. Libraries were prepared using the Visium Spatial Gene Expression workflow (10× Genomics Document CG000239). The quality and integrity of the library DNA fragments were assessed using an Agilent 2100 instrument. The qualified libraries were then subjected to sequencing on the Illumina NovaSeq 6000 platform.^25,26,31

SPT data processing

The raw images were exported as TIFF files using Aperio ImageScope Software (Version 12.3) with both low- and high-resolution settings. The areas containing spatial barcodes were identified and marked on Loupe Browser Software (Version 6.0), and the marked areas were exported as JSON files. The Visium spatial RNA-seq output FASTQ files, spatial barcodes and histology images were processed using Space Ranger software (Version 1.1.0) to generate feature-spot matrices.^25,26,31

For SPT data processing, the feature-spot matrices containing histology features for each sample (as obtained from Visium raw data processing) were loaded into R (Version 4.0.5) using the Seurat R package (Version 4.0.2).^25,26,31 The number of UMIs (nUMI) and genes (nGene) detected at each spot was calculated initially. Quality control (QC) was performed by filtering out genes expressed in less than three spots, mitochondrial genes and ribosomal genes. After QC, the IntergrateData function of the Seurat R package was utilised to integrate the expression data and remove batch effects (canonical correlation analysis [CCA] and RPCA, canonical correspondence analysis and reciprocal principal component analysis [PCA]) across different patients. The integrated expression matrix was then normalised, log-transformed, centred and scaled to calculate the common variable genes. PCA of the 2000 highly variable genes was conducted to obtain a low-dimensional space. The top five principal components (PCs) were selected as relevant dimensions for performing nearest-neighbours-based clustering and UMAP (uniform manifold approximation and projection) or tSNE visualisation analysis. Consensus subgroups of patients were defined based on RNA expression matrices, followed by analysis using individual sample RNA signatures combined with histological features.

Spatial-geneset-score processing

Spatial-geneset-score processing involved using cell-type-based signature genes as a gene set.^40,44 The spatial-geneset-score was calculated as the average log-transformed normalisation values of the signature genes. Since the current resolution of Visium SPT cannot assign specific cell types, the spatial-geneset-score of each spot was used as the score for defining the corresponding cell type.

Cell type identification at the transcriptomic and pathological levels

The differentially expressed genes of all identity clusters were identified by the FindAllMarkers function of the Seurat R package. The top 20 markers were used to identify the preliminary cell types at the transcriptome level.⁷ Because of the limitation of current Visium SPT resolution, each spot contains approximately 3−20 cells. We developed an improved signature-based strategy that combines the advantages of the transcriptome and histology approaches. To obtain an accurate definition of cell type at the pathological level, Prof. Liu, a clinical pathologist of PTC (KingMed Diagnostics, Guangzhou, China), was invited to distinguish and label the precise spots with pure tumour cells (enlarged/elongated/overlapped nuclei, transparent/marginalised/glassy chromatin, irregular nuclear outlines, nuclear grooves and nuclear pseudoinclusions), FCs (normal thyroid follicular epithelium cells), AFCs (slightly enlarged, overlapping, lightly stained nuclei, but no obvious nuclear grooves, nuclear inclusions and other PTC nuclear features) and immune cells (reactive hyperplastic lymphocytes).⁴⁵ Then, ‘pure’ spots that met both transcriptomic and pathological criteria were used to screen cell type-based signature genes of three cell types (tumour cells, FCs including AFCs, and immune cells) (Table S2). Next, the spatial-geneset-score of 3 cell types in each spot was calculated based on the top 20 cell-type-based signature genes, respectively. The spatial-geneset-score was calculated as the average log-transformed normalisation values of the signature genes. Then, a cutoff value of ‘0’ for the spatial-geneset-score was used to identify the final cell type. If the tumour calculation score of a spot is greater than ‘0’, the spot was defined as a ‘Tumour area’; if a spot had a FC or AFC calculation score greater than ‘0’ and a tumour calculation score less than ‘0’, the spot was defined as an ‘FC area’; and if the immune calculation score of a spot was greater than ‘0’, and the tumour & FC calculation score are less than ‘0’, the spot was defined as an ‘Immune area’; and if the calculation score of all cell in a spot are less than ‘0’ the spot was defined as an ‘NA area’ (Table S3).

Cell trajectory analysis

The monocle developmental trajectory was calculated on monocle 3 (the method used by SPATA2⁴⁶) or monocle 2 (monocle package, version 2.18.0).^7,46,47 Raw count served as input for Monocle, and pseudo-time was computed and graphed using the ‘reduceDimension (DDRTree)’, ‘orderCells’ and ‘plot_cell_trajectory’ functions. Trajectory patterns were determined by the ‘differentialGeneTest’ function.

Tumour foci heterogeneity measurement

To identify the heterogeneity of tumour foci, the Seurat object was transformed to SPATA data by the SPATA2 R package (version 4.0.2).^48,49 Tumour foci were the basic subunit of the tumour, with at least three gathered tumour spots. These foci were far away from other gathered tumour foci (at least 25 μm), and the other scattered spots were also treated as tumour foci for further analysis. Tumour foci were segmented in the SPATA2 environment, with the focus of the tumour foci heterogeneity on the common subpopulations among different patients. For this reason, we used the cell label transfer method (find anchors between a reference and query object, and use the anchors to transfer data from the reference to query object^50,51) to ensure that the Visium spatial samples shared the same subtypes of tumour area, FC area and immune area. We used patient 1 to explore the heterogeneity among different tumour foci because the tumour foci of patient 1 were highly heterogeneous. Then, the tumour foci were chosen to obtain differentially expressed genes, and the monocle developmental trajectory was calculated on monocle.^7,46,47 We then performed Gene Ontology (GO)/Kyoto Encyclopedia of Genes and Genomes pathway analysis and Gene Set Variation Analysis to identify the enriched key gene sets among the selected tumour foci, according to the enrichment score of gene sets from the Molecular Signatures Database. Survival analysis of signature genes on PTC were done on the Gepia2 database (based on the thyroid carcinoma data of The Cancer Genome Atlas [TCGA] Program).⁵²

External tumour microenvironment cell distribution measurement

The distances between cells and tumour foci were calculated to identify the cell distribution of the tumour microenvironment. The cell position information was extracted. The distance between selected cells and all tumour cells was calculated, and the shortest distance was used to classify cells surrounding the tumour foci (Get-nearby-cell function). One-way analysis of variance (ANOVA) was performed for the fractions of cell types surrounding all tumour foci to obtain the cell distribution in the external tumour microenvironment. Next, the surrounding cell features and the cluster interaction of tumour and microenvironment cells were identified through cluster trajectory analysis. In the SPATA2 environment, the cell trajectory from tumour to FCs passing by AFCs was drawn, and the biological features (genes or gene sets) were fitted according to the internal standard mathematical model (including linear ascending/descending, logarithmic ascending/descending, sinus and sinus (reversed)).^48,49

Contour map data processing

For each biological feature (genes or geneset-based score), we generated a contour map to show the value distribution on the spatial transcriptome plane.^53,54 We used the Python package Matplotlib to draw the contour map, which requires a 2-D numeric matrix as the input. The 10× Visium array arranges the data spots in a discontinuous style on the 2-D rectangular coordinate system. For a column with an odd/even X coordinate, there are only odd/even values for Y. Therefore, the actual data spots cover only half of the spots within the tissue region. For one biological feature, we filled the discontinuous spots within the tissue region with the mean value of 4 adjacent spots with actual data, and all the spots outside the tissue region were filled with the minimum value minus the standard deviation of the values of actual data spots. A contour map contains many contour line regions at several levels of the biological feature. To understand the connection between the level of biological features and the regions formed by one type of cell, we compared the overlap relationship between each contour line region and each cell type region. We defined a cell type region as the convex hull of a set of actual data spots with the same cell type. The spot set should have at least five spots with Chebyshev distances equal to or less than 2. This process was performed using the Python package NetworkX. We defined the match quality of a contour line region with a cell type region as the square of the overlap area divided by the product of the two region areas. For each cell type region, we sorted the overlapping contour line regions in descending order of the contour level and the match quality. The top genes of the region-matched results were considered specific to the cell type region.

Copy number variation analysis

The data extracted from SPT data were analysed with the infercnv package (version 1.3.3) for copy number variations (CNVs).^55,56 In the infercnv analysis, we configured the parameters by setting ‘denoise’, ‘default hidden Markov model (HMM)’ and ‘cluster_by_groups’ to TRUE. To analyse the gene CNV expression changes at various locations on the genome of tumour cell and AFC cells, the FCs were used as the reference cell, and its modified CNV expression value was defined as 1. Hence, the modified CNV level (relative CNV expression) of genes on each chromosome were summarised in the heatmap. These genes with modified CNV level greater than or less than 1 imply more gain or loss. And the top gain and loss genes were then extracted based on the modified CNV value and showed in the violin plots. On the other hand, we also evaluated and graphed the CNV values of some genes that are commonly genetic altered in PTC disease, including BRAF, RAS, AKT, PTEN, PIK3CA mutations, as well as PAX8–PPARG fusions. We calculated the CNV values of these genes and compared their values between AFCs and tumour cells. Gain of an oncogene may also represent an oncogenic stimulus, while loss can substitute for inactivating mutations of tumour suppressors.⁵⁷ Furthermore, the correlation of sample mutations in FCs, AFCs and tumour cells were calculated based on the relative expression of all genes.

Single-cell RNA sequencing data analysis

The data extracted from GSE184362, GSE158291 and GSE232237 were reanalysed with the Seurat package.^4,22 Specifically, cells from PTC patients were filtered based on the label in the meta data; cells with percent.mt < 40%, UMIs > 200 and UMIs < 6000 were filtered; and the top 2000 VariableFeatures were used for downstream data scale, ‘RunPCA’, ‘FindNeighbors’, unsupervised clustering and dimensional reduction analysis. The entire thyrocyte population (tumour cells and FCs) and immune populations (T, NK, B and myeloid) were isolated for further analysis. Then, we employed the cell label transfer function to identify anchors between a reference patient and query patients and then transferred the cell label data. Patient No. 1 in our work was selected as the reference patient, and the patients from scRNA-seq data (GSE184362, GSE158291 and GSE232237) were served as query patients, respectively. This analysis enabled the identification of both common and unique clusters of AFCs, FCs and tumour cells in the entire thyrocyte population of these scRNA-seq data (GSE184362, GSE158291 and GSE232237). Besides, the expression patterns of TGFB1 and CXCL13 in the scRNA-seq data (GSE184362) were analysed. The composition and percentage of clusters were further identified and calculated in R.

Multiplex immunohistochemistry

FFPE tissue sections were stained with the PANO 7-plex IHC kit (Panovue) according to the manufacturer's instructions, including deparaffinisation with xylene, hydration with gradient ethanol, antigen retrieval with sodium citrate buffer (0.01 M, pH = 6.0), removal of endogenous peroxidase with H₂O₂, blocking tissues with 10% goat serum and staining with antibodies and TSA-RM. The next cycle of staining started with antigen retrieval and was stopped with TSA labelling. The antibodies used included CD4 (ab133616, 1:500; Abcam) with Opal 620, CD8 (C8/144B, 1:200; CST) with Opal 690, and CD20 (ab78237, 1:2000; Abcam) with Opal 520. The antibodies used included CD45 (20103-1-AP, 1:2000; Proteintech) with Opal 570, SFRP4 (15328-1-AP, 1:200; Proteintech) with Opal 620, and IgG (EPR4421, 1:500; Abcam) with Opal 520. Then, all slides were scanned and analysed using the Vectra Automated Quantitative Pathology Imaging System (Vectra Polaris featuring MOTiF™). Finally, the percentage of SFRP4⁺ AFCs around the Tumour area was analysed in the HALO software (version 3.3). Specifically, the typical tumour area were identified based on DAPI staining, the distance and the area of 0–2, 2–4 and >4 mm from the tumour area were delineated, and the percentage of SFRP4⁺ AFCs in DAPI⁺ cells were calculated and statistical analysed.

Statistical analysis

All statistical data are presented as the mean ± standard deviation (SD). The normality test and homogeneity of variance test were conducted for unpaired samples (n ≥ 3) (GraphPad Prism, Version 9.0). For two sets of samples, if normality and variance homogeneity were met, a two-tailed unpaired Student's t-test was used; if not, the Mann–Whitney test was used. For multiple sets of samples, one-way ANOVA was used if the sets were normally distributed; otherwise, a nonparametric test was used. In the scRNA-seq data and spatial transcriptomes data, the difference of gene expression across the populations were calculated by wilcox.test function in R (Mann–Whitney’ test). A p value < .05 was considered statistically significant (*p < .05; **p < .01; ***p < .001).

RESULTS

Establishment of a spatial RNA-clinical analysis method for the PTC microenvironment

The tumour microenvironment of PTC consists of various cell types, including typical and atypical thyroid follicular cells (FCs and AFCs), tumour cells and immune cells. These cells aggregate within and around tumour foci, and their functional and distributional heterogeneity is greatly influenced by the disease state. Clinical pathologists usually accurately identify regions of FC, tumour cells and typical immune cells in intraoperative frozen PTC sections (Figure 1A, left). However, the presence of spatially atypical regions (Figure 1A, right) can lead to a certain rate of missed diagnoses, impeding the early detection and timely treatment of PTC patients.

[IMAGE OMITTED. SEE PDF]

In this study, we developed a novel spatial RNA-clinical combination analysis method to investigate the cellular heterogeneity of PTC in greater detail. By integrating SPT data with pathological diagnostic information, this method enables more accurate extraction of signature genes and delineation of regions containing typical and AFCs, tumour cells and immune cells within PTC. To overcome the pixel limitation associated with SPT-HE staining, we enhanced the HE staining method on a mirror slide of spatial slides (Figures 1B and S1A–S1D). The magnified HE view of the mirror slides exhibited improved clarity in identifying the ‘pure’ FC spot, tumour spot and immune spot compared with SPT-HE staining (Figures 1B and S1A–S1D). Concurrently, clinical pathologists identified the clinically ‘pure’ tumour cells, FCs, AFCs and immune cells (Figures 1C and S1E).

The SPT data from tissue sections belonging to four PTC patients were merged using CCA (Figure S2A). Subsequently, 14 common clusters were identified based on highly variable genes and top principal components, as depicted in the UMAP plot (Figure S2B) and spatial distribution plot (Figure S2C–S2D). The signature genes associated with these common clusters were extracted and utilised to categorise RNA clusters into five distinct cell types, namely, tumour, FC (normal follicular cells), FC2 (follicular cells with specific genes), Ant (adjacent noncancerous tissue, including AFCs, fibroblasts and blood vessel cells) and immune cells (Figures 1D and S2E). By intersecting the RNA clusters with the clinically determined ‘pure’ cell types, new RNA-clinical populations were established (Figures 1E and S2F). Cells that exhibited the same RNA clusters as the clinically ‘pure’ types were utilised to extract signature genes specific to tumour cells, immune cells, FCs (FC2) and AFCs (Figure 1F). Employing this methodology, we screened these cells across all four tissue sections (Figure S2G) and identified the top 20 signature genes for each cell type, which were defined as the spatial RNA-clinical signature genes (Table S2). The expression patterns of the top five signature genes among the 20 identified RNA-clinical signature genes exhibited distinct differences between tumour cells, FCs and immune cells, as illustrated in the dot plot of the merged samples from CCA analysis (Figure 1G). This observation underscores the accuracy and specificity of the spatial RNA-clinical combination analysis method in characterising the cellular heterogeneity within PTC.

Application of the RNA-clinical analysis method to the cell cluster pattern study of PTC

Using the spatial RNA-clinical analysis method developed above, we obtained the characteristic spatial RNA-clinical signature genes for tumour cells, FCs, AFCs and immune cells in PTC. These top 20 signature genes were then utilised to assess the spatial-geneset-score of each spot in the spatial transcriptome data using the previously described spatial grading system.⁴⁰ The resulting spatial feature plot depicted the regional distribution spatial-geneset-score of different cell types (tumour cells, FCs and immune cells) throughout the entire tissue section (Figure 2A). In addition, our RNA-clinical method was confirmed to be consistent with the results of the single-cell transcriptome in some studies²¹ based on the characteristic gene calculations, suggesting that our RNA-clinical method is accurate and specific (Figures 2A and S3A and B). The colour transition from blue to red indicates higher spatial-geneset-scores and the aggregation of specific cell types. The critical value for the cell spatial-geneset-score was defined as 0 (Table S3). Since spots in the SPT data are not directly interpreted as cells,⁵⁸ the rules were applied to accurately identify the regions of cell clusters (Table S3). If the spatial-geneset-score of a spot in the tumour area exceeded the critical value, it was marked as ‘T’; otherwise, it was marked as ‘nT’ (indicating no tumour). Similarly, spots in different cell areas were compared with the critical value and labelled ‘FC’/’nFC’, ‘AFC’/’nAFC’ and ‘I’/’nI’ to represent FCs, AFCs and immune cells, respectively. The final cell label depended on the presence of ‘T’ or a combination analysis result of ‘FC’, ‘AFC’ and ‘I’ (refer to Section 2 for details). Consequently, the FC area (including typical and AFCs), tumour area, immune area and nonassigned (NA) area in PTC tissue sections were comprehensively and accurately identified (Figure 2B).

[IMAGE OMITTED. SEE PDF]

By employing this method, we not only accurately identified the cell types defined by the pathologist (Figures 1A and 2A) but also screened all distribution regions of a particular cell type within the entire tissue section, facilitating a comprehensive analysis of the cellular heterogeneity among tumour cells, FCs, AFCs and immune cells. After defining the regions of tumour cells, FCs and immune cells, the overall distribution of these cells was examined (Figure 2B). The entire FC area consisted of the following cell clusters: AFC, AFCI, FC1, FC1I, FC2 and FC2I (Figure 2C). Monocle analysis revealed that immune cells and most typical FC cells were different from the tumour cells in the developmental state, while a portion of the FCs (predominantly AFCs) was closer to the tumour cells in the developmental state (Figures 2D and S3C). The simplified pseudotime plot further demonstrated that AFCs were the cell type closest to the tumour area, indicating that AFCs might serve as transition cells between FCs and tumour cells, potentially participating in the dedifferentiation of FCs into tumour cells (Figure S3D). Furthermore, 13 tumour foci were identified, grouped and distributed within the tumour area (Figure 2E). Among them, tumour foci No. 1 (orange box), tumour foci No. 2 (red box) and tumour foci No. 4 (blue box) were also identified by clinical pathologists (Figure 2E). Cell trajectory analysis further revealed that tumour foci No. 2 and No. 4 developed in different states, as determined by monocle analysis (Figure 2F). These findings provide insights into the cell heterogeneity of AFCs and different tumour foci in PTC. Moreover, the impact of immune cells on tumour heterogeneity will be further described below. To further investigate the characteristics of cell heterogeneity, we employed the cell label transfer function to identify anchors between a reference patient and query patients and then transferred the cell label data. Patient No. 1 was selected as the reference patient, and the other patients served as query patients. This analysis enabled the identification of tumour foci across all PTC tissue sections. The spatial dimplot revealed both common and unique clusters of AFCs, FCs, immune cells and tumour foci among the tissue sections of the four patients (Figure 2G), but little is known about the tissue location and characteristic genes of AFCs and tumour foci.^8–12

AFC represents a transition state from FCs to tumour cells

As shown above, the pseudotime plot revealed that AFCs exhibited closer proximity to tumour cells than FCs (Figure S3D). To gain a deeper understanding of the cell heterogeneity within AFCs and its relationship with FCs and tumour cells in PTC, we constructed UMAP plots and spatial-geneset-score plot to visualise the developmental trajectories of AFCs, FCs and tumour cells (Figure 3A). The UMAP plot illustrated a closer relationship between AFCs and tumour cells compared with FCs (Figure 3A). Moreover, the spatial-geneset-score plot of signature genes demonstrated that the expression pattern of certain AFCs resembled that of tumour cells (Figures 3B and S4A). Furthermore, we examined the gene expression levels of AFCs and tumour cells and observed similarities in the expression patterns of top five signature genes (Figure 3C). Notable genes displaying similar expression patterns in AFCs and tumour cells compared with FCs included TFF3 and TG (signature genes of FCs), IGHG4 and SFRP2 (signature genes of FCs) and FN1 and SFTPA1 (signature genes of tumour cells) (Figures 3C and D). Especially, the immunoglobin gene IGHG4 was highly expressed in the AFC area (Figure S4B). GO term enrichment analysis demonstrated a greater overlap in up-regulated GO pathways between tumour cells and AFCs than between tumour cells and FCs (Figure S4C). Similar trends were observed in the spatial-geneset-score analysis results of hallmark gene pathways in FCs, AFCs and tumour cells (Figure 3E). For instance, enrichment of the epithelial–mesenchymal transition (EMT) pathway was observed in AFCs and tumour cells but not in FCs, while the oxidative phosphorylation pathway exhibited higher enrichment in FCs than in AFCs and tumour cells (Figure 3F). Moreover, CNV⁵⁶ analysis of the key genes in PTC (such as KRAS, AKT1 and PIK3CA)⁵⁵ and the top gain and loss genes, and analysis of whole chromosome genes were conducted in FCs, AFCs and tumour cells (Figures 3G and S4D–S4F). The clustering results of cells (Figure S4D), the modified CNVs expression pattern of genes (Figure 3G, S4E) and the correlation analysis result of CNVs (Figure S4F) further supported the notion that AFCs represent a transition state that is more closely related to tumour cells. These findings indicate that AFCs represent a transition state from FCs to tumour cells.

[IMAGE OMITTED. SEE PDF]

Additionally, we evaluated the cellular composition of the external tumour microenvironment based on the distance from the tumour foci (Figures 3H, I and S4G, S4H). The analysis revealed a decrease in the number of AFCs and an increase in the number of FCs as the distance from the tumour foci increased in patient No. 1 (Figure 3I). The number of AFCs decreased across all patients (Figure S4H). We also conducted a multiplex immunohistochemical staining experiment on PTC tumour tissue and found that the number of AFCs decreased as the distance from the tumour foci increased (Figures 3J and S4I). Importantly, evidence from online single-cell RNA-seq data (GSE184362, GSE158291 and GSE232237)^4,22 also indicated the existence of AFCs and the similarity between tumour cells and AFCs compared with FCs (Figures 3K and S4J, S4K). Collectively, our work provides additional evidence highlighting the heterogeneity and transition state of AFCs from FCs to tumour cells.

The key ligand–receptors interactions between AFCs and tumour foci

To explore the effects of AFCs on tumour foci, we conducted an analysis of cell‒cell communication involving tumour cells, AFCs and FCs based on ligand‒receptor pathways in the PTC tumour microenvironment. Cell‒cell communication, particularly through ligand‒receptor pathways, plays a crucial role in modulating cellular functions and distribution within the tumour microenvironment. We observed that AFCs were enriched around tumour foci (Figures 3H–J). Initially, we manually drew a trajectory line from tumour cells and AFCs to FCs in the PTC SPT data using the SPATA2 package (Figures S5A–C). Ligand‒receptor pathway analysis revealed a decreasing intensity trend of the FN1–SDC4 and FN1–ITGA3 pathways along this trajectory (Figure S5C). However, manually drawn trajectory lines are limited by the contingency of trend line location and orientation and a one-sided perspective based on the trajectory direction rather than the entire tumour foci. To overcome these limitations, we introduced the concept of contour lines from mathematics and geography. Specifically, we projected the expression pattern of the FN1 feature plot into a ContourPlot, demonstrating high similarity and feasibility (Figure S5D). Subsequently, ligand‒receptor pairs were evaluated and projected into a contour plot using the contour plot analysis method. The top three ligand‒receptor pathways from tumour cells to AFCs to FCs were identified as FN1–SDC4, FN1–ITGA3 and LAMB3–ITGA2 (Figures S5E and F). Notably, high expression levels of FN1–SDC4, FN1–ITGA3 and LAMB3–ITGA2 were associated with a decreased relapse-free survival fraction of PTC based on gene expression profile analysis from TCGA database (Figure S5G). This suggests that FN1–SDC4, FN1–ITGA3 and LAMB3–ITGA2 may function as potential prognostic biomarkers for PTC. Overall, our analysis provides insights into the functions of AFCs and their role in cell‒cell communication within the PTC tumour microenvironment. The findings highlight the transition state of AFCs from FCs to tumour cells and their potential impact on the progression and heterogeneity of PTC. Additionally, the identified ligand‒receptor pathways involving AFCs may have implications for the prognosis and therapeutic targeting of PTC.

The heterogeneity of tumour foci in PTC

Besides the heterogeneity of FC area (including AFCs and FC), we also noted heterogeneity of tumour foci in the pseudotime plot analysis of tumour cells (Figure 2F) and clinical definition of tumour area (Figure 1A). The statistical analysis of all tumour foci in the four samples demonstrated that tumour foci Nos. 2, 3 and 4 were shared among all PTC patients, while tumour foci Nos. 1 and 10 were shared by two patients. Other tumour foci were observed in only one patient (Figure 4A). Therefore, our focus was primarily on tumour foci Nos. 2, 3 and 4. Among these, tumour foci No. 2 exhibited high expression of surfactant protein A2 (SFTPA2), surfactant protein A1 (SFTPA1), immunoglobulin lambda constant 7 (IGLC7), keratin 6A (KRT6A) and keratin 17 (KRT17) genes, whereas tumour foci No. 4 demonstrated high expression of cystatin SN (CST1), cartilage oligomeric matrix protein (COMP), collagen type 10 alpha 1 chain (COL10A1), matrix metallopeptidase 11 (MMP11) and periostin (POSTN) genes. Tumour foci No. 3 did not exhibit prominent expression of specific genes (Figure 4B).

[IMAGE OMITTED. SEE PDF]

GO enrichment analysis revealed that pathways such as ‘respiratory gaseous exchange by respiratory system’, ‘cornification’, ‘Toll-like receptor (TLR) signalling pathway’, ‘cell killing’ and ‘pattern recognition receptor signalling pathway’ were up-regulated in tumour cells from tumour foci No. 2 (Figure 4C). TLR signalling is known to promote tumour growth, while cornification represents a unique form of terminal differentiation and programmed cell death. However, in tumour cells from tumour foci No. 4, pathways related to ‘extracellular matrix (ECM) organisation’, ‘extracellular structure organisation’, ‘collagen fibril organisation’, ‘ossification’ and ‘cartilage development’ were up-regulated (Figure 4C). ECM effects and collagen formation play important roles in tumour cell metastasis. Spatial-geneset-score analysis of these pathways confirmed the enrichment of ‘kegg ECM_receptor_interaction’, ‘naba collagens’ and ‘reactome ECM_proteoglycans’ in tumour foci No. 4 (Figure 4D). Additionally, analysis using ContourPlot revealed that hallmark angiogenesis and EMT pathways were more enriched in tumour foci No. 4 than in No. 2, suggesting a potential facilitation of tumour cell metastasis (Figure 4E). Notably, tumour foci No. 2 exhibited a similar but different level of CNV across the entire genomic region compared with tumour foci No. 4 (Figures S6A–S6C). Moreover, high expression of the top seven genes from tumour foci No. 2 was associated with decreased relapse-free survival in patients (Figure 4F).

Tumour immune microenvironment affected tumour foci heterogeneity in PTC

The presence and activity of immune cells within the tumour microenvironment play a significant role in shaping tumour characteristics and tumour progression. Thus, we further investigated the impact of the tumour immune microenvironment on the heterogeneity of tumour foci, with a specific focus on tumour foci No. 2 and tumour foci No. 4. To explore the interplay between immune cells and tumour cells, we performed NicheNet analysis to identify specific ligand‒receptor interactions. Our results revealed a high interaction potential between TGFβ secreted from immune cells and TGFBR1 and TGFBR2 in tumour cells (Figure 5A). Notably, the response level of TGFβ signalling was lower in tumour foci No. 2 than in tumour foci No. 4 (Figure 5B and Table S4).⁵⁹ The TGFβ signalling pathways were positively correlated with the signature genes of tumour foci No. 4 only (Figure 5C). TGFβ signalling is known to have tumour-suppressive effects by inducing growth arrest and apoptosis. The differential activity of TGFβ signalling in different tumour foci suggests its involvement in the outcome of PTC, as tumour foci No. 2 were associated with decreased relapse-free survival. The results were further verified by the quantitative PCR results of TPC1 cells treated with or without TGFBR1 inhibitor (PF06952229 and AZ12601011) (Figure 5D) and cell proliferation results of TPC1 treated with TGFβ and TGFBR1 inhibitor PF and AZ (Figure 5E).

[IMAGE OMITTED. SEE PDF]

To determine the specific immune cell populations responsible for TGFβ and CXCL13 secretion within the tumour microenvironment, we analysed immune spots using multiple methods, including spatial-geneset-score analysis, CIBERSORT and MCPcounter analysis (Figures 5F and S7A). We identified five immune cell clusters, with clusters 1 and 3 showing high expression of B-cell-associated genes and cluster 3 showing high expression of T-cell-associated genes (Figures 5G and S7B and D). Spatial distribution analysis revealed that B-cell-enriched spots were primarily located at the centre, while T-cell-enriched spots were predominantly distributed in the periphery of the tumour (Figures 5H and S7E). Furthermore, we confirmed the presence of B cells, CD8⁺ T cells and CD4⁺ T cells (tertiary lymphoid structures [TLSs]) within typical immune areas using multiplex immunohistochemical staining in PTC tissue sections (Figures 5I and J). Additionally, we observed that TGFβ, without knowing its specific source, was distributed across all five immune cell clusters in the spatial transcriptome data (Figure 5K). We also focused on the tumour immune microenvironment surrounding the tumour foci besides TGFβ signalling. We examined the changes in features from immune cells to tumour cells. Interestingly, we observed a decrease in ligand‒receptor interactions, particularly in the CXCL13–CXCR5 and CXCL12–CXCR4 axes, along the trajectory direction (Figures S7E–S7H). These interactions are crucial for modulating the infiltration, activation and differentiation of immune cells in response to tumour antigens. Analysis of single-cell sequencing data indicated that myeloid cells were the primary source of TGFβ production (Figure S7I).

Taken together, these findings suggest that myeloid cell-derived TGFβ within the tumour immune microenvironment may contribute to the heterogeneity of tumour foci in PTC. The interaction between immune cells and tumour cells, particularly through ligand‒receptor pathways and the secretion of TGFβ, plays a crucial role in shaping the characteristics of different tumour foci and potentially influencing patient outcomes. Further understanding of the tumour immune microenvironment and its impact on tumour heterogeneity could have implications for the development of targeted therapies and personalised treatment strategies for PTC.

Validation of the cellular heterogeneity in the PTC microenvironment

In addition, we added new samples in the validation cohort with PTC tumour and paratumour tissue to further proved our findings, combining histopathology and RNA expression patterns according to protocols previously described (Figures 6A, B and S8A, S8B). We also found AFCs in the cohort samples, and the common tumour foci across samples were tumour foci Nos. 2, 3 and 4 (Figure 6C). There was a high correlation score between the experimental cohort and the validation cohort, suggesting that our analytical methods are suitable for more PTC samples (Figure 6D). And the immunoglobin gene IGHG4 was highly expressed in the AFCs area (Figure 6E). The multiplex immunohistochemistry (mIHC) results on the tumour tissue of PTC also showed that IgG was also expressed on other cells besides immune cells (CD45 as a marker), including AFCs (SFRP4 as a marker) (Figure 6F). The CNV analysis also showed higher similarity of AFCs and tumour cells (Figure 6G). And the survival related genes (CTSC and LUM) isolated from tumour foci were also different in tumour foci 2 and No. 4 (Figure 6H). The distribution of immune clusters also exited in the validation cohort (Figure S8C). These data further demonstrated the heterogeneity of tumour foci, and AFCs represent a transition state from FCs to tumour cells.

[IMAGE OMITTED. SEE PDF]

DISCUSSION

Alterations in cellular heterogeneity correlate with human PTC progression and prognosis. SPT is a powerful technology for studying cellular heterogeneity. It is generally recognised that SPT is predominantly used in combination with single-cell RNA sequencing (scRNA-seq). However, the high cost of scRNA-seq and the modification of the real cell state after cell suspension limit its widespread application. Moreover, the key clinical pathology information of PTC tumours provided by pathologists is usually overlooked. Pathologists are experts in distinguishing between different cells on tissues based on apparent cell morphology, differentiation status and tumour diameter. This information is considered a key determinant of clinical diagnosis. Here, we demonstrate the feasibility of a novel experimental and analytical approach to investigate the heterogeneity of PTC using SPT in combination with clinical pathologist identification. First, by combining SPT with mirror slice H&E staining, we obtained a higher resolution of the pathological structure and cell composition of PTC, enabling more accurate pathological assessments. Then, we identified RNA-clinical signature genes through the integration of RNA clustering and clinical information and employed these genes to analyse the heterogeneity of tumour foci and the transitional state of AFCs in the progression to tumour cells. Our results uncovered the prognosis-associated cellular heterogeneity of tumour cells, FCs and AFCs and immune cells in the PTC tumour microenvironment.

Generally, normal cells undergo a process to become transitional cells and then tumour cells. The presence of transitional cells during initial histological examination suggests the possibility of an early tumour, warranting further investigation through subsequent biopsies. Failure to identify these residual transitional cells could pose a risk. For example, low- and high-grade intraepithelial neoplasia were thought to be transitional cells from normal epithelium to invasive carcinoma in oesophageal squamous-cell carcinoma and were found to play a key role in promoting oesophageal cancer development.⁶⁰ In addition, recent studies proposed the term ‘AFCs’ or premalignant FCs in PTC, meaning a different state from normal FCs, but it did not elucidate the tissue localisation and characteristic genes of AFCs.^8–12 Here, by the spatial RNA-clinical analysis method, our findings demonstrated that AFCs represent a transitional state between FCs and tumour cells, as evidenced by developmental trajectory analysis,^7,47 spatial distribution patterns,^44,53 characteristic gene expression levels, gene enrichment pathways and CNVs.⁵⁵ In support of this, we also conducted a mIHC experiment and found that the number of AFCs decreased as the distance from the tumour foci increased. Our findings are consistent with scRNA-seq data from other studies that AFCs were a transitional state between FCs and tumour cells, exhibiting a higher resemblance to the latter.^4,22 To facilitate clinical translation, the specific diagnostic molecules associated with AFCs should be identified and validated in future studies. How AFC ultimately affects the outcome of human PTC is also worth investigating.

Recent insights have shown that the heterogeneity of tumour foci can significantly affect tumour prognosis. For example, distinct tumour foci exist in primary prostate cancers, and sequencing studies have shown that tumour foci heterogeneity poses obstacles for the diagnosis and treatment of prostate cancer.⁶¹ In contrast to some studies that used only scRNA-seq, our study focused on both the tumour cell level and the tissue and distribution heterogeneity of tumour foci. In our study, we successfully explored the heterogeneity of PTC tumour foci through signature gene extraction, pathway enrichment analysis, evolutionary state analysis and survival analysis.⁵² Three tumour foci were found to be shared among all patients out of the 13 observed, and these foci were focused on in the study of the relationship with PTC prognosis. Notably, tumour foci No. 2 displayed elevated expression levels of genes associated with lower relapse-free survival in PTC patients. Furthermore, our ContourPlot analysis model enabled the exploration of ligand‒receptor interactions within tumour cells, AFCs and FCs.⁵⁴ We identified several enriched pathways, including FN1–ITGA3, LAMB3–ITGA2 and FN1–SDC4, along the trajectory from tumour to AFC to fully transformed cells. Moreover, high expression of these ligand receptors was associated with decreased relapse-free survival in PTC patients. These interactions also exist in other cancers, such as pancreatic adenocarcinoma.⁶² However, here, their role in PTC has been proposed for the first time. Understanding how the tumour microenvironment shapes tumour cells and their distinct functional properties within different tumour foci, including ligand‒receptor interactions, distances from normal tissues, immune microenvironment and gene expression patterns, will be crucial for future investigations. Further investigation is warranted to elucidate the detailed mechanisms by which these ligand‒receptor interactions benefit tumour cells.

Multiple health outcomes are associated with the types and status of immune cells in the tumour microenvironment. For example, TLSs, including T and B cells within and outside the tumour, have different prognostic impacts on patients with intrahepatic cholangiocarcinoma.^63,64 In our study, reduced interaction between myeloid-derived TGFB1 and TGFBR1 in tumour foci No. 2 contributed to tumourigenesis and increased heterogeneity, which was support by the spatial-geneset-score and nichenet analysis, quantitate PCR and cell proliferation validation results. Additionally, we observed alterations in the CXCL13–CXCR5 and CXCL12–CXCR4 ligand‒receptor interactions between immune cells and tumour cells, which potentially play important roles in modulating the infiltration, activation and differentiation of immune cells in response to tumour antigens.^65–67 These findings reveal potential molecules and methods by which immune cells remodel tumour foci heterogeneity and thus affect the prognosis of PTC patients. Extensive clinical trials are needed to further verify the key effects of TGFB1 or CXCL13–CXCR5 and CXCL12–CXCR4 axes in PTC treatment.

One limitation of the 10× Visium SPT protocol is its low resolution, which can lead to suboptimal staining effects and affect sequencing quality. To overcome this limitation, we improved the imaging results of H&E staining on mirror slides, enabling pathologists to clearly define pathological areas. By combining the advantages of clinically identified area types with high-throughput gene expression data, we obtained pure RNA-clinical cells and signature genes, which were integrated with spatial transcriptome data. This approach allowed for spatial-geneset-score analysis, precise cell identification and comprehensive clinical analysis of the tumour microenvironment. Another limitation in this study is that pathologists did not find the typical cancer-associated fibroblasts (CAFs) area in patient No. 1, which is important in the PTC tumour microenvironment, so we did not further investigate the relationship between FCs, tumour cells and CAFs. The heterogeneity of thyrocytes in PTC is well documented.^4,6,7 Our investigation identified a distinct cluster of follicular cells, denoted as FC2, characterised by the expression of immunoglobulin genes, normal nuclei morphology identified by the clinical pathologist and typical FC gene expression (e.g., TFF3 and TPO). Different from FC2, the expression pattern of AFCs demonstrated greater similarity to tumour cells than FC cells. And it is particularly intriguing that the AFC spots highly express many immunoglobulin genes, including IGHG4. Recent findings have demonstrated that tumour cells and tissue cells can express immunoglobulin genes.^68–73 Therefore, as the transitional states of the tumour cells, AFC cells may also express high-level IGHG4-related IgG. So, we mainly focused on the immunoglobulin genes expression in AFCs. Both the mIHC data and the validation cohort data further proved that immune-related genes (IGHG4) were highly expressed in AFC cells, especially in tumour tissue compared with tumour adjacent tissue. However, the functional implications of immunoglobulin expression in AFCs within the context of PTC remain uncertain.

In summary, our study introduced a novel approach to examine the heterogeneity of PTC using SPT and mirror slice H&E staining. We demonstrated the utility of this approach in integrating pathological information with gene expression data, enabling comprehensive analyses of the tumour microenvironment and identification of clinically relevant features. Our findings elucidated the transitional state of AFCs, the role of immunoglobulin expression in PTC, and the importance of ligand‒receptor interactions in tumour progression and immune responses (Figure 7). By uncovering the heterogeneity of tumour foci, this study contributes to a more comprehensive understanding of PTC biology and provides valuable insights for future diagnostic and therapeutic strategies.

[IMAGE OMITTED. SEE PDF]

AUTHOR CONTRIBUTIONS

Conceptualisation: K. Y. and B. X. Methodology: K. Y., Q. Z. L. and R. R. H. Software: K. Y., Y. H. J. and Z. H. B. Validation: K. Y. and R. R. H. Formal analysis: K. Y. and Q. Z. L. Investigation: K. Y. Resources: K. Y. and S. J. L. Data curation: K. Y. Writing – original draft preparation: K. Y., R. R. H. and Y. H. Writing – review and editing: Z. X. L., Z. B. Z., L. L., Q. Z. L., Q. L. Z., H. X. G., R. R. H. and K. T. Visualisation: H. B. and Y. H. J. Supervision: Z. X. L., H. X. G. and B. X. Project administration: Z. X. L. and B. X. Funding acquisition: Z. X. L., H. X. G. and K. Y.

ACKNOWLEDGEMENTS

This work was supported by the Program for Guangdong Introducing Innovative and Entrepreneurial Teams (2017ZT07S054), National Natural Science Foundation of China (8197227 and 82173033), Key-Area Research and Development Program of Guangdong Province (2021B0101420005) and Guangzhou Municipal Science and Technology Bureau (2023A03J0967).

CONFLICT OF INTEREST STATEMENT

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

DATA AVAILABILITY STATEMENT

Raw sequencing data have been deposited in the Genome Sequence Archive (GSA-human: HRA003537) and are accessible at https://ngdc.cncb.ac.cn/gsa-human/s/O0X38Rcb. The major steps for processing the spatial transcriptome were integrated as an R package in GitHub (https://github.com/AlexyanKai/Spatial-pathology-pipeline). Supplementary Table 5 provides a summary of the source of the materials, instruments and analysis tools. All study data are included in the article and/or supplementary material.

ETHICS STATEMENT

All procedures and the use of tissue were performed according to recent national ethical guidelines. This study was conducted in accordance with the recommendations of the Guide for the Medical Ethics Committee of the Second Affiliated Hospital of the South China University of Technology (K-2019-185). Written informed consent was obtained from all participants included in the study.

References

McLeod DS, Sawka AM, Cooper DS. Controversies in primary treatment of low-risk papillary thyroid cancer. Lancet. 2013; 381 : 1046-1057.

Schneider DF, Chen H. New developments in the diagnosis and treatment of thyroid cancer. CA Cancer J Clin. 2013; 63 : 374-394.

Xing M, Haugen BR, Schlumberger M. Progress in molecular-based management of differentiated thyroid cancer. Lancet. 2013; 381 : 1058-1069.

Pu W, Shi X, Yu P, et al. Single-cell transcriptomic analysis of the tumor ecosystems underlying initiation and progression of papillary thyroid carcinoma. Nat Commun. 2021; 12 : 6058.

Hu A, Yang LY, Liang J, Lu D, Cao FF. Single-cell RNA sequencing reveals the regenerative potential of thyroid follicular epithelial cells in metastatic thyroid carcinoma. Biochem Biophys Res Commun. 2020; 531 : 552-558.

Pan J, Ye F, Yu C, et al. Papillary thyroid carcinoma landscape and its immunological link with hashimoto thyroiditis at single-cell resolution. Front Cell Dev Biol. 2021; 9 : [eLocator: 758339].

Luo H, Xia X, Kim GD, et al. Characterizing dedifferentiation of thyroid cancer by integrated analysis. Sci Adv. 2021; 7 .

Luu MH, Fischer AH, Stockl TJ, Pisharodi L, Owens CL. Atypical follicular cells with equivocal features of papillary thyroid carcinoma is not a low-risk cytologic diagnosis. Acta Cytol. 2011; 55 : 526-530.

VanderLaan PA, Marqusee E, Krane JF. Clinical outcome for atypia of undetermined significance in thyroid fine-needle aspirations: should repeated fna be the preferred initial approach? Am J Clin Pathol. 2011; 135 : 770-775.

Kim M, Kim JE, Kim HJ, Chung YR, Kwak Y, Park SY. Cytologic diagnosis of noninvasive follicular thyroid neoplasm with papillary-like nuclear features and its impact on the risk of malignancy in the bethesda system for reporting thyroid cytopathology: an institutional experience. J Pathol Transl Med. 2018; 52 : 171-178.

Ferraro V, Sgaramella LI, Di Meo G, et al. Current concepts in parathyroid carcinoma: a single centre experience. BMC Endocr Disord. 2019; 19 : 46.

Haddad RI, Bischoff L, Ball D, et al. Thyroid carcinoma, version 2.2022, NCCN clinical practice guidelines in oncology. J Natl Compr Canc Netw. 2022; 20 : 925-951.

Galdiero MR, Varricchi G, Marone G. The immune network in thyroid cancer. Oncoimmunology. 2016; 5 : [eLocator: e1168556].

Angell TE, Lechner MG, Jang JK, LoPresti JS, Epstein AL. MHC class I loss is a frequent mechanism of immune escape in papillary thyroid cancer that is reversed by interferon and selumetinib treatment in vitro. Clin Cancer Res. 2014; 20 : 6034-6044.

Hanahan D. Hallmarks of cancer: new dimensions. Cancer Discov. 2022; 12 : 31-46.

Schlumberger MJ. Papillary and follicular thyroid carcinoma. N Engl J Med. 1998; 338 : 297-306.

Coca-Pelaz A, Shah JP, Hernandez-Prera JC, et al. Papillary thyroid cancer-aggressive variants and impact on management: a narrative review. Adv Ther. 2020; 37 : 3112-3128.

Yan K, Da TT, Bian ZH, et al. Multi-omics analysis identifies FoxO1 as a regulator of macrophage function through metabolic reprogramming. Cell Death Dis. 2020; 11 : 800.

Liao J, Lu X, Shao X, Zhu L, Fan X. Uncovering an organ's molecular architecture at single-cell resolution by spatially resolved transcriptomics. Trends Biotechnol. 2021; 39 : 43-58.

Macerola E, Poma AM, Basolo F. NanoString in the screening of genetic abnormalities associated with thyroid cancer. Semin Cancer Biol. 2022; 79 : 132-140.

Yoo SK, Song YS, Lee EK, et al. Integrative analysis of genomic and transcriptomic characteristics associated with progression of aggressive thyroid cancer. Nat Commun. 2019; 10 : 2764.

Peng M, Wei G, Zhang Y, et al. Single-cell transcriptomic landscape reveals the differences in cell differentiation and immune microenvironment of papillary thyroid carcinoma between genders. Cell Biosci. 2021; 11 : 39.

Stahl PL, Salmen F, Vickovic S, et al. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics. Science. 2016; 353 : 78-82.

Maniatis S, Petrescu J, Phatnani H. Spatially resolved transcriptomics and its applications in cancer. Curr Opin Genet Dev. 2021; 66 : 70-77.

Bassiouni R, Gibbs LD, Craig DW, Carpten JD, McEachron TA. Applicability of spatial transcriptional profiling to cancer research. Mol Cell. 2021; 81 : 1631-1639.

Danaher P, Kim Y, Nelson B, et al. Advances in mixed cell deconvolution enable quantification of cell types in spatial transcriptomic data. Nat Commun. 2022; 13 : 385.

Asp M, Giacomello S, Larsson L, et al. A spatiotemporal organ-wide gene expression and cell atlas of the developing human heart. Cell. 2019; 179 : 1647-1660 e19.

Mantri M, Scuderi GJ, Abedini-Nassab R, et al. Spatiotemporal single-cell RNA sequencing of developing chicken hearts identifies interplay between cellular differentiation and morphogenesis. Nat Commun. 2021; 12 : 1771.

Fawkner-Corbett D, Antanaviciute A, Parikh K, et al. Spatiotemporal analysis of human intestinal development at single-cell resolution. Cell. 2021; 184 : 810-826 e23.

Ortiz C, Navarro JF, Jurek A, Martin A, Lundeberg J, Meletis K. Molecular atlas of the adult mouse brain. Sci Adv. 2020; 6 : [eLocator: eabb3446].

Joglekar A, Prjibelski A, Mahfouz A, et al. A spatially resolved brain region- and cell type-specific isoform atlas of the postnatal mouse brain. Nat Commun. 2021; 12 : 463.

Maynard KR, Collado-Torres L, Weber LM, et al. Transcriptome-scale spatial gene expression in the human dorsolateral prefrontal cortex. Nat Neurosci. 2021; 24 : 425-436.

Crosse EI, Gordon-Keylock S, Rybtsov S, et al. Multi-layered spatial transcriptomics identify secretory factors promoting human hematopoietic stem cell development. Cell Stem Cell. 2020; 27 : 822-839 e8.

Haas S, Baccin C, Al-Sabah J, Velten L, Lars S, Trumpp A. Combined single-cell and spatial transcriptomics to deconvolute the hematopoietic stem cell niche. Blood. 2018; 132 : 876.

Gao S, Shi Q, Zhang Y, et al. Identification of HSC/MPP expansion units in fetal liver by single-cell spatiotemporal transcriptomics. Cell Res. 2022; 32 : 38-53.

Ji AL, Rubin AJ, Thrane K, et al. Multimodal analysis of composition and spatial architecture in human squamous cell carcinoma. Cell. 2020; 182 : 497-514 e22.

Moncada R, Barkley D, Wagner F, et al. Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat Biotechnol. 2020; 38 : 333-342.

Nagasawa S, Kuze Y, Maeda I, et al. Genomic profiling reveals heterogeneous populations of ductal carcinoma in situ of the breast. Commun Biol. 2021; 4 : 438.

Andersson A, Larsson L, Stenbeck L, et al. Spatial deconvolution of HER2-positive breast cancer delineates tumor-associated cell type interactions. Nat Commun. 2021; 12 : 6012.

Wu R, Guo W, Qiu X, et al. Comprehensive analysis of spatial architecture in primary liver cancer. Sci Adv. 2021; 7 : [eLocator: eabg3750].

French JD, Bible K, Spitzweg C, Haugen BR, Ryder M. Leveraging the immune system to treat advanced thyroid cancers. Lancet Diabet Endocrinol. 2017; 5 : 469-481.

Cabanillas ME, McFadden DG, Durante C. Thyroid cancer. Lancet. 2016; 388 : 2783-2795.

Ahn S, Kim TH, Kim SW, et al. Comprehensive screening for PD-L1 expression in thyroid cancer. Endocr Relat Cancer. 2017; 24 : 97-106.

Longo SK, Guo MG, Ji AL, Khavari PA. Integrating single-cell and spatial transcriptomics to elucidate intercellular tissue dynamics. Nat Rev Genet. 2021; 22 : 627-644.

LiVolsi VA. Papillary thyroid carcinoma: an update. Mod Pathol. 2011; 24 (2): S1-S9. Suppl.

Kleino I, Frolovaite P, Suomi T, Elo LL. Computational solutions for spatial transcriptomics. Comput Struct Biotechnol J. 2022; 20 : 4870-4884.

Hong S, Xie Y, Cheng Z, et al. Distinct molecular subtypes of papillary thyroid carcinoma and gene signature with diagnostic capability. Oncogene. 2022; 41 : 5121-5132.

Ravi VM, Will P, Kueckelhaus J, et al. Spatially resolved multi-omics deciphers bidirectional tumor-host interdependence in glioblastoma. Cancer Cell. 2022; 40 : 639-655.

Ren Y, Huang Z, Zhou L, et al. Spatial transcriptomics reveals niche-specific enrichment and vulnerabilities of radial glial stem-like cells in malignant gliomas. Nat Commun. 2023; 14 : 1028.

Lotfollahi M, Naghipourfar M, Luecken MD, et al. Mapping single-cell data to reference atlases by transfer learning. Nat Biotechnol. 2022; 40 : 121-130.

Stuart T, Butler A, Hoffman P, et al. Comprehensive integration of single-cell data. Cell. 2019; 177 : 1888-1902 e21.

Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 2019; 47 : W556-W560.

Ligorio M, Sil S, Malagon-Lopez J, et al. Stromal microenvironment shapes the intratumoral architecture of pancreatic cancer. Cell. 2019; 178 : 160-175 e27.

Armingol E, Officer A, Harismendy O, Lewis NE. Deciphering cell-cell interactions and communication from gene expression. Nat Rev Genet. 2021; 22 : 71-88.

Morton LM, Karyadi DM, Stewart C, et al. Radiation-related genomic profile of papillary thyroid carcinoma after the Chernobyl accident. Science. 2021: 372.

Hovhannisyan G, Harutyunyan T, Aroutiounian R, Liehr T. DNA copy number variations as markers of mutagenic impact. Int J Mol Sci. 2019: 20.

Tanaka H, Watanabe T. Mechanisms underlying recurrent genomic amplification in human cancers. Trends Cancer. 2020; 6 : 462-477.

Ma Y, Zhou X. Spatially informed cell-type deconvolution for spatial transcriptomics. Nat Biotechnol. 2022; 40 : 1349-1359.

Verrecchia F, Chu ML, Mauviel A. Identification of novel TGF-beta /Smad gene targets in dermal fibroblasts using a combined cDNA microarray/promoter transactivation approach. J Biol Chem. 2001; 276 : 17058-17062.

Chen Y, Zhu S, Liu T, et al. Epithelial cells activate fibroblasts to promote esophageal cancer development. Cancer Cell. 2023; 41 : 903-918 e8.

Haffner MC, Zwart W, Roudier MP, et al. Genomic and phenotypic heterogeneity in prostate cancer. Nat Rev Urol. 2021; 18 : 79-92.

Hiroshima Y, Kasajima R, Kimura Y, et al. Novel targets identified by integrated cancer-stromal interactome analysis of pancreatic adenocarcinoma. Cancer Lett. 2020; 469 : 217-227.

Sautes-Fridman C, Petitprez F, Calderaro J, Fridman WH. Tertiary lymphoid structures in the era of cancer immunotherapy. Nat Rev Cancer. 2019; 19 : 307-325.

Ding GY, Ma JQ, Yun JP, et al. Distribution and density of tertiary lymphoid structures predict clinical outcome in intrahepatic cholangiocarcinoma. J Hepatol. 2022; 76 : 608-618.

Feng W, Huang W, Chen J, et al. CXCL12-mediated HOXB5 overexpression facilitates colorectal cancer metastasis through transactivating CXCR4 and ITGB3. Theranostics. 2021; 11 : 2612-2633.

Aoki T, Chong LC, Takata K, et al. Single-cell profiling reveals the importance of CXCL13/CXCR5 axis biology in lymphocyte-rich classic Hodgkin lymphoma. Proc Natl Acad Sci USA. 2021; 118 .

Cabrita R, Lauss M, Sanna A, et al. Tertiary lymphoid structures improve immunotherapy and survival in melanoma. Nature. 2020; 577 : 561-565.

Wang J, Lin D, Peng H, Huang Y, Huang J, Gu J. Cancer-derived immunoglobulin G promotes tumor cell growth and proliferation through inducing production of reactive oxygen species. Cell Death Dis. 2013; 4 : [eLocator: e945].

Wang Z, Geng Z, Shao W, et al. Cancer-derived sialylated IgG promotes tumor immune escape by binding to Siglecs on effector T cells. Cell Mol Immunol. 2020; 17 : 1148-1162.

Zhao J, Peng H, Gao J, et al. Current insights into the expression and functions of tumor-derived immunoglobulins. Cell Death Discov. 2021; 7 : 148.

Zhao Y, Liu Y, Chen Z, Korteweg C, Gu J. Immunoglobulin g (IgG) expression in human umbilical cord endothelial cells. J Histochem Cytochem. 2011; 59 : 474-488.

Shao W, Zhang C, Liu E, et al. Identification of liver epithelial cell-derived Ig expression in mu chain-deficient mice. Sci Rep. 2016; 6 : [eLocator: 23669].

Shao W, Hu F, Ma J, et al. Epithelial cells are a source of natural IgM that contribute to innate immune responses. Int J Biochem Cell Biol. 2016; 73 : 19-29.