Integrated analysis of tumor microenvironment features to establish a diagnostic model for papillary thyroid cancer using bulk and single-cell RNA sequencing technology

Article type: Research Article

Keywords: Papillary thyroid cancer, Biomarker, Diagnostic model, Single-cell RNA sequencing, Tumor microenvironment

Authors: Yizeng Wang, Wenbin Song, Yingxi Li, Zhaoyi Liu, Ke Zhao, Lanning Jia, Xiaoning Wang, Ruoyu Jiang, Yao Tian, Xianghui He

Affiliations: https://ror.org/003sav965grid.412645.00000 0004 1757 9434Department of General Surgery, Tianjin Medical University General Hospital, 154 Anshan Road, Heping District, Tianjin, 300052 People’s Republic of China; https://ror.org/02mh8wx89grid.265021.20000 0000 9792 1228Key Laboratory of Immune Microenvironment and Disease (Ministry of Education), Tianjin Medical University, Tianjin, 300070 People’s Republic of China

License: © The Author(s) 2023 CC BY 4.0 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Article links: DOI: 10.1007/s00432-023-05420-8  |  PubMed: 37733241  |  PMC: PMC10645658

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (8.7 MB)

Introduction

Thyroid cancer (THCA) is one of the most common endocrine tumors globally (Xia et al. ref. 2022; Siegel et al. ref. 2023). Papillary thyroid cancer (PTC)—a pathological type of THCA—accounts for approximately 90% of THCA patients (Fagin and Wells ref. 2016; Chen et al. ref. 2023). Ultrasound is currently the standard technique for evaluating the characteristics of thyroid nodules. Existing risk stratification systems classify the risk of malignant tumors according to the ultrasonic findings of thyroid nodules (Kobaly et al. ref. 2022). Although 10–15% of nodules eventually prove malignant, 5% of patients have distant metastases at diagnosis (Alexander and Cibas ref. 2022). About 15% of patients treated with surgery relapse during postoperative follow-up, leading to a poor prognosis (Ito et al. ref. 2018; Alexander and Cibas ref. 2022; Fallahi et al. ref. 2022). Even after postoperative radioactive iodine therapy, some patients develop refractory PTC. With a better understanding of the molecular progression of THCA, several drugs have been developed for inhibiting oncogenic kinases or signaling kinases (RET/PTC, BRAF (V600E)), like those related to vascular endothelial growth factor receptor platelet-derived growth factor receptor and platelet-derived growth factor receptor (Fallahi et al. ref. 2022). BRAF (V600E) mutations have been confirmed to be associated with poor prognosis of PTC patients (Costa et al. ref. 2008; Song et al. ref. 2021). Unfortunately, PTC patients without BRAF mutation could not benefit from these newly developed drugs. The effect of BRAF inhibitors in non-BRAF mutated cancers has reportedly been controversial (Agianian and Gavathiotis ref. 2018). Therefore, future research must delve into the molecular heterogeneity of PTC to find novel diagnostic biomarkers and provide patients with individualized treatment.

With the development of sequencing technology, researchers are able to observe the cellular level changes in tumor tissue though single-cell RNA sequencing (scRNA-seq) technology, which brings new strategies for tumor diagnosis and treatment (Hwang et al. ref. 2018). The tumor microenvironment, a daedal ecosystem comprise stromal, epithelial, and immune cells (Xiao and Yu ref. 2021). Different immune infiltrating state in the tumor microenvironment is related to distinct prognostic outcomes in tumor patients. For instance, the presence of tumor-associated fibroblast is correlated with worse outcomes in gastric, bladder, and breast cancers (Bartoschek et al. ref. 2018; Chen et al. ref. 2020; Li et al. ref. 2022). New markers and therapeutic targets have been developed by further exploring the tumor microenvironment. Sui et al. (ref. 2023) revealed the role of the CCL18/PITPNM3 ligand-receptor pair in the interaction between tumor cells and macrophages by single-cell data of esophageal squamous cell carcinoma, and the receptor could be used as a potential therapeutic target for esophageal cancer. Ma et al. (ref. 2020) found the heterogeneity of prostate tumor cells, and HPN could be used as an early-stage diagnostic marker for prostate cancer by analyzing single-cell data of prostate cancer. Therefore, a deep understanding of the tumor microenvironment of PTC helps to elucidate the underlying mechanisms of tumor progression and occurrence, as well as to search for new biomarkers and potential therapeutic targets.

