The role of lactate metabolism-related LncRNAs in the prognosis, mutation, and tumor microenvironment of papillary thyroid cancer

Article type: Research Article

Keywords: papillary thyroid cancer, long non-coding RNAs, lactate metabolism, risk model, immune microenvironment

Authors: Minqi Xia, Shuo Wang, Li Wang, Yingna Mei, Yi Tu, Ling Gao

Affiliations: Department of Endocrinology & Metabolism, Renmin Hospital of Wuhan University, Wuhan, China; Department of Infection Prevention and Control Office, Renmin Hospital of Wuhan University, Wuhan, China; Department of Breast & Thyroid Surgery, Renmin Hospital of Wuhan University, Wuhan, China

License: Copyright © 2023 Xia, Wang, Wang, Mei, Tu and Gao CC BY 4.0 This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Article links: DOI: 10.3389/fendo.2023.1062317  |  PubMed: 37025405  |  PMC: PMC10070953

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (249 KB)

Introduction

Thyroid carcer is the most prevalent endocrine tumor in US (ref. 1). It has an increasing incidence at an average rate of 4.5% per year (ref. 2), which is faster than any other type of cancer. There are four subtypes of thyroid cancer: papillary thyroid cancer (PTC), follicular thyroid cancer (FTC), medullary thyroid cancer (MTC), and anaplastic thyroid cancer (ATC) (ref. 3). Among these, the PTC accounts for the majority of the growing incidence, which has more than tripled over the previous four decades (ref. 4, ref. 5). Risk stratification is an important component of the PTC treatment strategy, but so far only age at diagnosis is considered the most important risk factor (ref. 6). The widespread use of thyroid ultrasound in physical examination has led to the increased detection rate of thyroid cancer annually (ref. 7). Fine needle aspiration (FNA) is currently the most effective method of identifying PTC, and it has been used extensively in clinical practice (ref. 8). Unfortunately, although FNA generally proves to be effective, there are still many missed diagnoses. In addition, while there have been some advances in diagnostics and treatment, prompt diagnosis and predict prognosis of patient remain inadequate. Therefore, molecular methods are crucial for the early diagnosis of PTC.

Metabolic reprogramming can be caused by mutations in critical metabolic enzymes and is one of the hallmarks of cancer. Targeting cellular metabolism is a potential cancer therapy strategy (ref. 9, ref. 10). The Warburg effect refers to the tendency of cancer cells to change their gluconeogenic mode from oxidative phosphorylation to glycolysis and to produce lactate under hypoxic conditions (ref. 11). The Warburg effect shields tumor cells from hypoxic conditions and provides a plethora of precursors for the production of nucleic acids (ref. 12), fatty acids, and phospholipids that sustain tumor cell proliferation (ref. 13).

LncRNAs are non-coding RNAs that exceed 200 nucleotides in length. Recent studies have demonstrated that lncRNAs play a crucial role in carcinogenesis and metastasis, and abnormal expression of lncRNAs have been identified in a range of malignancies. These dysregulated lncRNAs are regarded to be potential biomarkers or therapeutic targets (ref. 14–ref. 16). However, the role of lactate metabolism -related lncRNAs (LRLs) in PTC development remains to be discovered. In this study, for the first time, a risk model containing six LRLs was constructed and verified with an area under the curve (AUC) value of 0.852-0.979. Moreover, the lncRNA-based risk model was tested with other malignant characteristics such as immune infiltration and tumor somatic mutation. Overall, our work shows that the six-LRLs risk model may serve as a useful prognostic tool for the assessment of PTC patients, and the identified six LRLs can be further studied for possible targeted therapy applications.

Methods

Ethical declaration

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.

