In the present study, the specimens from two female patients with PTC who underwent surgery at Jining First People's Hospital (Jining, China) and whose diagnosis was confirmed by histology were included. A metastatic PTC (PTC-M) tumor and its adjacent NT tissues were collected from a 44-year-old patient in June 2022, whose PTC metastasized to LNs in the III/IV/V/VI area of the left neck and VI regions of the right neck. Another tumor sample was taken from a 64-year-old patient with PTC without any metastasis in the LNs in the bilateral cervical area in June 2022. Both patients carried a BRAF gene mutation (V600E) that was revealed by PCR sequencing (data not shown) (11). The study was approved by the Ethics Committee of Jining First People's Hospital (Jining, China). All procedures carried out in this study were in line with ethical standards of the Helsinki Declaration (revised in 2013) and written informed consent was obtained from all participants. Although patient samples were collected after ethical approval was granted, the analysis performed was retrospective in nature, as the samples were collected after diagnosis and did not involve patient follow-up.

The fresh tissues were stored in sCelLive^® Tissue Preservation Solution (Singleron Biotechnologies) on ice after surgical removal within 30 min. The specimens were washed with Hank's Balanced Salt Solution three times, minced into small pieces and then digested with 3 ml sCelLive^® Tissue Dissociation Solution (Singleron Biotechnologies) using the PythoN^® Tissue Dissociation System (Singleron Biotechnologies) at 37°C for 15 min. The cell suspension was collected and filtered through a 40-micron sterile strainer. The mixture was then centrifuged at 300 × g and 4°C for 5 min, and following removal of the supernatant, the pellet was gently suspended with PBS. Finally, the samples were stained with 0.4% Trypan blue, mixed thoroughly at room temperature and immediately evaluated for cell viability under a microscope.

Single-cell suspensions (2×10⁵ cells/ml) in PBS were loaded onto a microwell chip using the Matrix^® Single Cell Processing System (Singleron Biotechnologies). The scRNA-seq libraries were constructed according to the protocol of the GEXSCOPE^® Single Cell RNA Library Kits (Singleron Biotechnologies) as previously described (12). Briefly, Barcoding Beads are collected from the microwell chip, followed by reverse transcription of the mRNA captured by the Barcoding Beads and to obtain cDNA, PCR amplification. The amplified cDNA was then fragmented and ligated with sequencing adapters. Individual libraries were diluted to 4 nM, pooled and sequenced on a Novaseq 6,000 (Illumina, Inc.) with 150 bp paired-end reads.