The present study extensively characterized the tumor microenvironment of PTC and established a diagnostic model-based six genes. Our findings revealed extensive FN1/ITGB1 communication between the T/NK and other cells in PTC compared to normal thyroid tissue. Furthermore, we found dynamic changes in the evolution of thyroid cells into tumor cells and a high degree of heterogeneity among tumor cells. MDK/ALK/ALP and GAS/MERTK/AXL interactions between tumor cells and endothelial cells/fibroblasts may be potential therapeutic targets for PTC. Finally, we established a diagnostic model with good diagnostic efficacy for PTC by combining bulk RNA-seq. Our findings provided novel insights into the diagnosis of PTC.

Materials and methods

Data collection and data preprocessing

A total of 665 samples were included in the present study: 7 scRNA-seq samples (1 normal thyroid sample and 3 pairs of bilateral PTC) from the Gene Expression Omnibus (GEO) cohort (GSE191288, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE191288); 564 RNA-seq data (505 PTC and 59 normal samples) from The Cancer Genome Atlas (TCGA) cohort; 94 RNA-seq data (49 PTC and 45 normal samples) from GEO cohort. (GSE33630, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE33630). The gene expression profiles of the TCGA dataset were downloaded using the TCGAbiolinks package (v2.26.0), then converted to TPM format and standardized with log2.

The “Seurat” R package (v4.3.0) was used to conduct quality control procedures and downstream bioinformatics analyses for scRNA-seq datasets. We used the following criteria to obtain high-quality cells and filter out low-quality cells: the proportion of mitochondrial genes counts ≤ 15%; the proportion of erythrocyte genes counts ≤ 3%; UMIs ≥ 500; genes detected per cell ≥ 200. The DoubletFinder (McGinnis et al. ref. 2019) package was used to remove the doublets in each sample. After removing the doublets, 28,587 cells were included in further research. Subsequently, we normalized the scRNA-seq expression through the glmGamPoi (Ahlmann-Eltze and Huber ref. 2021) package. The top 3000 highly variable genes were identified using the SCTransform method and then used to calculate principal components. The TSNE and UMAP methods in Seurat were performed for cell clustering. The Clustree package was utilized to find the optimal cluster resolution. The markers used for cell identity were obtained from published research and the CellMarker database (Zhang et al. ref. 2019; Pu et al. ref. 2021).

Cell–cell interaction analysis

The CellChat (Jin et al. ref. 2021) package was used to evaluate the difference in cell–cell interaction between the normal and THCA samples. The two data sets were normalized separately, the high-variable genes were identified, and PCA analysis was performed in each of the two gene sets using the high-variable genes. Anchors were identified with the findinintegrationanchors() function, and the IntegrateData() function was used to combine two data sets. Additionally, the cell–cell interaction analysis among 16 cell types was also performed according to standard procedures. P value < 0.05 was considered statistically significant.

CNV analysis

SCEVAN—an R package—can infer non-malignant and malignant cells in the tumor microenvironment by calculating the raw count matrix of scRNA data. The thyrocyte cell expression matrix of six tumor samples (N1L, N1R, N2L, N2R, N3L, and N3R) were extracted and performed copy number variation (CNV) analysis to determine malignant thyrocyte cells with the SCEVAN package.

Single-cell downstream analysis

The SCP package (https://github.com/zhanghao-njmu/SCP) provides a comprehensive set of tools for single-cell data downstream analysis. The present study used SCP to perform multiple single-cell downstream analyses, such as pseudo-time trajectory analysis, enrichment analysis, GSEA analysis, and differential expression analysis. In detail, the trajectory analysis was performed by RunSlingshot() function in SCP to infer the evolution of thyrocytes into tumor cells. The tumor cells were extracted and regrouped into eight groups. Enrichment analysis between tumor groups was performed using the RunGSEA() function in SCP, and an adjusted P value < 0.05 was considered statistically significant. Furthermore, RunDEtest() function in SCP was used to identify the over-expression genes in eight distinct groups of tumor cells. Genes with log2FC > 1 and P value < 0.05 were considered over-expressed genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to annotate the biological processes over-expressed genes are involved in.

Clinical relevance

Differential expression analysis between thyrocytes and tumor cells was performed using the RunDEtest() function in SCP. We obtained 22 markers with log2FC > 1 and P value < 0.05. A total of 1186 genes up-regulated in PTC tissue were calculated by edgeR package. In detail, differential methods were set as ANOVA, and genes with log2FC > 1 and P value < 0.05 were considered up-regulated genes. The intersection of the aforementioned genes led to the identification of eleven genes. Boruta and LASSO method were utilized to identify the most important genes.

Logistic regression was used to construct a diagnostic model with important genes to better predict the thyrocyte tissue type. Receiver operating characteristic (ROC) curves were utilized to assess the discriminative performance. The calibration curves were applied to estimate the predictive accuracy of the model using the bootstrap method with 1000 re-samplings. To better explore model fitting, Hosmer–Lemeshow (HL) tests were performed, and P value > 0.05 was recognized as a good model fitting. The clinical applicability of the diagnostic model was estimated by the decision curve analysis (DCA). The TCGA dataset was used as training cohort, and GEO dataset was used as testing cohort.

RT-qPCR

Total RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) from frozen THCA tissue and corresponding control tissue resected surgically according to the manufacturer’s recommendations. RNA quality and concentration were measured using the NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). Real-time polymerase chain reaction (PCR) reactions were performed using SYBR Green PCR Master Mix (TransGen Biotech, Beijing, China) on the 7500 Real-Time PCR System (Applied Biosystems, Waltham, MA, USA). The specific sequences of primers are listed in Additional file 1: Table S1.

