All subject information was collected from the public database and informed consent was gained by default. This study was approved (WDRY2022-K037) by the Ethics Committee of the Renmin Hospital of Wuhan University, China. Ten PTC tissues and paired paracancer normal tissues were used for in vitro validation. Each patient signed an informed consent form.

The RNA expression data in counts format, accompanying clinical information, and mutation data of PTC patients (511 cancer and 59 normal) were collected from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) (17). Gene expression data of normal thyroid tissues (279 normal tissues) were obtained from the Genotype-Tissue Expression (GTEx) database (https://gtexportal.org/home/) (18). The clinical information of 503 PTC patients is summarized in Supplementary Table 1 . A total of 260 genes that associated with lactate metabolism and enrichment analysis data set was obtained from MsigDB (http://www.gsea-msigdb.org/gsea/downloads.jsp) (19). Moreover, 300 chemokines genes and 149 immune checkpoint genes were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov/) (20).

We used the R package limma (21) and the normalizeBetweenArrays function to reduce the batch effects that may exist between or within the two cohorts and merge the TCGA and GTEx RNA expression data (22). Using the limma package, we evaluated the differential expression of 205 lactate metabolism-related genes (LRGs) in PTC vs. normal patients (|Log₂FC|>1, FDR<0.05). The 205 LRGs that were differentially expressed in the TCGA-THCA (23) cohort and the 16,877 lncRNAs were then analyzed for relationships using Pearson correlation analysis (P<0.05, correlation coefficient>0.4), from which 4,032 LRLs were identified. Subsequently, the limma package was also used to evaluate the differential expression of LRLs in PTC vs. normal patients (|Log₂FC|>1, FDR<0.05). Following this, 1,118 differentially expressed LRLs were screened for additional bioinformatic analyses.

The ConsensusClusterPlus (24) package in R was used to cluster the data based on the expression of 1,118 LRLs in TCGA-THCA cohort, using a consensus matrix and a cumulative distribution function (CDF) to calculate the optimal number of clusters. The survival analyses on the three clusters were performed using the “survival” package.

Prognostic LRLs were screened in the TCGA-THCA cohort using univariate Cox regression analysis (P<0.05). To avoid overfitting, a least absolute shrinkage and selection operator (Lasso) regression analysis using the glmnet package (25) was performed. The risk score was calculated using the following formula:

where Coef(i) and x(i) represent the respective LRLs coefficient and expression level.

The aforementioned algorithm was used to calculate the risk score for each PTC patient. The characteristics of the risk model were assessed through AUC evaluation and Kaplan-Meier (KM) survival analysis. A risk score was calculated for each patient and the median risk score for the TCGA-THCA cohort was used to define the high- and low-risk groups. The R package scatterplot3d (26) was used to perform principal component analysis (PCA) and construct the relevant plots. The correlation analyses between clinical traits were performed using the survival (27), survminer (28), timeROC (29), and ggplot2 (30) packages in R, combined with the Cox algorithm, KM survival analysis, univariate Cox analysis, multivariate Cox analysis, and ROC curve. Finally, a nomogram was constructed based on the coefficients of the multivariate Cox regression using the rms (31) package.

Using seven different algorithms (XCELL (32), TIMER (33), QUANTISEQ (34), MCPCOUNTER (35), EPIC, CIBERSORT-ABS, and CIBERSORT (36)), the relationship between risk scores and immune cells was estimated. Using single sample gene set enrichment analysis (ssGSEA) and the estimate algorithm, we further compared the immune function and immunological pathways between the high- and low-risk groups. Using the CIBERSORT algorithm, we evaluated the cellular composition of each sample from the high- and low-risk groups to determine the differences in immune cell infiltration between these groups. In addition, the immune escape and tumor microenvironment (TME) between the two groups were compared using the ggpubr (37) and limma packages. Finally, the relationship between immune checkpoints and chemokines expression, as well as the association between risk scores and DNA Stemness Scores (DNAss)/RNA Stemness Scores (RNAss) in the two risk groups were explored using the limma package.

The TCGA database was used to extract the somatic mutation profiles of 485 TCGA-THCA patients. Using the maftools (38) package, the mutations in both high- and low-risk groups were analyzed. In this study, the tumor mutation burden (TMB) score was computed by dividing the number of mutations by the length of the exon (30Mb). Based on the median TMB score, all samples with somatic mutations were split into two groups: high and low TMB groups. Variance analysis was performed to assess the differences in TMB between the high- and low-risk groups. Using a correlation scatter plot, the association between risk score and TMB was represented. KM analysis was performed to determine survival probability differences between the high and low TMB groups.

Gene Ontology (GO) (39), Kyoto Encyclopedia of Genes and Genomes (KEGG) (40), and Gene Set Enrichment Analysis (GSEA) (41) analyses were performed on the risk model using the clusterProfile (42) package. A p-value<0.05 and FDR<0.05 were considered statistically significant in the GO and KEGG analyses. In the GSEA, the GSEA software (version 4.2.3) was obtained from the GSEA website. Subsequently, the samples were divided into two groups based on risk score, and a subset of c2.cp.kegg.v7.5.1.symbols.gmt was obtained from MsigDB to assess the relevant pathways and molecular mechanisms. The minimum gene set was adjusted to 15 and the maximum gene set to 500, with resampling of 1000x. A P-value<0.05, and an FDR<0.25 were considered statistically significant.

Ten PTC tissues and paired paracancer normal tissues were acquired from PTC patients undergoing tumor excision. Fresh tumor and non-tumor tissues were flash-frozen in liquid nitrogen. Total RNA was obtained using an extraction kit (Servicebio, China). Complementary DNA (cDNA) was obtained using the Servicebio^® RT First Strand cDNA Synthesis Kit. The SYBR Green qPCR Master Mix (Servicebio) was used for RT-PCR. As an internal control, all expression data were standardized to GAPDH using the 2^−ΔΔCt approach. All primers used were outsourced from Servicebio Biotechnology Company (Wuhan, China). Primers for these lncRNAs are included in Supplementary Table 2 .

The R software (v.4.2.0) and GraphPad Prism (9.0.0) software were used to analyze all the data. GSEA software (version 4.2.3) was used for the GSEA enrichment analysis. To analyze the differences between the two groups, the paired sample t-test was used. Statistical significance was evaluated based on p-values less than 0.05.

The study design is illustrated in Figure 1 . Using gene expression data from the TCGA and GTEx datasets, 205 differentially expressed LRGs between thyroid normal and PTC tissues were identified. Moreover, 1,118 differentially expressed LRLs obtained from the correlation analysis (P<0.05, correlation coefficient>0.4) and differential analysis (|Log₂FC|>1, FDR<0.05). Cox regression analysis and Lasso regression analysis were used to obtain the risk prognostic model, which was further used for subsequent clinical, TME, tumor mutation, and enrichment analyses. The expression levels of six LRLs in PTC and normal tissues were verified by RT-PCR.

An evaluation of the expression of 260 LRGs in normal and PTC samples led to the identification of 205 LRGs specific to the TCGA-THCA cohort (|Log₂FC|>1 and P<0.001). Subsequently, 18 significantly different LRGs were selected to shown in the heatmap ( Figure 2A ). In contrast to normal samples, 104 LRGs were upregulated, and 101 LRGs were downregulated in the PTC samples ( Figure 2B ). To further explore the correlation between LRGs, we performed Pearson correlation analyses. The analyses showed that MDH2 exhibited the largest positive association with GOT2 (r=0.85), while NFS1 had the strongest negative correlation with MRPS14 (r=-0.55) ( Figure 2C ). Using the Pearson correlation analyses, 4,032 LRLs linked with LRGs were identified in the TCGA cohort (correlation coefficient>0.4, P<0.05) ( Figure 2D ).

We obtained 1,118 differentially expressed LRLs (|Log₂FC|>1, FDR<0.05) in the TCGA-THCA cohort by evaluating the expression levels of 4,032 LRLs in normal and PTC samples using the limma package ( Figure 3A ). The heatmap shows the expression levels of 20 of these LRLs ( Figure 3B ). To further understand the overall role of LRLs in PTC patients, a consensus clustering analysis was performed on the TCGA-THCA cohort. Duplicate samples were removed, and similar samples were grouped into the same category. To determine the optimal number of LRL clusters, we variously tested the total number of LRL clusters from 2 to 9 and examined the CDF curves of the consensus matrix separately. The values with the largest CDF, but with little correlation change were selected ( Figure 3C, D ). Finally, K=3 was used as the optimal number of clusters based on the consensus matrix ( Figure 3E ). KM curves showed a significantly lower survival probability in cluster 2 compared to that of cluster 1 (P<0.05) ( Figure 3F ).

To further determine the LRLs associated with prognosis among the identified 1,118 differentially expressed LRLs, a univariate Cox analysis was performed. A total of 303 LRLs were screened for further investigation using the limma package ( Figure 4A ). To reduce model overfitting and increase the model’s discriminative ability, the Lasso regression analysis was further applied ( Figure 4B, C ). Subsequently, six LRLs were extracted from the Lasso regression analysis and were used to develop a risk model. The expression and regression coefficients of the six LRLs were then combined to determine the risk scores of the TCGA-THCA patients using the following formula: risk score = (0.455*AC133785.1) + (0.0244*AL138781.1) + (0.1385*AC084871.3) + (0.2779*AL008733.1) + (1.3443*AC245014.3) + (0.1387*AC124276.2) ( Figure 4D ). Based on the median risk score (0.39717-0.39771) of the TCGA-THCA cohort, the patients were separated into the high- and low-risk groups. PCA plots were used to clarify the distribution of LRLs, showing that the median risk score could significantly distinguish the high- and low-risk patients in the TCGA-THCA cohort ( Figure 4E ). The six LRLs that were used to construct the risk model and their related LRGs were then subjected to a Pearson correlation analysis (correlation coefficient>0.4, P<0.05) ( Figure 4F ).

KM survival analysis showed that the overall survival (OS) of the high-risk group was significantly lower than that of the low-risk group (P<0.01) ( Figure 5A ). Furthermore, the risk scores, survival status, and risk genes expression level were also compared between the high- and low-risk groups in the TCGA-THCA cohort ( Figure 5B ). In addition, the AUC for the risk model were 0.979, 0.826, and 0.852 at 1, 3, and 5 years, respectively ( Figure 5C ). Finally, an examination of the ROC curves showed that risk score had the highest AUC value when compared to other clinical traits, showing the best predictive potential ( Figures 5D–F ).

To evaluate whether the risk score and clinical traits were independent prognostic factors in PTC patients, univariate and multivariate Cox analyses were performed. The results of the univariate Cox regression analysis showed that the risk score was substantially linked to prognosis in the TCGA cohort (TCGA cohort: HR=6.265, 95% CI=3.796-10.340). Meanwhile, the multivariate Cox regression analysis showed that the risk score remained an independent predictor of prognosis after adjustment for other covariates (TCGA cohort: HR=174.454, 95% CI=3.799-8011.106) ( Figure 6A, B ). In addition, after deleting the missing data in various clinical traits, we integrated data on seven clinical traits, including age, gender, stage, risk score, Tumor, Node and Metastasis (TNM) stage to construct a nomogram to assess the prognostic significance of these traits in 277 PTC patients. The overall C-index of the model was 0.9713, 95% CI (0.9224-1), p-value=1.6427e-79 ( Figure 6C ). Moreover, the prediction curve for PTC patients was close to the standard curve, indicating that the expected survival at 1-, 3-, and 5-years was similar to the actual survival at these time points ( Figure 6D ).

To further explore the relationship between each clinical trait, survival probability, and risk score, clinical data such as age, gender, stage, and TNM stage were extracted separately. Samples with missing information were removed, and the remaining samples were subjected to Cox regression analysis to determine the correlation of clinical traits with survival probability and risk score, respectively. Our results showed that survival probability was significantly lower in patients >65 years of age compared to those<=65 years (P<0.001) ( Figure 7A ). Moreover, stage I-II, T1-T2, and M0 patients have a significantly higher survival probability than stage III-IV, T3-T4, and M1 patients (P<0.05), respectively. In contrast, gender and N stage did not seem to have a significant effect on survival ( Figures 7B–F ). Furthermore, significant differences in risk scores between the high- and low-risk groups in terms of gender, age, and N stage were observed (P<0.05). More specifically, patients with age >65 years, men, and patients with N1 stage had higher risk scores than patients with age<=65 years, women, and patients with N0 stage, respectively (P<0.05). Finally, stage I and II, T1-T2 stage, and T3-T4 stage did not seem to have a significant effect on risk scores ( Figures 7G–L ).

We evaluated the relationship between immunity and risk score. First, using seven different algorithms (XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, CIBERSORT), we compared the correlation between risk score and various immune cells expression level. Our analyses revealed that endothelial cells, cancer-associated fibroblasts, monocytes, and the remaining CD4+ T cells were positively correlated with risk score ( Figure 8A ). The generated box plots show that the high-risk group had an immune cell proportion that was significantly lower than that of the low-risk group when 16 immune cell types and their associated functions were analyzed using the ssGSEA algorithms ( Figure 8B, C ). Next, the various cellular composition of each sample in the high- and low-risk groups were measured using the CIBERSORE algorithm ( Figure 8D ). Our results showed that the immune score was significantly lower in the high-risk group than in the low-risk group (P<0.05) ( Figure 8E ). Finally, the Tumor Immune Dysfunction and Exclusion (TIDE) results showed that there was no significant difference in immune escape between the high- and low-risk groups ( Figure 8F ). Taken together, the results showed that the low-risk group had higher immune cell infiltration.

Considering the relevance of human chemokines and checkpoint immunotherapy, it is also crucial to note that there were considerable differences in the expression of immunological checkpoints and chemokines between the different risk groups ( Figures 9A–D ). Increasing evidence reveals a link between the increased expression of stemness-related markers in tumor cells and drug resistance, cancer recurrence, and tumor development (43). Our results demonstrated an inverse relationship between risk scores and DNAss/RNAss (P<0.01) ( Figures 9E, F ).

Since mutations play critical roles in cancer development, we focused on the distribution of somatic mutations in the two risk groups. The 15 most commonly mutated genes in both populations were identified ( Figure 10A, B ). No significant difference in TMB was observed between the high- and low-risk groups, and no strong correlation between TMB and risk score was seen ( Figure 10C, D ). KM analysis showed a better survival probability in the low TMB group than in the high TMB group (P<0.001) ( Figure 10E ). Furthermore, patients with low risk and low TMB scores had the best prognosis (P<0.001) ( Figure 10F ), indicating that survival probability is negatively correlated with risk scores and TMB. Overall, immunotherapy for PTC patients with high TMB appears to be ineffective. This may be related to the limited data on tumor mutations in PTC patients and warrants further studies to explore the relationship between TMB and immunotherapy.

To elucidate the underlying biological functions and pathways associated with the six LRLs, a GSEA enrichment analysis was performed. The results showed that the high-risk group was mainly enriched in amino acid metabolism and insulin metabolism pathway, while the low-risk group was not enriched in these pathways ( Figure 11A ). Based on the different gene expression of high- and low-risk group data (|Log₂FC|>1 and P<0.05), 418 risk differentially expressed genes (DEGs) were identified ( Figure 11B ). KEGG enrichment analysis showed that risk DEGs were mainly enriched in amino acid metabolism, glycerol phospholipid metabolism, and thyroid hormone synthesis ( Figure 11C ). Moreover, GO enrichment analysis showed that these risk DEGs were mainly related to extracellular matrix and ion transport ( Figure 11D ). Taken together, the six LRLs were found to be mainly related to amino acid metabolism and extracellular matrix.

To further illustrate the viability of our prognostic model, we analyzed the expression levels of the six LRLs in TCGA-THCA cohort. Consistent with the RNA-Seq data, the results showed that AC084871.3 and AL008733.1 were upregulated, while AC133785.1, AL138781.1, AC245014.3, and AC124276.2 were significantly downregulated in PTC samples compared with normal counterparts. These findings collectively strengthen the validity and dependability of the risk model ( Figure 12A, B ).