All data are available to the public and were primarily sourced from The Cancer Genome Atlas (TCGA, http://portal,gdc.cancer.gov/). TC whole genome-wide expression profiles in transcripts per kilobase per million format, clinical annotations and simple nucleotide variation estimated by VarScan2 Variant Aggregation and Masking tool were retrospectively downloaded utilizing the TCGAbiolinks (version 2.25.0) package (20) from the TCGA database. Samples from TCGA-thyroid carcinoma (THCA; n=572) were included, with 513 TC samples and 59 healthy control samples. The current research respected the data access guidelines of every database. The entire research workflow, including data collection, processing, and analysis, is outlined in Figure 1.

The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo) contains abundant data from scRNA-seq. A GEO dataset (accession no. GSE184362; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE184362), consisting of 7 TC tissue samples was downloaded using scRNA-seq. Through the Seurat package from the R (version 4.2.0; http://www.R–project.org/) (21), the raw data from the dataset were imported as single-cell data. Initially, cells and genes were filtered out based on the following criteria: i) Genes expressed in <1 cell type were removed; ii) cells were deleted where <200 genes were expressed; iii) cells with gene expression between 200–3,000, mitochondrial gene content <20% and unique molecular identifiers of 500 and 12,000 were kept. Data were normalized using the normalizedata function (R Seurat package; version 5.2.1; http://satijalab.org/seurat). Differentially expressed genes (DEGs) in single cells were pinpointed through the connection between average expression and dispersion. Next, principal components analysis was conducted for graph-based clustering. The harmony method was used to eliminate the batch effect of various samples. The FindClusters function (R Seurat package; version 5.2.1) was adopted for shared nearest neighbor modularity optimization according to the clustering algorithm on 22 principal components with 0.2 resolutions, leading to 13 clusters. The RunUMAP function was adopted for Uniform Manifold Approximation and Projection (UMAP). UMAP-1 and UMAP-2 were used to showcase cell groupings. The FindAllMakers function (R Seurat package; version 5.2.1) was used to detect DEGs within cell clusters. Subsequently, cell clusters were distinguished using specific biomarkers for each cell type and the ratio of each cell type was assessed.

The AUCell package (version 1.18.0; http://bioconductor.org/packages/release/bioc/html/AUCell.html) (22) was adopted to score pathways for all cells based on gene set enrichment analysis (GSEA). According to the area under the curve (AUC) values of 134 genes related to glutamine metabolism from MSigDB (Table SI), gene expression was ranked for each cell type to determine the proportion of upregulated DEGs in all cells. Higher AUC values were observed in cells that expressed several genes from the gene set. The AUCell_exploreThresholds function (R AUCell package; version 1.18.0) was used to establish the threshold for detecting active cells. Next, the UMAP embedding was used to visualize the active clusters by mapping the AUC value of all cells with the ggplot2 package (version 3.3.5 http://ggplot2.tidyverse.org).

Pseudotime assay was performed based on Monocle 2 (23). To conduct the pseudotime assay on epithelial cells, the initial count data were adjusted by computing size elements for trajectory inference. Genes that were both widely spread out and highly active (empirical dispersion/dispersion fit ≥1 and mean expression ≥0.1) were used to construct the pseudotime trajectory (24). The DDRTree (version 0.1.5; http://CRAN.R-project.org/package=DDRTree) algorithm parameters were initialized with default values. To further examine these diverging occurrences, branched expression analysis modeling was employed within the framework of Monocle2 (version 2; http://cole-trapnell-lab.github.io/monocle-release) to discern all genes with branch-dependent expression (23). A heatmap was plotted using Monocle2 to visualize the branch-dependent expression patterns.

A cellular communications assay was adopted to evaluate the levels of ligand-receptor pairs in cells and exposed definite pathways (25). CellChat (version 1.1.3; http://github.com/sqjin/CellChat) assay was used to uncover the incoming and outgoing communication patterns of various cell types, measure different pathways of cellular communication and determine the flow of information through signaling pathways or the likelihood of cellular communication (26). CellChat was used to examine the single-cell specimens. Cellular communication in all TC samples was analyzed using CellChat; the processed scRNA-seq information from the Seurat package was imported into CellChat for analysis. The communication signals within TC cells were analyzed to determine the strength of all pathways and identify definite pathways. Default parameters were used in CellChat analyses. P<0.05 was considered to indicate a statistically significant difference and the P-value was adjusted using the Benjamini-Hochberg mean.

GO (27) analysis was performed to reveal the enriched biological process, molecular function and cellular component. KEGG (28) is a valuable bioinformatics tool used to discern metabolic pathways that are markedly related to the DEGs. The clusterProfiler (version 4.2.2) package (29) was used for GO and KEGG enrichment analyses (P<0.05) on glutamine metabolism-related DEGs in TC.

Prognostic genes were screened among DEGs using Kaplan-Meier (K-M) curves. The genes, meeting P<0.05, were selected for additional examinations. Tumor tissues with clinical information were assigned to a training set (n=348) and a validation set (n=154). Boruta algorithm (R package Boruta; version 8.0.0; http://gitlab.com/mbq/Boruta/) (30) and multivariate Cox regression were then adopted to reduce candidates and develop a risk model. According to the minimum criteria, the penalty parameter (λ) was determined. The risk score was reckoned as follows:

riskscore = Σⁿ_i=1 Coef(gene_i) * Expression(gene_i)

(Coef (genei): coefficients, Expression (genei): gene level).

The training cohort was classified into low- and high-risk subgroups. The K-M technique was used to plot survival curves and log-rank tests were conducted to check statistical significance. Receiver operating characteristic (ROC) curves were used to appraise the effectiveness of the prediction model. AUC >0.6 indicated satisfactory prediction performance. During the validation process, the validation cohorts were further categorized into subgroups based on their risk levels and the model accuracy was compared.

The clinical information (sex, age and tumor stage) of individuals was acquired from the TCGA database. Univariate Cox regression models were employed to analyze these features in conjunction with the risk score. Additionally, a nomogram was developed to forecast OS at 1, 3 and 5 years, with the risk score as a key prognostic indicator. A nomogram was developed by integrating predictive and clinical features utilizing the R package RMS (version 7.0-0; http://hbiostat.org/R/rms/). ROC curves were used to assess the nomogram.

GSEA (31) was used to assess whether a predefined geneset exhibited important, consistent variations between two groups. The limma (version 3.50.0) package (32) was applied for gene set variation analysis (GSVA). Gene expression difference between risk cohorts was computed. The clusterProfiler (version 4.2.2) package was applied for GSEA on ranked genes according to log2FC values. For each analysis, gene set permutations were conducted 1,000 times. In the Molecular Signatures Database (MSigDB) Collections, which is A Broad Institute resource for genomics and bioinformatics that organizes annotated gene sets, such as gene sharing functions, locations or regulatory patterns, to support gene expression analysis, pathway studies and omics data interpretation (31,33,34), the gene set was found using c2.cp.kegg.v7.5.1.symbols. Genes with an adjusted P<0.05 were reckoned as prominently enriched.

The GSVA (35) package was used to compare the biological function using c2.cp.kegg.v7.5.1.symobols. The pheatmap package (version 1.0.12; http://CRAN.R-project.org/package=pheatmap) was used to show the results.

The single-sample GSEA (ssGSEA) (36) is an expansion of GSEA that reckons individual enrichment scores for every combination of a sample and gene set. The ssGSEA result revealed how closely the overall upregulated or downregulated genes were in a definite gene set. A total of 28 immune cells were downloaded from the TISIDB (http://cis.hku.hk/TISIDB/index.php) (37) and their relative enrichment was calculated from all gene expression in tumor samples. Immune cell infiltration between risk cohorts was assessed through the ggplot2 (version 3.3.6) (38) package.

The sensitivity of TC patients in two risk subgroups to potential therapeutic drugs was predicted using the oncoPredict (version 0.2) (39) package, based on the IC₅₀ from GDSC (https://www.cnacerrxgene.org/) (40) and clinical gene expression.

The mutation data illustrated the genomic variations landscape. Different clusters were analyzed using the maftools software (version 2.12.0, http://github.com/PoisonAlien/maftools) for somatic variations, such as tumor mutation burden (TMB) (41). The primary driver genes for tumors were typically those frequently mutated genes with the highest mutation frequency in the top 20 (42).

Normal thyroid cell lines (Nthy-ori3-1) and TC cell lines (TPC-1 and BCPAP) were obtained from the medical research building, Fenglin campus, Fudan University (Jiangsu, China). All cells were verified through Short Tandem Repeat (STR) profiling and all STR identification materials have been placed in the supplementary materials. All cells were cultured in RPMI 1640 medium (cat. no. PM150110; Wuhan Pricella Biotechnology Co., Ltd.) or DMEM (cat. no. ZQ-100; Shanghai Zhongqiao Xinzhou Biotechnology Co., Ltd.) with 10% fetal bovine serum (cat. no. F103-01; Vazyme Biotech Co., Ltd.) and 1% penicillin/streptomycin solution (cat. no. CSP006; Shanghai Zhongqiao Xinzhou Biotechnology Co., Ltd.) at 37°C with 5% CO₂. RNA was extracted from the cell line using the FastPure Cell/Tissue Total RNA Isolation Kit V2 (cat. no. R223-01; Vazyme Biotechnology). HiScript II qPT SuperMix II reagents were added according to the manufacturer's manual to reverse transcribe RNA into cDNA (cat. no. Q331-01; Vazyme Biotechnology), then the ChamQ SYBR qPCR Master Mix reagent (cat. no. Q331-02; Vazyme Biotechnology) was added. RT-qPCR was performed using a quantitative fluorescence PCR machine (cat. no. ABI 7500; Applied Biosystems; Thermo Fisher Scientific, Inc.). The thermocycling protocol used was as follows: 95°C for 30 sec; followed by 40 cycles of 95°C for 10 sec, 60°C for 30 sec, 95°C for 15 sec, 60°C for 60 sec and 95°C for 15 sec. After obtaining the CT values, comparisons between multiple groups were made using the 2^−ΔΔCq method (43). All experiments were repeated three times. The primer sequences are shown in Table SII.

SIPA1L2 small interfering RNA (siRNA) was designed and synthesized by Generay Biotech Co., Ltd (Table SII). The logarithmic-phase TPC-1 and BCPAP cells were digested with trypsin for 2 min at 37°C, and the resulting cell suspension was resuspended in antibiotic-free medium and seeded into a 6-well plate at a density of 8×10⁵ cells per well. When the cell density reached ~70%, the Hieff Trans^® liposomal transfection reagent (cat. no. 40802ES03; Shanghai Yesen Biotechnology Co., Ltd.) and RNAi complexes were added for transfection. The procedure involved mixing 5 µl of transfection reagent with 250 µl of serum-free medium or combining 100 pmol RNAi with 250 µl of serum-free medium, then combining these solutions and incubating at room temperature for 20 min. The mixture was then added to the cells, and the cells were incubated at 37°C with 5% CO₂ for 24 h to complete the transfection.

The transfected cells were resuspended and seeded in a 12-well plate with coverslips (3×10⁵ cells/well). After incubating at 37°C with 5% CO₂ for 24 h, the EdU probe (cat. no. C0071S; Beyotime Institute of Biotechnology) was added and cells were incubated at 37°C for 2 h. The cells were then fixed at room temperature with 4% paraformaldehyde for 15 min, permeabilized with 0.5% Triton-X-100 at room temperature for 15 min, stained with Azide 488 at room temperature in the dark for 1 h and subsequently stained with Hoechst 33342 at room temperature in the dark for 30 min. Cell proliferation was observed using a fluorescence microscope (Olympus BX43; Olympus Corportation) using a ×10 objective. The coverslip was divided into four areas: Upper left, lower left, upper right and lower right, and images were captured from these fields.

Transfected cells were resuspended and cultured in a 12-well plate. After 24 h, when cells reached 95% confluence, a cross-shaped scratch was created using a sterile 200 µl pipette tip. After PBS washing, the cells were cultured in a serum-free medium and photographed at 0 and 24 h.

Transfected cells were resuspended in a serum-free medium at 2×10⁵ cells per well. The upper chamber was pre-coated with Matrigel solution (Beyotime Institute of Biotechnology) at 37°C for 3 h. Subsequently, 200 µl of the cell suspension was added in the upper chamber and 600 µl of complete medium was added in the lower chamber. The cells were incubated at 37°C for 24 h. After 24 h, the cells were fixed with 4% paraformaldehyde at room temperature for 30 min and stained with 0.1% crystal violet (Beyotime Institute of Biotechnology) at room temperature for 30 min. Samples were removed from the chamber, placed on a slide and observed using an inverted DMI3000B fluorescence microscope (Leica Microsystems GmbH) with a ×10 objective to assess cell invasion. The chamber was divided into four areas: Upper left, lower left, upper right, and lower right, and images were taken from these fields.

R (version 4.1.2; http://www.R-project.org/) was applied for data analysis and statistical calculations. K-M curves and log-rank tests were adopted for survival analysis. Prognostic variables were assessed using Cox regression analysis. Factors with a greater impact on results were identified using the Boruta algorithm. R program ggplot2 was used to visualize the data and the R survival package was used to calculate OS and risk scores. Significant quantitative differences in variables with normal distribution were analyzed utilizing two-tailed t-tests or one-way ANOVA, with Tukey's correction for multiple comparisons. Wilcoxon and Kruskal-Wallis tests were adopted to assess notable variances in data without normal distribution. P<0.05 was considered to indicate a statistically significant difference.

Cell populations of TC were examined to investigate highly expressed genes using the GSE184362 dataset. After an initial evaluation of quality control and removal of duplicates, single-cell transcriptomes of 56,580 cells were obtained. All cells were clustered into 13 distinct groups (Fig. 2A). Distinct cell types were discovered based on the genetic features of each group using cell type-specific biomarkers listed in Table SIII. A total of seven cell types, containing epithelial cells, macrophages and T cells are displayed in Fig. 2B. Marker gene levels of all groups are shown in the heatmap (Fig. 2C). Distinct genes associated with individual cell types are visualized in dot plots (Fig. 2D). The distribution of cell types in all samples is shown in Fig. 2E.

To explore the primary sites of glutamine metabolism in TC, the present study identified active cell subpopulations at the single-cell level using the expression patterns of glutamine metabolism-related genes. Active cells were determined by the optimal threshold, with cell populations having an AUC value >0.12 considered as high glutamine metabolism activity groups and those with an AUC value <0.12 as low glutamine metabolism activity groups. The results showed that there were 1,751 glutamine metabolism -active cells (Fig. 3A). The AUCell scores of glutamine metabolism for each cell type was further calculated and their differences before and after the occurrence of TC compared. Fig. 3B displays the t-SNE plot of active cells, where epithelial cells showed the highest glutamine metabolism scores among all cells, thus the present study regarded them as the main subjects for subsequent analysis. A transcriptional trajectory was constructed using definitive epithelial cells to identify crucial gene expression patterns that govern TC progression. Differentiated paths were revealed in the trajectory's transcriptional states. Epithelial cells showed five states, in which epithelial cells from the GSM5585102, GSM5585104 and GSM5585121 datasets were enriched in state 5 and the remaining epithelial cells were enriched in states 1–4 (Fig. 3C-F). The genes responsible for the first branching point of TC were investigated to understand the molecular basis. Highly expressed genes in pre-branch cells were mainly enriched in pathways related to inhibiting endopeptidase activity, peptidase activity and receptor-mediated endocytosis. In addition, genes enriched in cytoplasmic translation, adjustment of leukocyte proliferation and myeloid leukocyte migration showed high expression in cell fate 1 (Fig. 3G).

The interactions between these types of cells are shown in Fig. 4A and B. Subsequently, the potential signals and specific pairs of molecules in these seven types of cells were analyzed. Macrophages were the primary source of signals, while epithelial cells were the main recipients (Fig. 4C and D). Epithelial cells potentially communicated through the MK, TWEAK, MIF, ANNEXIN and VEGF signaling pathways. Consequently, the specific signal pairs in epithelial cells were examined. The findings indicated that epithelial cells primarily communicated with B cells through the MIF-(CD74+CXCR4) signaling pathway, while fibroblasts communicated with themselves through the MDK-SDC4 and MDK-NCL signaling pathways (Fig. S1A-C). The ligand MIF was primarily expressed in epithelial cells (Fig. S1D). CD74 and CXCR4 receptors were expressed in B cells. The ligand MDK was primarily expressed in epithelial cells and fibroblasts (Fig. S1E). SDC2 and SDC4 receptors were mostly expressed in epithelial cells. NCL receptor was expressed in various cells.

A total of 2,098 DEGs were uncovered between epithelial cells and other cell types, with P<0.05 and |Log2 fold change|>0.25 (Table SIV and Fig. 5A). A total of 5,825 DEGs were found between TC issues and control thyroid issues (P<0.05, |Log2 fold change|>0.5). The top 10 upregulated and downregulated genes in TC samples are displayed in a heatmap (Table SV, Fig. 5B). The intersection between the two types of DEGs yielded 633 hub genes (Table SVI; Fig. 5C). To explore the biological roles of the central genes, the GO (Table SVII) and KEGG analyses (Table SVIII) were conducted. These genes were enriched in pathways linked to T cell (GO: 0042110), lymphocyte activation (GO: 0051251) and cell adhesion (GO: 0045785; Fig. 5D). The enhanced KEGG pathways comprised the synthesis of thyroid hormone (hsa04918), infection by human T-cell leukemia (hsa05166) and Epstein-Barr virus (hsa05169; Fig. 5E).

DEGs in epithelial cells were compared with those in TC and control groups. Signature genes in epithelial cells were distinguished utilizing the K-M curve with P<0.05. In Table SIX, 63 TC-associated genes were ultimately pinpointed. Boruta was further screened among 63 targets and 17 genes were obtained. The model was constructed by multivariate COX and the genes with a linear relationship with each other were removed by collinearity analysis and six genes were obtained after optimizing the model, with the seed value set to 78. A total of six genes (S100A13, APOE, TUSC3, SNX22, SIPA1L2 and PSMB8) associated with TC prognosis were discovered. 70% of TC samples (n=502) were chosen for the training set (n=348), while 30% of samples were allocated to the validation set (n=154). The model constructed based on these six genes was evaluated by dividing the samples into high- and low-risk categories according to the median score. Risk scores, survival time and levels of six genes in various groups from the training set are displayed in Fig. 6A. The high-risk population exhibited a greater density of death states and increased gene levels. Survival curves were generated for patients in various groups in the training (Fig. 6B) and validation sets (Fig. 6C). High-risk populations had notably worse prognoses than low-risk populations across all cohorts. The effectiveness of the model in forecasting patient outcomes was assessed using ROC curves. The training set showed 0.7, 0.781 and 0.823 AUC values for 1-, 3- and 5-year survival rates, respectively (Fig. 6D). The validation cohort showed AUC values of 0.835, 0.722 and 0.707, respectively (Fig. 6E). To explore the mechanisms underlying these DEGs, GSEA using MSigDB Collection identified that LYSINE_DEGRADATION, PROPANOATE_METABOLISM and ASCORBATE_AND_ALDARATE_METABOLISM were enriched in the high-risk cohorts, while TYPE_I_DIABETES_MELLITUS, ALLOGRAFT_REJECTION and SYSTEMIC_LUPUS_ERYTHEMATOSUS were enriched in the low-risk cohorts (Fig. S2A-F; Table SX). Additionally, GSVA highlighted five pathways with the most notable differences between the high-risk and low-risk populations (Fig. S2G; Table SXI).

The infiltration of 28 immune cells was evaluated by the ssGSEA technique (Table SXII). The proportions of 28 cells were manifested in a stacked bar plot (Fig. 7A). In addition, cell division differed among TC patients. Most immune cells showed positive correlations with others (Fig. 7B). Monocytes were negatively linked with myeloid-derived suppressor cells (Fig. 7B). Noticeable differences were observed between the high-risk and low-risk populations in Activated B cell, Activated CD4 T cell and Activated CD8 T cell (P<0.05; Fig. 7C).

To verify whether the risk score can serve as an independent prognostic factor, univariate Cox regression analysis was performed on the patients' clinical characteristics, such as age, stage and sex. All clinical data were sourced from the TCGA database. The results of the univariate Cox regression analysis indicated that stage, age and risk score were all associated with the patients' prognostic risk (Fig. 8A). Based on the prognostic risk factors obtained from the R package RMS and one-way COX, the present study constructed a nomogram of prognostic diagnosis (Fig. 8D) and detected its diagnostic value by ROC curve and calibration curve. AUC values for survival rates at 1, 3 and 5 years were 0.979, 0.941 and 0.939, respectively (Fig. 8B). The results indicated that the prognostic model had good diagnostic value for the prognosis of patients. The calibration curve also shows minimal deviation from the ideal model, suggesting good predictive performance of the model (Fig. 8C).

The specific genetic mutation in TC was evaluated and the top 20 driver genes with frequent mutations were screened. BRAF was the most frequent mutation, followed by NRAS (Fig. 9A and B). TMB is a vital element in determining response to immunotherapy. Consequently, genetic mutations related to TC were examined. The high-risk population displayed markedly greater TMB than the low-risk population (P<0.05; Fig. 9C). The present study focused on determining the effectiveness of risk scores in predicting chemotherapy response in TC patients. Camptothecin_1003, Dactinomycin_1811, Dactinomycin_1911, Dactolisib_1057, Epirubicin_1511, MG-132_1862, Podophyllotoxin bromide_1825 and Sabutoclax_1849 were investigated for the clinical efficiency of TC treatment (Table SXIII). Low-risk patients may exhibit high sensitivity to Camptothecin_1003 (Fig. 9D), Dactinomycin_1811 (Fig. 9E), Dactinomycin_1911 (Fig. 9F), Dactolisib_1057 (Fig. 9G), Epirubicin_1511 (Fig. 9H), MG-132_1862 (Fig. 9I), Podophyllotoxin bromide_1825 (Fig. 9J) and Sabutoclax_1849 (Fig. 9K). These results suggest that chemotherapy could be a favorable option for low-risk individuals.

The mRNA levels of S100A13, APOE, TUSC3, SNX22, SIPA1L2 and PSMB8 differ between normal thyroid cells (Nthy-ori3-1) and TC cell lines (TPC-1, BCPAP). S100A13, APOE, TUSC3, SNX22 and SIPA1L2 have higher mRNA expression levels in the BCPAP cell line compared to normal thyroid cells, while APOE, TUSC3, SNX22 and SIPA1L2 have higher expression levels in the TPC-1 cell line compared to normal thyroid cells. However, the expression level changes of S100A13 and PSMB8 in the TPC-1 cell line were not statistically significant (Fig. 10A).

Existing literature reports that the TUSC3 and APOE genes play a role in promoting disease progression in TC (44,45). The present study found that the changes in mRNA expression levels of PSMB8 and S100A13 genes in TPC-1 cells were not statistically significant. SNX22 itself is expressed at low levels in TC and is not easily knocked down. Therefore, the present study chose to validate the role of the SIPA1L2 gene in TC cells. After knocking down the SIPA1L2 gene, mRNA level detection revealed that it could be knocked down by more than 50% (Fig. 10B), followed by functional validation. After SIPA1L2 knockdown, EDU proliferation experiments detected enhanced proliferation ability in BCPAP and TPC-1 cell lines (Fig. 10C and D), scratch experiments detected enhanced migration ability in BCPAP and TPC-1 cell lines (Fig. 10E and F) and Transwell chamber experiments found enhanced invasion ability in BCPAP and TPC-1 cell lines (Fig. 10G and H). This indicated that SIPA1L2, as a protective gene, can inhibit the proliferation, migration and invasion abilities of TC cells.