Statistical analyses

The RT-qPCR results were statistically analyzed using Prism 8 (Graph pad Software, CA) and are presented as mean ± standard deviation (SD) for at least three individual experiments. The statistical significance of differences was determined with the unpaired, two-tailed student t test, and P < 0.05 was considered statistically significant. Other statistical analyses in this study were performed using R studio software 4.2.2

Results

Identification of cell types

A total of 28,587 cells were obtained after quality control procedures through Seurat. All cells in the scRNA-seq dataset were clustered into 38 clusters using PCA with a resolution of 1. As shown in Fig. 1A, two-dimensionality reduction methods—TSNE and UMAP—revealed that different cell clusters were separated. We further annotated the cell types with cell markers from the Cell Marker database and published research. All cells were further annotated into eight cell types: T/NK cells, B cells, endothelial cells, myeloid cells, fibroblasts, mast cells, pericytes, and thyrocytes (Fig. 1B). The used cell markers are shown in Fig. 1C. For example, CD3D, CD3E, CD3G, and CD247 are markers of T/NK cells, and nine-cell clusters (2, 8, 13, 17, 20, 21, 23, 25, and 29) are identified as T/NK cells. These results indicated that the used markers could clearly distinguish different cell types. Furthermore, Fig. 1D shows the proportion of different cell types: thyrocytes (40.8%), T/NK cells (21.3%), pericytes (11.2%), myeloid cells (5%), endothelial cells (12.5%), fibroblasts (3.1%), B cells (4%), and mast cells (2%). Figure 1E depicts the proportion of cell types in normal and tumor samples. The proportion of non-epithelial cells in tumor samples is higher than that in normal samples, indicating a more complex microenvironment of tumor samples. These results revealed the characteristics of the data set and the differences in cell proportions between the samples. Taken together, we identified eight main cell types for further exploration.

Communication network difference between the normal and thyroid cancer tissue

Cell communication analysis can infer cell–cell interactions to further explore the changes in the tumor microenvironment. To analyze the changes in cell–cell interactions between normal and PTC tissues, we performed the cell communication analysis among seven cell types, excluding thyrocytes, using CellChat. Figure 2A depicts the number of interactions among seven cell types between normal (NT) and tumor samples, with cell–cell interactions in tumor tissues more than that in normal tissue. Subsequently, we compared the number and intensity of interactions between tumors and normal tissue (Fig. 2B). Compared to the normal sample, the number and intensity of interactions among endothelial cells, myeloid cells, and fibroblasts were higher in tumor tissues (N1L, N1R, N2L, N2R, N3L, N3R), suggesting the crucial role of three cell types in the formation and development of the PTC tumor microenvironment.

Considering the significance of T/NK cells in the TME, we further analyzed the interactions between T/NK and other cells, including themselves. For ligand-receptor pairs such as FN1–ITGB1, the probability of T/NK cells interacting with endothelial cells, fibroblasts, and myelocytes in PTC tissue was significantly higher than that in non-tumor tissue (Fig. 2C). This result suggested that FNI/ITGB1 could play a critical role in the tumor microenvironment of PTC. These findings reveal the potential role of extensive cellular communication in promoting tumor microenvironment formation in PTC tissue.

The cell trajectory of PTC epithelial cells

As not all epithelial cells in tumor tissues are malignant, SCEVAN was used to infer the benign and malignant nature of epithelial cells in tumor tissues. SCEVAN was used to analyze copy number variation (CNV) on 9577 thyrocytes from tumor tissues (N1L, N1R, N2L, N2R, N3L, and N3R). As shown in Fig. 3A, 5720 tumor cells were identified with abnormal levels of CNV compared to thyrocytes. Subsequently, we extracted 11,655 thyrocytes and tumor cells in 7 samples (1 NT and 6 tumor samples) for further analysis (Fig. 3B).