Raw reads from scRNA-seq were processed to generate gene expression matrixes using the CeleScope (https://github.com/singleron-RD/CeleScope) v1.9.0 pipeline. In brief, raw reads were first treated with CeleScope to remove low-quality reads with Cutadapt v.1.17 (13) to trim poly-A tail and adapter sequences. The cell barcode and Unique Molecular Indentifier (UMI) were extracted. Subsequently, STAR v2.6.1a (14) was used to map reads to the reference genome GRCh38 (Ensembl v92 annotation) (Index of/pub/release-92/fasta/homo_sapiens). The counts of UMI and genes in each cell were acquired with feature Counts v2.0.1 (15) software and used to generate expression matrix files for subsequent analysis.

Cells were filtered by gene counts <200 and the top 2% gene counts and the top 2% UMI counts. Cells with >30% mitochondrial content were removed. After filtering, 23,949 cells were retained for the downstream analyses, with on average 2,905 genes and 9,155 UMIs per cell. Functions from Seurat v3.1.2 (16) were used for dimension reduction and clustering. Subsequently, the NormalizeData and ScaleData functions were employed to normalize and scale all gene expression and the top 2000 variable genes were selected with the FindVariableFeautres function for principal component analysis. Using the top 20 principal components, cells were separated into multiple clusters with FindClusters. Cell clusters were visualized using t-Distributed Stochastic Neighbor Embedding (t-SNE) or Uniform Manifold Approximation and Projection (UMAP) with Seurat functions RunTSNE and RunUMAP.

To identify DEGs, the Seurat FindMarkers function was used based on the Wilcox likelihood-ratio test with default parameters and the genes expressed in >10% of the cells in a cluster and with an average log (fold change) value >0.25 were selected as DEGs. For the cell type annotation of each cluster, the expression of canonical markers found in the DEGs was combined with information from the literature and the expression of markers of each cell type was displayed using heatmaps/dot plots/violin plots that were generated with Seurat's DoHeatmap/DotPlot/Vlnplot function. Doublet cells were identified as expressing markers for different cell types and removed manually.

The cell type identity of each cluster was determined based on the expression of canonical markers found in the DEGs using the SynEcoSys^® database (https://singleron.bio). Heatmaps/dot plots/violin plots displaying the expression of markers used to identify each cell type were generated by Seurat v3.1.2 DoHeatmap/DotPlot/Vlnplot.

To investigate the potential functions of DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used with the ‘clusterProfiler’ R package v4.0.2 (17). Pathways with an adjusted P<0.05 were considered as significantly enriched. GO gene sets, including molecular function, biological process and cellular component categories, were used as references.

For Gene Set Variation Analysis (GSVA) pathway enrichment analysis, the average gene expression of each cell type was used as input data using the GSVA package (18).

The cell differentiation trajectory was reconstructed with Monocle2 (19). DEGs were used to sort cells in an order of spatial-temporal differentiation. The DDRTree function was used to perform FindVairableFeatures and dimension-reduction. Finally, the trajectory was visualized using the plot_cell_trajectory function. Next, CytoTRACE (v0.3.3) (20) (a computational method that predicts the differentiation state of cells from scRNA-seq data using gene counts and expression) was used to predict the differentiation potential of monocyte subpopulations.

The cell-cell interaction analysis was performed using CellPhoneDB v2.1.7 (21) based on known receptor-ligand interactions between two cell types/subtypes. Cluster labels of all cells were randomly permuted 1,000 times to calculate the null distribution of average ligand-receptor expression levels of the interacting clusters. The individual ligand or receptor expression was thresholded with a cut-off value based on the average log gene expression distribution for all genes across all of the cell types. The significant cell-cell interactions were defined as having P<0.05 and only receptors and ligands expressed in a threshold percentage of >0.1 of the cells in the specific cluster.

The Opal protocol staining approach was utilized to conduct mIHC staining, as previously described (22). mIHC was performed by several rounds of staining, each including a protein block with 1% BSA followed by addition of primary antibody and corresponding secondary horseradish peroxidase-conjugated antibody against mouse or rabbit immunoglobulins (ARH1001EA; Akoya Biosciences). The slides were then incubated in different Opal fluorophores (1:100) diluted in 1X Plus Amplification Diluent (ARD1001EA; Akoya Biosciences). After tyramide signal amplification and covalent linkage of the individual Opal fluorophores (NEL861001KT; Akoya Biosciences) to the relevant epitope or epitopes, the primary and secondary antibodies were removed via antigen retrieval, as previously mentioned, and the next cycle of immunostaining was initiated. The primary antibodies and Opal fluorophores were as follows: Lysosomal associated membrane protein 3 (LAMP3; cat. no. ab134045; diluted 1:100; Abcam) was labeled using Akoya Opal fluorophores 520. C-C motif chemokine receptor 4 (CCR4; cat. no. NBPI-86584; diluted 1:100; Novus Biologicals) was labeled using Akoya Opal fluorophores 620. CD4 (cat. no. ab133616; diluted 1:100; Abcam) was labeled using Akoya Opal fluorophores 690. The nucleus was labeled with DAPI (diluted 1:100; Akoya). The cover-slipping of all sections was performed using Anti-Fade Fluorescence Mounting Medium (cat. no. ab104135; Abcam). Multichannel imaging was performed on a PANNORAMIC SCAN II Imaging System (3Dhistech Kft.). Image analysis was performed in QuPath V.0.5.1 (Queen's University) (23).

Single-cell RNA-seq analyse was performed on both metastatic tumor and adjacent NT tissues from a PTC case with aggressive metastasis (PTC-M). Additionally, the tumor tissue collected from a patient with non-metastatic PTC served as a control for comparison with the metastatic tumor. The clinical information of the patients was collected at the time of recruitment and was described in the methods section.

After stringent quality filtering, a total of 23,401 single cells were included in the analysis. Transcriptional data from all of the cells were integrated and low-resolution t-distributed stochastic neighbor embedding clustering was used to identify eight major cell populations. These populations were labeled as T cells (CD3D, CD3E, CD2, TRAC, TRBC1 and TRBC2), B cells (MS4A1, CD79A and CD79B), plasma cells (JCHAIN, MZB1 and IGHG1), mononuclear phagocytes (MPs) (LYZ, CD14, C1QC, MRC1, CD68, CD163, CD1C and LAMP3), endothelial cells (PECAM1, VWF, CLDN5 and CDH5), mural cells (ACTA2, TAGLN, MYLK, MYL9, PDGFRB and NOTCH3), plasmacytoid dendritic cells (pDCs) (IL3RA, CLEC4C, LILRB4 and LILRA4) and epithelial cells (EPCAM, TG, DIO2, TACSTD2 and KRT19) (Fig. 1). A notable infiltration of MPs was observed in the PTC-M sample, where MPs accounted for 18.97% of total cells (998 cells), compared to only 1.92% in NT (148 cells) and 4.54% in PTC tissues (476 cells) (Fig. 1D). Fig. 1B displays the most distinct DEGs for each cell population.

MPs in PTC displayed increased expression of inflammation-related genes, including immune proteins (LYZ) and inflammatory cytokines (CCL4, CCL3, CXCL8, CCL4L2 and IL1B), as well as human leukocyte antigens (HLA) (HLA-DRA, HLA-DPA1 and HLA-DRB1) (Fig. 1B). DEG analysis revealed that the MPs in the PTC-M tumor tissue expressed a higher level of immune-suppression-related genes, such as MRC1, CD163, CCL18, SPP1, CCL2, CCR2, CD40, CCL17 and CCL22, whereas exhibiting lower levels of TNF, IL1B and RGS1 (Fig. S1A) when compared to the MPs from NT. The GO and KEGG analyses showed that the upregulated genes in MPs from PTC-M were mostly enriched in pathways related to the immune-inflammatory response, including ‘neutrophil activation’, ‘neutrophil degranulation’, ‘Staphylococcus aureus infection’ and ‘tuberculosis’ (Fig. S1B and C).

By conducting UMAP analysis, four distinct MP subgroups infiltrating the tumors were identified: Monocytes (LYZ, CD14, FCN1, VCAN and FCGR3A), macrophages (LYZ, MRC1, CD68, CD163, C1QC, APOE and MARCO), cDC2s (CD1C, CD1E and FCER1A) and mature DCs (LYZ, LAMP3 and CCR7) (Fig. 2A). Further analysis demonstrated that cDC2s and macrophages made up 42.67 and 40.45% of MPs in the PTC-M sample, respectively. In the NT sample, monocytes were the most common MPs (52.89% of MPs; Fig. 2C and D). Although a similar percentage of macrophage cells and cDC2s was found in PTC-M and PTC, the number of both subtypes of MP cell in PTC-M was ~2-fold higher than that in PTC (Fig. 2D). The top DEGs of each MP subtype are listed in Fig. 2B.

Macrophages were then annotated using UMAP and they were classified into 7 different clusters (Fig. 3A). Of note, the macrophage c3 cluster was solely enriched in PTC and PTC-M tumor tissues (Fig. 3A and B), demonstrating the unique tumor-associated function of immune-regulating macrophages. In addition, the macrophage c3 and c5 clusters expressed high levels of MRC1, CD163 and TGFB and low levels of IL1B, TNF and RGS1, which correspond to the alternatively-activated M2 subtype (Fig. 3C and D and Table SI) (24–26). It was also observed that CCL18, SPP1, GPNMB, CCL2, CCL8 and C1QA were upregulated in the macrophage c3 cluster (MRC1⁺CCL18⁺ M2), while in the macrophage c5 cluster, the expression of C1QB, C1QC and SLC40A1 was elevated (Fig. 3C; Table SI).

To analyze the TMEs potentially mediated by unique macrophage subtypes, DEG analysis was performed between the macrophage c1 and c3 clusters. Fig. 3D shows that in the macrophage c1 cluster, IL1B, TNF and RGS1, which are markers of classically-activated M1-type macrophages, were highly expressed (26), whereas the macrophage c3 subtype has a broad spectrum of molecules that are critical for pro-M2 conversions, such as GPNMB, or for inflammation, such as CTSL, APOE and PLA2G7 (27–30). Furthermore, GO annotation analysis revealed ‘positive regulation of cytokine production’, ‘leucocyte cell-cell adhesion’ and ‘T-cell activation’ in c1 macrophages, while ‘neutrophil activation’ and ‘neutrophil-mediated immunity’ were more likely activated in c3 macrophages. In addition, KEGG enrichment analysis showed that c3 macrophage-expressed genes were enriched in the ‘lysosome’, ‘phagosome’, ‘oxidative phosphorylation’ and ‘cholesterol metabolism pathways’ (Fig. S2). Expression of classical macrophage M1 and M2 marker genes is shown in Fig. 3E. Since c1 was annotated as M1-like, while c3 and c5 were denoted as M2-like macrophages, the data revealed that the ratio of M2-like/M1-like macrophages substantially increased in the progressively advanced stage (38/77 in PTC vs. 113/56 in PTC-M) of PTC (Fig. 3B).

UMAP analysis was used to sub-classify the monocytes into 8 clusters. However, it was not possible to tell the difference between functional differences or differences in the make-up of groups in the PTC, PTC-M and NT samples. A cell trajectory analysis was then performed to further investigate the potential transition between cell types. The pseudotime trajectory axis derived from Monocle indicates that monocytes could transdifferentiate into macrophage clusters, particularly c1 and c3 macrophages. The results illustrated that most monocytes remained quiescent in the NT sample, while they preferentially differentiated into the M1-like c1 macrophages in the PTC sample and more cells were directed towards the M2-like c3 macrophages in the PTC-M sample (Fig. 4A). Pseudotemporal expression patterns of representative genes also indicated the transition of monocytes into various macrophage clusters (Fig. 4B). These findings clearly outline the potential paths of differentiation for monocytes at a signal-cell level and suggest that they have the ability to give rise to a metastatic environment.

DCs have a crucial function in presenting antigens, initiating T-cell responses and migrating to lymph nodes (31). By applying a UMAP analysis, 6 clusters of cDC2s were identified (Fig. 5A and C). The cDC2 c1 cluster was the predominant cluster in the NT tissue and exhibited significant expression of IL1B, CCL3, CCL4, ATF3 and CXCL8 (Fig. 5B). The cDC2 clusters c2, c4 and c6 showed significant enrichment in the PTC-M sample. The cDC2 c4 cluster was specifically studied because it shared similarities with the subtype of ‘mature DCs enriched in immunoregulatory molecules’ (mregDCs) that have the ability to migrate towards lymph nodes (32–34). The mregDCs-like cluster exhibited simultaneous expression of type 2 T-helper cell responsive genes (CCL22 and CCL17) and other immune-regulatory genes (C1QA and C1QC), along with the DC maturation genes LAMP3, CD40 and CD80 (Fig. 5B; Table SII). By utilizing the feature plot and violin plot, it was possible to identify that the cDC2 c4 cells formed a distinct cluster characterized by elevated expression of immune-suppressive genes LAMP3 and CCL22 (LAMP3⁺CCL22⁺ DC) (Fig. 5D). Based on these observations, it may be proposed that the cDC2 c4 cluster likely functions as mregDCs, contributing to migration and the establishment of an immune-suppressive microenvironment in PTC-M.

To elucidate the interactions between different immune cell subtypes in the PTC microenvironment, a CellPhoneDB analysis was performed, focusing on ligand-receptor (L-R) interactions among macrophages, DCs and T cells (21). Based on the prior analysis, it was found that the majority of the projected crosstalk took place between macrophages c3, c5, cDCs c2 and T cells (Fig. 6A and B). Treg markers, such as TIGIT and CTLA-4, were expected to facilitate various immune-suppressive interactions between macrophages and DCs (35). In addition, it was projected that macrophage c3 and cDCs c2 could interact with Treg cells through the CCL18-CCR8 complex and CCL22-CCR4 complex, respectively (Fig. 6A) (36). The presence of cell-cell contacts was predominantly observed in the PTC-M sample, while being rare in the NT sample (Fig. 6B and C). It is worth mentioning that there are likely communications between cDCs c2 and macrophages in the tumor-associated environment and these communications are projected to be facilitated by the CXCL2/3/8-CXCR2 complex (Fig. 6A; Table SIII). In summary, the present data showed that the cDCs c2 subtypes had the strongest association with CD4⁺ T cells in terms of L-R pair and were possibly linked to the infiltration or dysfunction of T-cell subtypes.

In order to comprehend the function of T cells in the microenvironment, UMAP was used to categorize T cells into four distinct subclasses, including CD4⁺ naive T cells (CD4, CCR7, LEF1, SELL, TCF7 and IL7R), Treg cells (CD3D, FOXP3, IL2RA, CTLA4 and IKZF2), CD8+ effector T cells (CD3D, CD8A/B, NKG7, GZMA and GNLY) and proliferating T cells (TOP2A, MKI67 and CD3D) (Fig. 7A and C). Unsurprisingly, PTC-M tissue exhibited over four times the amount of T-cell infiltration and a greater proportion of Tregs compared to NT and PTC tissues, indicating a more immune-suppressive milieu in PTC-M (Fig. 7C). Fig. 7B displays the manifestation of distinct T-cell marker genes. Based on the examination of cell-cell communication, it was observed that the majority of interactions took place between the cDC c2 cluster and T cells (Fig. 6B; Table SIII). Upon further examination, the application of mIHC staining on PTC-M tumor sections allowed for the clear observation of the close proximity of LAMP3+ mregDCs and CCR4+ T cells (Fig. 7D). The ligand-receptor pair showed a specific enrichment in the PTC-M tissue, indicating its connection with the development of tumors and the metastasis of cancer in this particular malignancy.