Datasets and preprocessing of data

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/) (ref. 17). Gene expression data of normal thyroid tissues (279 normal tissues) were obtained from the Genotype-Tissue Expression (GTEx) database (https://gtexportal.org/home/) (ref. 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) (ref. 19). Moreover, 300 chemokines genes and 149 immune checkpoint genes were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov/) (ref. 20).

Extraction of LRLs

We used the R package limma (ref. 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 (ref. 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 (ref. 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.

Clustering analysis of LRLs

The ConsensusClusterPlus (ref. 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.

Development of a risk prognosis model

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 (ref. 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 (ref. 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 (ref. 27), survminer (ref. 28), timeROC (ref. 29), and ggplot2 (ref. 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 (ref. 31) package.

Investigating immune infiltration

Using seven different algorithms (XCELL (ref. 32), TIMER (ref. 33), QUANTISEQ (ref. 34), MCPCOUNTER (ref. 35), EPIC, CIBERSORT-ABS, and CIBERSORT (ref. 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 (ref. 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.

Evaluation of somatic mutations

The TCGA database was used to extract the somatic mutation profiles of 485 TCGA-THCA patients. Using the maftools (ref. 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.

GO, KEGG, and GSEA enrichment analyses

Gene Ontology (GO) (ref. 39), Kyoto Encyclopedia of Genes and Genomes (KEGG) (ref. 40), and Gene Set Enrichment Analysis (GSEA) (ref. 41) analyses were performed on the risk model using the clusterProfile (ref. 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.

Real-time quantitative polymerase chain reaction (RT-PCR)

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.

Statistical analyses

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.

Results

Study design

The study design is illustrated in f1. 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.

LRGs expression levels in the TCGA and GTEx cohorts

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 (f2). In contrast to normal samples, 104 LRGs were upregulated, and 101 LRGs were downregulated in the PTC samples (f2). 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) (f2). Using the Pearson correlation analyses, 4,032 LRLs linked with LRGs were identified in the TCGA cohort (correlation coefficient>0.4, P<0.05) (f2).

Cluster analysis of LRLs

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 (f3). The heatmap shows the expression levels of 20 of these LRLs (f3). 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 (f3). Finally, K=3 was used as the optimal number of clusters based on the consensus matrix (f3). KM curves showed a significantly lower survival probability in cluster 2 compared to that of cluster 1 (P<0.05) (f3).

Evaluation of LRLs and construction of a risk prognostic model

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 (f4). To reduce model overfitting and increase the model’s discriminative ability, the Lasso regression analysis was further applied (f4). 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) (f4). 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 (f4). 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) (f4).

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) (f5). 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 (f5). In addition, the AUC for the risk model were 0.979, 0.826, and 0.852 at 1, 3, and 5 years, respectively (f5). 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 (f5).

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) (f6). 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 (f6). 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 (f6).

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) (f7). 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 (f7). 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 (f7).

Tumor-infiltrating immune cell profiles

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 (f8). 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 (f8). Next, the various cellular composition of each sample in the high- and low-risk groups were measured using the CIBERSORE algorithm (f8). Our results showed that the immune score was significantly lower in the high-risk group than in the low-risk group (P<0.05) (f8). 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 (f8). Taken together, the results showed that the low-risk group had higher immune cell infiltration.

Analysis of immune checkpoints and chemokines

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 (f9). Increasing evidence reveals a link between the increased expression of stemness-related markers in tumor cells and drug resistance, cancer recurrence, and tumor development (ref. 43). Our results demonstrated an inverse relationship between risk scores and DNAss/RNAss (P<0.01) (f9).

Evaluation of somatic mutations in tumors

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 (f10). 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 (f10). KM analysis showed a better survival probability in the low TMB group than in the high TMB group (P<0.001) (f10). Furthermore, patients with low risk and low TMB scores had the best prognosis (P<0.001) (f10), 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.

Enrichment analysis

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 (f11). 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 (f11). KEGG enrichment analysis showed that risk DEGs were mainly enriched in amino acid metabolism, glycerol phospholipid metabolism, and thyroid hormone synthesis (f11). Moreover, GO enrichment analysis showed that these risk DEGs were mainly related to extracellular matrix and ion transport (f11). Taken together, the six LRLs were found to be mainly related to amino acid metabolism and extracellular matrix.

Expression analysis of LRLs

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 (f12).

Discussion

To our knowledge, this is the first study to investigate the relationship between LRLs and PTC. In this study, we used multiple public databases to explore the predictive potential of novel biomarkers in PTC. A total of 1,118 LRLs were initially identified by differential expression analysis on PTC and normal samples. The top six LRLs (AC084871.3, AC133785.1, AL138781.1, AL008733.1, AC245014.3, and AC124276.2) were identified using Lasso and Cox regression analyses. The risk model for the six-LRLs signature accurately predicted the OS of PTC patients, demonstrating its potential as an independent prognostic model. Meanwhile, the AUC value of the risk model reached 0.979, which was significantly higher than the prognostic model constructed in previous studies (ref. 44–ref. 46), indicating that our six-LRLs signature is a better prediction model. Based on the risk model, differences in enrichment pathways, immune function, immune escape and mutations were compared between the high- and low-risk groups. Our findings can contribute to the design of small molecule targeting drugs based on the differences in the expression of immune infiltration, immune checkpoints, and chemokines between the two risk groups, with the potential to improve chemotherapy, radiotherapy, and even immunotherapy. Our model also provides a new strategy for risk score-based PTC prognostic assessment.

In this study, we identified three clusters based on the expression of 1,181 LRLs. Cluster 2 had a worser OS compared to cluster 1, and there was not significance difference in the effect on survival probability between the cluster1 and cluster3, cluster2 and cluster3. Unfortunately, the results showed that the cluster analysis did not seem to separate the three clusters, hence further risk modeling was carried out. Our risk model was composed of six LRLs (AC084871.3, AC133785.1, AL138781.1, AL008733.1, AC245014.3, and AC124276.2), and the risk score of each patient was calculated based on the expression level of these LRLs. However, the roles of these LRLs in other diseases or cancer types have not yet been reported. Early studies have shown that an increasing number of LRLs are associated with tumor growth and immune escape. For instance, decreases lactate generation by inhibiting PKM2 oligomerization, hence lowering cancer growth and macrophage polarization (ref. 47).The lncNSPL facilitated influenza virus a immunological escape by inhibiting trim25-mediated RIG-I ubiquitination of the k63 linkage (ref. 48). Another study has found that during oral carcinogenesis, LPS/TLR4-activated lncRNA IFITM4P can promote immune escape by upregulating PD-L1 through a dual mechanism (ref. 49). Therefore, the six LRLs we found are promising biomarkers in the future. Next, our further research showed that risk score was considered as a risk factor for the prognosis of TCGA-THCA patients. By combining our analysis of clinical traits, we also found that age and risk score were independent risk factors for the prognosis of PTC patients (HR>1, P<0.01). Moreover, age and risk scores were significantly associated with survival probability. This is consistent with the findings of previous studies (ref. 50).

In addition to examining the impact of risk model on the survival prognosis of PTC patients, we also investigated the effect on the TME of PTC patients. In our research, ssGSEA analysis revealed that 12 types of immune cells were significantly reduced, while 13 types of immune signals were significantly suppressed in the high-risk group. Consistent with the results of our immunological studies, there are considerable disparities in the immune responses between the high- and low-risk groups. One of the most prevalent innate immune cells that act as an antitumor response is the macrophage (ref. 51). Macrophages play a crucial part in the progression or resolution of PTC, which is often surrounded by a significant number of immunological “reactive” cells (ref. 52–ref. 54). Traditionally, macrophage can become M1-phenotype and can enhance both innate and adaptive immunity when activated by a variety of environmental stimuli like bacterial LPS and IFN-γ (ref. 51). M2-like macrophages suppress the immune system by producing immunosuppressive molecules such as HLA-G, IL-10, and TGF-β (ref. 55). Specifically, the M2/repair macrophages promote tumor development, while the M1/kill macrophages block or delay tumor growth (ref. 56). Our study showed that macrophage expression was negatively correlated with risk score, meaning that as risk score increased, macrophage expression gradually decreased, resulting in a significantly lower macrophage immune score in the high-risk group than in the low-risk group. Therefore, we infer that the significantly lower survival probability of patients in the high-risk group compared to those in the low-risk group may be related to the significantly lower expression of M1/kill macrophages. Studies have shown that pro-tumor M2 macrophages can be changed into anti-tumor M1 macrophages using some immunotherapeutic techniques (ref. 57, ref. 58). The macrophage immune score in the high-risk group was lower than that in the low-risk group, which may be related to the significantly decreased M1 macrophage expression in the macrophages, resulting in a higher mortality in the high-risk group. Our work may give a fresh perspective on PTC function and immune mechanism processes.

Some of the upregulated immune checkpoints, including CXCR4, are expected to become new immunotherapeutic targets, which can achieve anti-tumor effects together with PD-1/PD-L1 inhibitors. CXCR4 is overexpressed in more than 20 cancer types (ref. 59–ref. 62). It has been reported that CXCL12 can mobilize cancer cells via CXCR4 (ref. 63). NF-κB stimulates the expression of the chemokine receptor CXCR4 to encourage the migration and metastasis of breast cancer cells (ref. 64). When CXCR4 receptors are present in large numbers on cancer cells, it can ensure the migration of cancer cells, thus laying the foundation for cancer metastasis (ref. 65). Tailored drug strategies for cancer patients may prove to be important. These results imply that this prognostic characteristic may aid in elucidating the regulatory mechanisms of tumor immunity and provide novel insights for future TME research.

Tumor mutation analysis showed that gene mutations in the high- and low-risk groups were mainly concentrated in missense mutations of the BRAF gene. Furthermore, there was no significant difference in the TMB score between high- and low-risk groups. TMB values can vary considerably between tumor types as well as within the same tumor type. In the TCGA-THCA cohort, there was little difference in TMB between the high- and low-risk groups, suggesting that there may be insignificant effects of single immunosuppressive therapy on the treatment of the 2 groups. TMB may have its greatest utility when combined with other biomarkers such as PD-L1 and T-cell inflammatory markers.

The GSEA showed that the main functions of the LRLs in the high-risk group were mainly enriched in amino acid metabolism and insulin metabolism. In addition, KEGG analysis revealed that the distinctions between the high- and low-risk groups centered on amino acid metabolism, glycerol phospholipid metabolism, and thyroid hormone synthesis. These were also true for the GO enrichment analysis, with an additional concentration on the extracellular matrix. The amino acid metabolism is distinct since it connects both carbon and nitrogen metabolism, and amino acids are important intermediates in many metabolic processes (ref. 66). It is established that most amino acids undergo a series of transformations to produce sugars, which in turn can produce lactic acid under anaerobic conditions. In other words, amino acid metabolism is closely related to lactate metabolism. Meanwhile, the extracellular matrix invasion is crucial in the early stage of the metastatic cascade. Matrix metalloproteinase inhibition can slow the deterioration of the extracellular matrix and prevents the spread of primary tumor cells (ref. 67, ref. 68). It has been reported that extracellular matrix signaling protein ADAM22 can facilitate distant disease spread in vivo and is a crucial mechanism in the metastasis of breast cancer (ref. 69). High expression of ECM-related gene COL6A1 predicts poor prognosis and poor immunotherapy response in bladder cancer (ref. 70). Overall, there are altered expression patterns of extracellular matrix-associated genes in metastatic cancers, and extracellular matrix-associated genes are expected to be markers of aggressive, growing tumors.

PTC continues to be treated with surgery as the first-line therapy, but immunotherapy is gaining importance for refractory thyroid cancer (ref. 71, ref. 72). In the United States, there were about 62,450 cases of differentiated thyroid cancer diagnosed in 2015 (ref. 73). Moreover, it was estimated that about 25-50% of patients with locally advanced or metastatic disease would have a refractory response to radioiodine (ref. 74). Although well-differentiated thyroid cancer responds well to radioiodine therapy and often has favorable treatment outcomes, the management of refractory differentiated thyroid cancer remains a major challenge and is a leading cause of death in thyroid cancer patients, with no effective treatment (ref. 75, ref. 76). To our knowledge, this is the first study to construct a predictive risk model for PTC patients based on six LRLs. We also systematically evaluated the functional enrichment, immune cell infiltration, and tumor somatic mutations of the model, adding rigor to this study. Despite its numerous benefits, the present study has several drawbacks. First, all data used are from retrospective research data, and additional prospective, large-scale, multicenter trials are required to give more compelling evidence for the therapeutic implications of our study. Moreover, since the LRLs screened in this work have not yet been subjected to relevant carcinogenic mechanisms, additional research is required to determine their regulatory processes.

Conclusion

The six-LRLs risk model developed in this work is a consistent and strong predictor of PTC patient survival outcomes. This characteristic was also closely associated with TMB and immune checkpoints can be used as a treatment target for future immunotherapy or drug development.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement

The studies involving human participants were reviewed and approved by Ethics Committee of the Renmin Hospital of Wuhan University. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin. The animal study was reviewed and approved by Ethics Committee of the Renmin Hospital of Wuhan University.

Author contributions

MX and YM performed the data analysis and wrote the manuscript. SW,LW and LG participated in the study design,construction of figures, and manuscript writing and revision. YT provided the clinical samples for gene expression analysis and helped to prepare the manuscript. All authors contributed to the article and approved the submitted version.

References

1. SI Sherman. Thyroid carcinoma.. Lancet (, 2003. [DOI]
2. N Pozdeyev, LM Gay, ES Sokol, R Hartmaier, KE Deaver, S Davis. Genetic analysis of 779 advanced differentiated and anaplastic thyroid cancers.. Clin Cancer Res (, 2018. [DOI]
3. V Subbiah, RJ Kreitman, ZA Wainberg, JY Cho, JHM Schellens, JC Soria. Dabrafenib and trametinib treatment in patients with locally advanced or metastatic BRAF V600-mutant anaplastic thyroid cancer.. J Clin Oncol (, 2018. [DOI | PubMed]
4. 4 SEER Program (National Cancer Institute (U.S.))National Center for Health Statistics (U.S.)National Cancer Institute (U.S.)Surveillance ProgramNational Cancer Institute (U.S.)Cancer Statistics Branch. SEER cancer statistics review. NIH publication (1993). Bethesda, Md.: U.S. Dept. of Health and Human Services, Public Health Service, National Institutes of Health, National Cancer Institute.
5. L Davies, HG Welch. Current thyroid cancer trends in the united states.. JAMA Otolaryngol Head Neck Surg (, 2014. [DOI]
6. X Shen, G Zhu, R Liu, D Viola, R Elisei, E Puxeddu. Patient age-associated mortality risk is differentiated by BRAF V600E status in papillary thyroid cancer.. J Clin Oncol (, 2018. [DOI]
7. WJ Moon, JH Baek, SL Jung, DW Kim, EK Kim, JY Kim. Ultrasonography and the ultrasound-based management of thyroid nodules: Consensus statement and recommendations.. Korean J Radiol (, 2011. [DOI | PubMed]
8. M Bongiovanni, A Spitale, WC Faquin, L Mazzucchelli, ZW Baloch. The Bethesda system for reporting thyroid cytopathology: A meta-analysis.. Acta Cytol (, 2012. [DOI]
9. V Ruiz-Rodado, TM Malta, T Seki, A Lita, T Dowdy, O Celiku. Metabolic reprogramming associated with aggressiveness occurs in the G-CIMP-high molecular subtypes of IDH1mut lower grade gliomas.. Neuro Oncol (, 2020. [DOI]
10. PH Chen, L Cai, K Huffman, C Yang, J Kim, B Faubert. Metabolic diversity in human non-small cell lung cancer cells.. Mol Cell (, 2019. [DOI | PubMed]
11. B Faubert, A Solmonson, RJ DeBerardinis. Metabolic reprogramming and cancer progression.. Science (, 2020. [DOI]
12. MG Vander Heiden, LC Cantley, CB Thompson. Understanding the warburg effect: the metabolic requirements of cell proliferation.. Science (, 2009. [DOI]
13. A Bononi, H Yang, C Giorgi, S Patergnani, L Pellegrini, M Su. Germline BAP1 mutations induce a warburg effect.. Cell Death Differ (, 2017. [DOI]
14. WJ Cao, HL Wu, BS He, YS Zhang, ZY Zhang. Analysis of long non-coding RNA expression profiles in gastric cancer.. World J Gastroenterol (, 2013. [DOI]
15. G Yu, W Yao, J Wang, X Ma, W Xiao, H Li. LncRNAs expression signatures of renal clear cell carcinoma revealed by microarray.. PloS One (, 2012. [DOI | PubMed]
16. X Ji, Z Liu, J Gao, X Bing, D He, W Liu. N(6)-methyladenosine-modified lncRNA LINREP promotes glioblastoma progression by recruiting the PTBP1/HuR complex.. Cell Death Differ (, 2023. [DOI]
17. Comprehensive genomic characterization defines human glioblastoma genes and core pathways.. Nature (, 2008. [DOI]
18. C Huang, H Yi, Y Zhou, Q Zhang, X Yao. Pan-cancer analysis reveals SH3TC2 as an oncogene for colorectal cancer and promotes tumorigenesis via the MAPK pathway.. Cancers (Basel) (, 2022. [DOI]
19. W Ning, A Acharya, S Li, G Schmalz, S Huang. Identification of key pyroptosis-related genes and distinct pyroptosis-related clusters in periodontitis.. Front Immunol (, 2022. [DOI | PubMed]
20. GG Nair, JS Liu, HA Russ, S Tran, MS Saxton, R Chen. Author correction: Recapitulating endocrine cell clustering in culture promotes maturation of human stem-cell-derived beta cells.. Nat Cell Biol (, 2019. [DOI | PubMed]
21. ME Ritchie, B Phipson, D Wu, Y Hu, CW Law, W Shi. Limma powers differential expression analyses for RNA-sequencing and microarray studies.. Nucleic Acids Res (, 2015. [DOI | PubMed]
22. M Zhang, Y Wang, Y Wang, L Jiang, X Li, H Gao. Integrative analysis of DNA methylation and gene expression to determine specific diagnostic biomarkers and prognostic biomarkers of breast cancer.. Front Cell Dev Biol (, 2020. [DOI | PubMed]
23. R Robb, L Yang, C Shen, AR Wolfe, A Webb, X Zhang. Inhibiting BRAF oncogene-mediated radioresistance effectively radiosensitizes BRAF(V600E)-mutant thyroid cancer cells by constraining DNA double-strand break repair.. Clin Cancer Res (, 2019. [DOI]
24. MD Wilkerson, DN Hayes. ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking.. Bioinformatics (, 2010. [DOI]
25. S Abelson, G Collord, SWK Ng, O Weissbrod, N Mendelson Cohen, E Niemeyer. Prediction of acute myeloid leukaemia risk in healthy individuals.. Nature (, 2018. [DOI]
26. J Yang, DJ Ryan, W Wang, JC Tsang, G Lan, H Masaki. Establishment of mouse expanded potential stem cells.. Nature (, 2017. [DOI]
27. C Engblom, C Pfirschke, R Zilionis, J Da Silva Martins, SA Bos, G Courties. Osteoblasts remotely supply lung tumors with cancer-promoting SiglecF(high) neutrophils.. Science (, 2017. [DOI]
28. HI Scher, RP Graf, NA Schreiber, A Jayaram, E Winquist, B McLaughlin. Assessment of the validity of nuclear-localized androgen receptor splice variant 7 in circulating tumor cells as a predictive biomarker for castration-resistant prostate cancer.. JAMA Oncol (, 2018. [DOI]
29. P Blanche, JF Dartigues, H Jacqmin-Gadda. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks.. Stat Med (, 2013. [DOI]
30. DM Popescu, RA Botting, E Stephenson, K Green, S Webb, L Jardine. Decoding human fetal liver haematopoiesis.. Nature (, 2019. [DOI]
31. EM Thompson, T Hielscher, E Bouffet, M Remke, B Luu, S Gururangan. Prognostic value of medulloblastoma extent of resection after accounting for molecular subgroup: A retrospective integrated clinical and molecular analysis.. Lancet Oncol (, 2016. [DOI]
32. D Aran, Z Hu, AJ Butte. xCell: Digitally portraying the tissue cellular heterogeneity landscape.. Genome Biol (, 2017. [DOI | PubMed]
33. M Uhlen, P Oksvold, L Fagerberg, E Lundberg, K Jonasson, M Forsberg. Towards a knowledge-based human protein atlas.. Nat Biotechnol (, 2010. [DOI]
34. L Wang, RP Sebra, JP Sfakianos, K Allette, W Wang, S Yoo. A reference profile-free deconvolution method to infer cancer cell-intrinsic subtypes and tumor-type-specific stromal profiles.. Genome Med (, 2020. [DOI | PubMed]
35. E Becht, NA Giraldo, L Lacroix, B Buttard, N Elarouci, F Petitprez. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression.. Genome Biol (, 2016. [DOI | PubMed]
36. X Tekpli, T Lien, AH Rossevold, D Nebdal, E Borgen, HO Ohnstad. An independent poor-prognosis subtype of breast cancer defined by a distinct tumor immune microenvironment.. Nat Commun (, 2019. [DOI | PubMed]
37. L Ferrafiat, D Pflieger, J Singh, M Thieme, M Bohrer, C Himber. The NRPD1 n-terminus contains a pol IV-specific motif that is critical for genome surveillance in arabidopsis.. Nucleic Acids Res (, 2019. [DOI]
38. A Mayakonda, DC Lin, Y Assenov, C Plass, HP Koeffler. Maftools: Efficient and comprehensive analysis of somatic variants in cancer.. Genome Res (, 2018. [DOI]
39. N Kourtis, C Lazaris, K Hockemeyer, JC Balandran, AR Jimenez, J Mullenders. Oncogenic hijacking of the stress response machinery in T cell acute lymphoblastic leukemia.. Nat Med (, 2018. [DOI]
40. S Rao, L Mondragon, B Pranjic, T Hanada, G Stoll, T Kocher. AIF-regulated oxidative phosphorylation supports lung cancer development.. Cell Res (, 2019. [DOI]
41. S Pei, M Minhajuddin, B Adane, N Khan, BM Stevens, SC Mack. AMPK/FIS1-mediated mitophagy is required for self-renewal of human AML stem cells.. Cell Stem Cell (, 2018. [DOI | PubMed]
42. FF Hoyer, K Naxerova, MJ Schloss, M Hulsmans, AV Nair, P Dutta. Tissue-specific macrophage responses to remote injury impact the outcome of subsequent local immune challenge.. Immunity (, 2019. [DOI | PubMed]
43. W Ouyang, Y Jiang, S Bu, T Tang, L Huang, M Chen. A prognostic risk score based on hypoxia-, immunity-, and epithelialto-mesenchymal transition-related genes for the prognosis and immunotherapy response of lung adenocarcinoma.. Front Cell Dev Biol (, 2021. [DOI | PubMed]
44. ZX Cao, X Weng, JS Huang, X Long. Receptor-ligand pair typing and prognostic risk model for papillary thyroid carcinoma based on single-cell sequencing.. Front Immunol (, 2022. [DOI | PubMed]
45. B Han, X Yang, DK Hosseini, P Luo, M Liu, X Xu. Development and validation of a survival model for thyroid carcinoma based on autophagy-associated genes.. Aging (Albany NY) (, 2020. [DOI]
46. Y Lin, F Gan, X He, H Deng, Y Li. Identification of ferroptosis-associated long noncoding RNA prognostic model and tumor immune microenvironment in thyroid cancer.. J Immunol Res (, 2022. [DOI | PubMed]
47. K Zhao, X Wang, D Zhao, Q Lin, Y Zhang, Y Hu. lncRNA HITT inhibits lactate production by repressing PKM2 oligomerization to reduce tumor growth and macrophage polarization.. Res (Wash D C) (, 2022. [DOI]
48. J Jiang, Y Li, Z Sun, L Gong, X Li, F Shi. LncNSPL facilitates influenza a viral immune escape by restricting TRIM25-mediated K63-linked RIG-I ubiquitination.. iScience (, 2022. [DOI | PubMed]
49. L Shi, Y Yang, M Li, C Li, Z Zhou, G Tang. LncRNA IFITM4P promotes immune escape by up-regulating PD-L1 via dual mechanism in oral carcinogenesis.. Mol Ther (, 2022. [DOI]
50. ZH Yu, ST Feng, D Zhang, XC Cao, Y Yu, X Wang. The functions and prognostic values of m6A RNA methylation regulators in thyroid carcinoma.. Cancer Cell Int (, 2021. [DOI | PubMed]
51. W Chen, J Jiang, W Xia, J Huang. Tumor-related exosomes contribute to tumor-promoting microenvironment: An immunological perspective.. J Immunol Res (, 2017. [DOI | PubMed]
52. A Fiumara, A Belfiore, G Russo, E Salomone, GM Santonocito, O Ippolito. In situ evidence of neoplastic cell phagocytosis by macrophages in papillary thyroid cancer.. J Clin Endocrinol Metab (, 1997. [DOI]
53. M Ryder, RA Ghossein, JC Ricarte-Filho, JA Knauf, JA Fagin. Increased density of tumor-associated macrophages is associated with decreased survival in advanced thyroid cancer.. Endocr Relat Cancer (, 2008. [DOI]
54. S Scarpino, A Stoppacciaro, F Ballerini, M Marchesi, M Prat, MC Stella. Papillary carcinoma of the thyroid: Hepatocyte growth factor (HGF) stimulates tumor cells to release chemokines active in recruiting dendritic cells.. Am J Pathol (, 2000. [DOI]
55. CY Wang, CY Li, HP Hsu, CY Cho, MC Yen, TY Weng. PSMB5 plays a dual role in cancer development and immunosuppression.. Am J Cancer Res (, 2017
56. S Imam, P Dar, R Paparodis, K Almotah, A Al-Khudhair, SA Hasan. Nature of coexisting thyroid autoimmune disease determines success or failure of tumor immunity in thyroid cancer.. J Immunother Cancer (, 2019. [DOI | PubMed]
57. C Guiducci, AP Vicari, S Sangaletti, G Trinchieri, MP Colombo. Redirecting in vivo elicited tumor infiltrating macrophages and dendritic cells towards tumor rejection.. Cancer Res (, 2005. [DOI]
58. P Sinha, VK Clements, S Ostrand-Rosenberg. Reduction of myeloid-derived suppressor cells and induction of M1 macrophages facilitate the rejection of established metastatic disease.. J Immunol (, 2005. [DOI]
59. X Sun, G Cheng, M Hao, J Zheng, X Zhou, J Zhang. CXCL12 / CXCR4 / CXCR7 chemokine axis and cancer progression.. Cancer Metastasis Rev (, 2010. [DOI]
60. D Mukherjee, J Zhao. The role of chemokine receptor CXCR4 in breast cancer metastasis.. Am J Cancer Res (, 2013. [PubMed]
61. P Kukreja, AB Abdel-Mageed, D Mondal, K Liu, KC Agrawal. Up-regulation of CXCR4 expression in PC-3 cells by stromal-derived factor-1alpha (CXCL12) increases endothelial adhesion and transendothelial migration: role of MEK/ERK signaling pathway-dependent NF-kappaB activation.. Cancer Res (, 2005. [DOI]
62. B Furusato, A Mohamed, M Uhlen, JS Rhim. CXCR4 and cancer.. Pathol Int (, 2010. [DOI | PubMed]
63. A Muller, B Homey, H Soto, N Ge, D Catron, ME Buchanan. Involvement of chemokine receptors in breast cancer metastasis.. Nature (, 2001. [DOI]
64. G Helbig, KW Christopherson, P Bhat-Nakshatri, S Kumar, H Kishimoto, KD Miller. NF-kappaB promotes breast cancer cell migration and metastasis by inducing the expression of the chemokine receptor CXCR4.. J Biol Chem (, 2003. [DOI]
65. M Azizi, H Dianat-Moghadam, R Salehi, M Farshbaf, D Iyengar, S Sau. Interactions between tumor biology and targeted nanoplatforms for imaging applications.. Adv Funct Mater (, 2020. [DOI]
66. D Schulze-Siebert, D Heineke, H Scharf, G Schultz. Pyruvate-derived amino acids in spinach chloroplasts: Synthesis and regulation during photosynthetic carbon metabolism.. Plant Physiol (, 1984. [DOI]
67. WG Stetler-Stevenson. Tissue inhibitors of metalloproteinases in cell signaling: Metalloproteinase-independent biological activities.. Sci Signal (, 2008. [DOI | PubMed]
68. K Kessenbrock, V Plaks, Z Werb. Matrix metalloproteinases: regulators of the tumor microenvironment.. Cell (, 2010. [DOI | PubMed]
69. S Charmsaz, B Doherty, S Cocchiglia, D Vareslija, A Marino, N Cosgrove. ADAM22/LGI1 complex as a new actionable target for breast cancer brain metastasis.. BMC Med (, 2020. [DOI | PubMed]
70. X Zhang, J Liu, X Yang, W Jiao, C Shen, X Zhao. High expression of COL6A1 predicts poor prognosis and response to immunotherapy in bladder cancer.. Cell Cycle (, 2023. [DOI]
71. ZT Al-Salama, YY Syed, LJ Scott. Lenvatinib: A review in hepatocellular carcinoma.. Drugs (, 2019. [DOI]
72. BC Bertol, ES Bales, JD Calhoun, A Mayberry, ML Ledezma, SB Sams. Lenvatinib plus anti-PD-1 combination therapy for advanced cancers: Defining mechanisms of resistance in an inducible transgenic model of thyroid cancer.. Thyroid (, 2022. [DOI]
73. RL Siegel, KD Miller, A Jemal. Cancer statistics, 2015.. CA Cancer J Clin (, 2015. [DOI | PubMed]
74. SP Hodak, SE Carty. Radioiodine-resistant differentiated thyroid cancer: hope for the future.. Oncol (Williston Park) (, 2009
75. RL Brown, JA de Souza, EE Cohen. Thyroid cancer: Burden of illness and management of disease.. J Cancer (, 2011. [DOI]
76. EL Mazzaferri, SM Jhiang. Long-term impact of initial surgical and medical therapy on papillary and follicular thyroid cancer.. Am J Med (, 1994. [DOI]