To investigate the dynamic development of epithelial cells in PTC microenvironment, we performed the cell trajectory analysis with the slingshot method in the SCP package to infer the trajectory from thyrocytes to tumor cells. As shown in Figs. 3C, D, thyrocytes were shown at the beginning of the trajectory, whereas tumor cells were located in the end of trajectory. To further explore the critical genes and biological progress in PTC progression, we analyzed the dynamic changes of genes in the trajectory “lineage 1”. As shown in Fig. 3E, the most important genes involved in PTC progression were identified: IGFBP7, MGP, PLVAP, IFI27, CA4, SRGN, A2M, RAMP3, RAMP4, TFF3, TPO, CRYAB, CYR61, SERPINA1, TIMP1, FN1, S100A6, TG, APOE, and APOC1. Additionally, the biological progress, including dysregulation of steroid metabolic process, dysregulation of thyroid hormone metabolic process, and positive regulation of fibroblast proliferation, were potentially involved in PTC progression. These findings reveal the dynamics of the biological processes that occur on the trajectory from thyrocytes to tumor cells. Furthermore, these genes closely related to tumor cell development may serve as potential biomarkers for PTC progression.

The heterogeneity between thyroid tumor cells

To further explore the heterogeneity among thyroid tumor cells, we re-clustered 5720 tumor cells using the same procedures and the SCTransform method. As shown in Fig. 4A, B, tumor cells were divided into eight groups with PCA. We performed gene set enrichment analysis (GSEA) based on each group compared to other groups. Compared with other groups, subgroup 0 of tumor cells was related to enzyme-linked receptor protein signaling pathway, organ growth, and DNA-templated transcription elongation (Fig. 4C). Subgroup 1 of tumor cells was mainly involved in tube morphogenesis, angiogenesis, blood vessel development, and morphogenesis, suggesting that this subgroup may be involved in angiogenesis in the tumor microenvironment. Subgroups 2 and 3 were associated with cytoplasmic translation, oxidative phosphorylation, and ribosome biogenesis compared with other groups (Fig. 4C). Similarly, Subgroups 4 and 5 were correlated with cell adhesion, cell motility, and cell migration (Fig. 4C), indicating the potential role of Subgroups 4 and 5 in tumor metastasis. In addition, Subgroups 6 and 7 were also involved in different biological processes (Fig. 4C). These findings indicate the characteristics of tumor subgroups at the level of the biological processes involved and the role of each subgroup in the TME.

We further analyzed the differential genes of each subpopulation compared to other subpopulations (Fig. 4D). The genes with logFC > 1 and P value < 0.05 were considered markers of each group (Table S2). As shown in Fig. 4E, the up-regulated genes of subgroup 0 were related to purine metabolism. The up-regulated genes of Subgroups 1 and 2 were involved in estrogen signaling pathways. Additionally, antigen processing and presentation were up-regulated in subgroups 2 and 5. Additionality, subgroup 5 was also associated with autoimmune thyroid disease. These findings reveal the heterogeneity of thyroid tumor cells, including intercellular metabolism, immunity, and tumor-related signaling pathways.

The crosstalk of tumor cells with other cells in the TME

To further explore the role of tumor cells in the TME, we analyzed all cell–cell interactions using CellChat. As shown in Fig. 5A, the interaction strength between tumor cell subgroups 0–7 and endothelial cells was highest, implying that the endothelial cells may be crucial in promoting tumor cell growth. Additionally, the interaction strength of fibroblasts on tumor cells was significantly enhanced compared with thyrocytes (Fig. 5A). Therefore, we further analyzed the ligand–receptor pair interaction between endothelial cells, thyrocytes, and tumor cells. Ligand–receptor pairs, such as PROS1-AXL, NAMPT-INSR, MDK-LRP1, MDK-ALK, GAS6-MERTK, and GAS6-AXL, between endothelial and tumor cells exhibited a higher communication possibility than that between endothelial cells and thyrocytes (Fig. 5B). These results revealed the potential interactions between tumor cells and endothelial cells in PTC.

Notably, MDK/LRP1, MDK/ALK, GAS6/MERTK, and GAS6/AXL between fibroblasts and tumor cells exhibited a higher communication possibility than that between fibroblasts and thyrocytes (Fig. 5C). In view of our previous findings in cell trajectory analysis that genes that positively regulate fibroblast growth were up-regulated in tumor cells, we analyzed the interactions of thyroid cells and tumor cells on fibroblasts. The probability of interaction between tumor cells and fibroblast receptor-ligand pairs was higher than that of thyroid cells, especially in TGFB1/ACVR1/TGFBR1, FGF18/FGFR1, BTC/EGFR, BMP8A/BMPR1A/BMPR2, and BMP8A/BMPR1A/BMPR2 (Fig. 5D). Consistent with the results of cell trajectory analysis, these results suggested that tumor cells could also have a regulatory effect on fibroblast growth via the aforementioned ligand-receptor pairs. Overall, these results revealed the potential role of endothelial cells and fibroblasts in regulating tumor cell growth in the tumor microenvironment.

Identification of critical diagnostic genes for PTC

To identify the critical diagnostic genes for PTC, we performed the differential expression gene analysis between thyrocytes and tumor cells in the scRNA-seq dataset, and 22 markers of tumor cells were identified (Table S3). Next, we performed the differential expression gene analysis between TCGA-THCA (PTC samples) and normal samples, and 1186 up-regulated genes in PTC were obtained (Table S4). After intersection of these two gene sets, we firstly identified eleven genes (TRPC5, TENM1, PDZRN4, NELL2, FRMD3, APOE, ARMCX3, DMD, APOC1, SLC35F3, AUTS2) highly expressed at both the tumor cell and tissue levels (Fig. 6A). Boruta package was used to identify importance of genes, and all eleven genes were confirmed as critical genes (Fig. 6B). We next performed LASSO method to identify six critical genes (TRPC5, TENM1, NELL2, DMD, SLC35F3, and AUTS2; Fig. 6C, D). As shown in Fig. 6E, these candidate diagnostic genes were obviously overexpressed in the tumor cells. Finally, the RT-qPCR analysis demonstrated that TRPC5, TENM1, NELL2, DMD, SLC35F3, and AUTS2 were highly expressed in PTC tissue than in normal tissue (Fig. 6F). In summary, TRPC5, TENM1, NELL2, DMD, SLC35F3, and AUTS2 could be potential biomarkers for PTC diagnosis.

Construction of a diagnostic model for PTC

To construct a diagnostic model for diagnosis for PTC, logistic regression was performed to construct a nomogram to predict the risk of PTC based on six critical genes (Fig. 7A). ROC R package was utilized to assess the discriminatory ability of the diagnostic model. The AUCs were 0.948 (Fig. 7B) and 0.965 (Fig. 7C) in the train and test sets, respectively, suggesting the effective role of the model to distinguish the PTC and normal tissues. In addition, the calibration curves in the train and test set suggested good predictive accuracy (Fig. 7D, E), with a P value for HL tests were 0.94 and 0.14 in train and test sets. Finally, the DCA were performed to analyze the clinical benefit of the model we constructed. As shown in Fig. 7F, G, the model we established could confirm the good benefits of clinical intervention. In a word, our model with good discriminatory ability could be served as a potential tool for PTC diagnosis.

Discussion

A large number of transcriptomic and genomic studies have been performed to identify diagnostic and prognostic biomarkers for THCA. Several effective diagnostic and prognostic biomarkers are prepared to apply in clinical practice for THCA, including somatic mutations and other molecular changes (Nikiforov and Nikiforova ref. 2011). However, PTC, as a main type of THCA, is a highly heterogeneous cancer. Previous studies have explored RNA expression in tumor tissue based on transcriptomics, which could not accurately reflect differences in gene expression between different cells within the tissue. Single-cell sequencing has emerged to identify individual-cell-level changes in gene expression. Therefore, with the progress of sequencing technology, we can better characterize the tumor microenvironment of PTC and find more specific tumor-related biomarkers for precise diagnosis and individual treatment.

Here, we comprehensively analyzed the tumor microenvironment of PTC. Our results revealed the potential aberrant interactions between the T/NK and other cells (endothelial cells, fibroblasts, and myeloid cells) through FN1/ITGB1. Fibronectin is confirmed to be a critical role of wound healing, cellular differentiation and growth, migration, and adhesion (Rick et al. ref. 2019). And integrin‐β (ITGB), a member of the integrin superfamily, is essential for cell proliferation, adhesion, and differentiation (Miranti and Brugge ref. 2002). Gu et al. (ref. 2023) found that overexpression of ITGB1 in gastric cancer was related to a poor prognosis and immune suppression. Also, increased ITGB1 expression was confirmed to be associated with poor prognosis and increased fibroblast infiltration in pancreatic ductal adenocarcinoma (Benesch et al. ref. 2022). Consistently with ITGB1, high FN1 expression was also related to M2 macrophage infiltration and poor prognosis in THCA (Geng et al. ref. 2021). These findings suggested the immune suppression role of FN1/ITGB1 in tumor microenvironment. Taken together, T/NK cells in PTC may form the immunosuppressive tumor microenvironment via FN1/ITGB1 interaction with other cells. Therefore, FN1/ITGB1 may be a potential therapeutic target to target the tumor microenvironment in PTC.

Through the CNV analysis, we identified thyroid tumor cells and analyzed the potential trajectories of thyroid cell differentiation into tumor cells. Several genes related to tumor cell differentiation were identified: CYR61, SERPINA1, TIMP1, FN1, S100A6, APOE, and APOC1. Studies have confirmed the prognostic values of CYR61 (Ren et al. ref. 2021), SERPINA1 (Wu et al. ref. 2021), and FN1 in THCA. By further subdividing the tumor cells, eight distinct tumor subgroups were identified. Tumor subgroups were involved in different biological processes, indicating a high degree of heterogeneity among tumor subgroups. The cell–cell interaction analysis revealed a strong interaction possibility between tumor cells and fibroblast/endothelial cells. Cancer-associated fibroblasts (CAFs) play essential roles in tumor development. They secrete extracellular matrix proteins, inflammatory ligands, and growth factors which promote cancer cell proliferation, immune exclusion, and therapy resistance (Biffi and Tuveson ref. 2021). Fibroblasts and endothelial cells are the primary sources of CAFs. Endothelial cells are involved in intravasation, which allows invasive cancer cells to translocate into the blood vessel lumen (Sobierajska et al. ref. 2020). These studies have revealed the critical role of endothelial cells and fibroblasts in TME. Further, our results identified potential ligand-receptor pairs (MDK/LRP1, MDK/ALK, GAS6/MERTK, and GAS6/AXL) for cell communication between these two types of cells and tumor cells. The MDK and GAS signaling pathway has been confirmed to regulate several biological processes in cells, including proliferation, survival and migration in tumor microenvironment by binding to their receptors (Wu et al. ref. 2018; Filippou et al. ref. 2020). Therefore, endothelial cells and fibroblasts may regulate the development of tumor cells and promote tumor microenvironment though these pairs.

To identify effective diagnostic markers, we intersected the up-regulated genes in tumor cells and up-regulated genes in PTC tissues. TRPC5, TENM1, NELL2, DMD, SLC35F3, and AUTS2 were identified as diagnostic biomarkers for THCA. To validate our findings, we performed RT-qPCR on 20 cases of paired PTC and para-cancer tissue. These gene expression levels confirmed the accuracy of our analysis. Although further study and clinical experiment are needed, the potential of the six genes to be successful diagnostic biomarkers for PTC has been consolidated in published studies. For example, Kechin et al. (ref. 2022) found that AUTS2 might be as one of important genes which can classify PTC in relation to the presence of driver NTRK-chimeric TRK genes. In addition, AUTS2 was also reported to be positively associated with other cancer progression via TGF-beta pathway activation, HEDGEHOG and WNT signaling pathway (Han et al. ref. 2015). These pathways play an essential role in tumor metabolism and immunity (Hanna and Shevde ref. 2016; Zhao et al. ref. 2020; Zou and Park ref. 2023). Cheng et al. (ref. 2017) reported that TEMN1 was highly expressed in cancerous thyroid tissues, and TENM1 expression in PTC was related to an advanced stage, BRAF V600E mutation, extra-thyroidal invasion and the classical subtype. These findings demonstrate the diagnostic efficacy and potential of candidate genes in differentiating benign and malignant thyroid tissue. In addition, the diagnostic model we constructed shows good differentiation and accuracy in both train and external test set, which indicates that the model is stable and practical. However, our study has some limitations. First, the scRNA-seq data in our study was limited. More single-cell data must be collected to reveal the prevalence of microenvironmental features in PTC. Additionally, in vitro and in vivo experiments are needed to further verify the mechanism of the receptor–ligand we screened for in PTC. Finally, we just analyzed and validated the mRNA levels of six candidate biomarkers; their application must be confirmed through further experiments and clinical studies.

In conclusion, we characterized the tumor microenvironment of thyroid papillary carcinoma and identified essential receptor and ligand pairs. Notably, we comprehensively explored the heterogeneity of thyroid tumor cells. Further, we established a diagnostic model based on six candidate diagnostic markers for PTC. Our findings provide new insights into the heterogeneity of PTC and the theoretical basis for diagnosing PTC.

References

1. B Agianian, E Gavathiotis. Current insights of BRAF inhibitors in cancer. J Med Chem, 2018. [PubMed]
2. C Ahlmann-Eltze, W Huber. glmGamPoi: fitting Gamma-Poisson generalized linear models on single cell count data. Bioinformatics, 2021. [PubMed]
3. EK Alexander, ES Cibas. Diagnosis of thyroid nodules. Lancet Diabetes Endocrinol, 2022. [PubMed]
4. M Bartoschek, N Oskolkov, M Bocci, J Lovrot, C Larsson, M Sommarin, CD Madsen, D Lindgren, G Pekar, G Karlsson, M Ringner, J Bergh, A Bjorklund, K Pietras. Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing. Nat Commun, 2018. [PubMed]
5. MG Benesch, R Wu, G Menon, K Takabe. High beta integrin expression is differentially associated with worsened pancreatic ductal adenocarcinoma outcomes. Am J Cancer Res, 2022. [PubMed]
6. G Biffi, DA Tuveson. Diversity and biology of cancer-associated fibroblasts. Physiol Rev, 2021. [PubMed]
7. Z Chen, L Zhou, L Liu, Y Hou, M Xiong, Y Yang, J Hu, K Chen. Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma. Nat Commun, 2020. [PubMed]
8. DW Chen, BHH Lang, DSA McLeod, K Newbold, MR Haymart. Thyroid cancer. Lancet, 2023. [PubMed]
9. SP Cheng, MJ Chen, MN Chien, CH Lin, JJ Lee, CL Liu. Overexpression of teneurin transmembrane protein 1 is a potential marker of disease progression in papillary thyroid carcinoma. Clin Exp Med, 2017. [PubMed]
10. AM Costa, A Herrero, MF Fresno, J Heymann, JA Alvarez, J Cameselle-Teijeiro, G Garcia-Rostan. BRAF mutation associated with other genetic events identifies a subset of aggressive papillary thyroid carcinoma. Clin Endocrinol (oxf), 2008. [PubMed]
11. JA Fagin, SA Wells. Biologic and clinical perspectives on thyroid cancer. N Engl J Med, 2016. [PubMed]
12. P Fallahi, SM Ferrari, MR Galdiero, G Varricchi, G Elia, F Ragusa, SR Paparo, S Benvenga, A Antonelli. Molecular targets of tyrosine kinase inhibitors in thyroid cancer. Semin Cancer Biol, 2022. [PubMed]
13. PS Filippou, GS Karagiannis, A Constantinidou. Midkine (MDK) growth factor: a key player in cancer progression and a promising therapeutic target. Oncogene, 2020. [PubMed]
14. QS Geng, T Huang, LF Li, ZB Shen, WH Xue, J Zhao. Over-expression and prognostic significance of FN1, correlating with immune infiltrates in thyroid cancer. Front Med (lausanne), 2021. [PubMed]
15. W Gu, H Sun, M Zhang, S Mo, C Tan, S Ni, Z Yang, Y Wang, W Sheng, L Wang. ITGB1 as a prognostic biomarker correlated with immune suppression in gastric cancer. Cancer Med, 2023. [PubMed]
16. Y Han, GQ Ru, X Mou, HJ Wang, Y Ma, XL He, Z Yan, D Huang. AUTS2 is a potential therapeutic target for pancreatic cancer patients with liver metastases. Med Hypotheses, 2015. [PubMed]
17. A Hanna, LA Shevde. Hedgehog signaling: modulation of cancer properies and tumor mircroenvironment. Mol Cancer, 2016. [PubMed]
18. B Hwang, JH Lee, D Bang. Single-cell RNA sequencing technologies and bioinformatics pipelines. Exp Mol Med, 2018. [PubMed]
19. Y Ito, A Miyauchi, M Kihara, M Fukushima, T Higashiyama, A Miya. Overall survival of papillary thyroid carcinoma patients: a single-institution long-term follow-up of 5897 patients. World J Surg, 2018. [PubMed]
20. S Jin, CF Guerrero-Juarez, L Zhang, I Chang, R Ramos, CH Kuan, P Myung, MV Plikus, Q Nie. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun, 2021. [PubMed]
21. AA Kechin, AA Ivanov, AE Kel, AS Kalmykov, IP Oskorbin, UA Boyarskikh, EA Kharpov, SY Bakharev, NA Oskina, OV Samuilenkova, IV Vikhlyanov, NE Kushlinskii, ML Filipenko. Prediction of EVT6-NTRK3-dependent papillary thyroid cancer using minor expression profile. Bull Exp Biol Med, 2022. [PubMed]
22. K Kobaly, CS Kim, SJ Mandel. Contemporary management of thyroid nodules. Annu Rev Med, 2022. [PubMed]
23. X Li, Z Sun, G Peng, Y Xiao, J Guo, B Wu, X Li, W Zhou, J Li, Z Li, C Bai, L Zhao, Q Han, RC Zhao, X Wang. Single-cell RNA sequencing reveals a pro-invasive cancer-associated fibroblast subgroup associated with poor clinical outcomes in patients with gastric cancer. Theranostics, 2022. [PubMed]
24. X Ma, J Guo, K Liu, L Chen, D Liu, S Dong, J Xia, Q Long, Y Yue, P Zhao, F Hu, Z Xiao, X Pan, K Xiao, Z Cheng, Z Ke, ZS Chen, C Zou. Identification of a distinct luminal subgroup diagnosing and stratifying early stage prostate cancer by tissue-based single-cell RNA sequencing. Mol Cancer, 2020. [PubMed]
25. CS McGinnis, LM Murrow, ZJ Gartner. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst, 2019. [PubMed]
26. CK Miranti, JS Brugge. Sensing the environment: a historical perspective on integrin signal transduction. Nat Cell Biol, 2002. [PubMed]
27. YE Nikiforov, MN Nikiforova. Molecular genetics and diagnosis of thyroid cancer. Nat Rev Endocrinol, 2011. [PubMed]
28. W Pu, X Shi, P Yu, M Zhang, Z Liu, L Tan, P Han, Y Wang, D Ji, H Gan, W Wei, Z Lu, N Qu, J Hu, X Hu, Z Luo, H Li, Q Ji, J Wang, X Zhang, YL Wang. Single-cell transcriptomic analysis of the tumor ecosystems underlying initiation and progression of papillary thyroid carcinoma. Nat Commun, 2021. [PubMed]
29. H Ren, X Liu, F Li, X He, N Zhao. Identification of a six gene prognosis signature for papillary thyroid cancer using multi-omics methods and bioinformatics analysis. Front Oncol, 2021. [PubMed]
30. JW Rick, A Chandra, C Dalle Ore, AT Nguyen, G Yagnik, MK Aghi. Fibronectin in malignancy: cancer-specific alterations, protumoral effects, and therapeutic implications. Semin Oncol, 2019. [PubMed]
31. RL Siegel, KD Miller, NS Wagle, A Jemal. Cancer statistics, 2023. CA Cancer J Clin, 2023. [PubMed]
32. K Sobierajska, WM Ciszewski, I Sacewicz-Hofman, J Niewiarowska. Endothelial cells in the tumor microenvironment. Adv Exp Med Biol, 2020. [PubMed]
33. E Song, M Jin, A Jang, MJ Jeon, DE Song, HJ Yoo, WB Kim, YK Shong, WG Kim. Mutation in genes encoding key functional groups additively increase mortality in patients with BRAF(V600E)-mutant advanced papillary thyroid carcinoma. Cancers (basel), 2021. [PubMed]
34. X Sui, C Chen, X Zhou, X Wen, C Shi, G Chen, J Liu, Z He, Y Yao, Y Li, Y Gao. Integrative analysis of bulk and single-cell gene expression profiles to identify tumor-associated macrophage-derived CCL18 as a therapeutic target of esophageal squamous cell carcinoma. J Exp Clin Cancer Res, 2023. [PubMed]
35. G Wu, Z Ma, Y Cheng, W Hu, C Deng, S Jiang, T Li, F Chen, Y Yang. Targeting Gas6/TAM in cancer cells and tumor microenvironment. Mol Cancer, 2018. [PubMed]
36. L Wu, Y Zhou, Y Guan, R Xiao, J Cai, W Chen, M Zheng, K Sun, C Chen, G Huang, X Zhang, Z Qian, S Shen. Seven genes associated with lymphatic metastasis in thyroid cancer that is linked to tumor immune cell infiltration. Front Oncol, 2021. [PubMed]
37. C Xia, X Dong, H Li, M Cao, D Sun, S He, F Yang, X Yan, S Zhang, N Li, W Chen. Cancer statistics in China and United States, 2022: profiles, trends, and determinants. Chin Med J (engl), 2022. [PubMed]
38. Y Xiao, D Yu. Tumor microenvironment as a therapeutic target in cancer. Pharmacol Ther, 2021. [PubMed]
39. X Zhang, Y Lan, J Xu, F Quan, E Zhao, C Deng, T Luo, L Xu, G Liao, M Yan, Y Ping, F Li, A Shi, J Bai, T Zhao, X Li, Y Xiao. Cell marker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res, 2019. [PubMed]
40. H Zhao, J Wei, J Sun. Roles of TGF-beta signaling pathway in tumor microenvirionment and cancer therapy. Int Immunopharmacol, 2020. [PubMed]
41. G Zou, JI Park. Wnt signaling in liver regeneration, disease, and cancer. Clin Mol Hepatol, 2023. [PubMed]