HTseq-FPKM and HTseq-count data of TC samples were obtained from TCGA (https://cancergenome.nih.gov/, accessed January 24, 2022). After removing TC samples of other pathological types, a total of 502 PTC samples and 58 adjacent normal samples as well as the corresponding clinical information were obtained from TCGA.

The CCRG sets were extracted from MsigDB (http://www.gsea-msigdb.org/gsea/index.jsp, accessed January 24, 2022). Differentially expressed genes (DEGs) between PTC and normal tissues were screened using DEseq2 (|log Fold Change (FC)| > 1, p.adjust< 0.05). Next, 266 differentially expressed CCRGs were identified from the intersection of CCRGs and DEGs. CCR-lncRNAs were acquired by Pearson correlation coefficient analysis between the 266 CCRGs and differentially expressed lncRNAs (DElncRNAs) in the PTC samples, and the Benjamini–Hochberg (BH) method was used for p-value correction (|R|> 0.4, BH< 0.05).

Post-operative recurrence and metastasis are the main factors leading to the poor prognosis of PTC; accordingly, the PFI was chosen as the endpoint of this study, instead of overall survival. A total of 498 PTC samples with complete PFI information were included after removing samples with incomplete clinical information and prognosis less than 30 days. The samples were categorized into training (249) and test (249) cohorts at a 1:1 ratio by the “caret” package. Univariate Cox regression analyses were primarily used to screen CCR-lncRNAs that were associated with PFI in the training cohort (p< 0.05). To reduce noise caused by gene interactions and co-expression patterns, we applied LASSO-Cox regression to filter lncRNAs; the formula is as follows:

where the β is the regression coefficient, x is the gene expression level, and π is an adjusted parameter decided by 10-fold cross validation. Parameter optimization was repeated 1000 times for all lncRNAs. These lncRNAs were utilized to construct a Cox proportional hazards model. The Akaike information criterion was adopted to select the appropriate model. The risk score for each patient (P_RS) was calculated with the following formula:

Thereafter, the cut-off value for the risk score was determined using the “survminer” package and patients were separated into high- and low-risk groups. The optimal cut-off values for the training, test and entire cohort were 2.026, 1.998 and 2.038, respectively. Kaplan–Meier (K–M) survival curves was plotted using the “survival” package. Time-dependent receiver operating characteristic (ROC) curves was generated using the “timeROC” package to assess the predictive accuracy of the model. In addition, clinicopathological differences between the high- and low-risk groups were analyzed using the Wilcoxon rank sum test.

To validate the CCRLSig as an independent prognostic indicator for PFI in PTC, the signature risk score and clinical parameters (age, gender, tumor size, T stage, N stage, AJCC stage, multifocality, aggressiveness, anatomic site of the tumor, and BRAFV600E mutation status) were subjected to univariate and multivariate Cox regression analyses.

A nomogram incorporating clinical and pathological factors and the risk score was established using the “survival” and “rms” R packages to predict 1-, 3- and 5-year PFI in the entire cohort. The calibration curves and ROC curves were used to estimate the predictive accuracy of the nomogram.

DEGs between the high- and low-risk groups based on the CCRLSig were identified using DEseq2 with a |log FC| > 0.5 and p.adjust< 0.05 as thresholds, after which Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and gene set enrichment analysis (GSEA) of these DEGs were conducted using the R “clusterProfiler” package.

The somatic mutation data of PTC patients were acquired from TCGA. The TMB was compared between the high- and low-risk groups using the “maftools” package. The TMB was measured as follows: (total mutations/total covered bases) × 10^6 for each patient with PTC.

To further explore if the risk signature was associated with immune cell infiltration in PTC, immunity-associated data were extracted from xCell (http://xCell.ucsf.edu/) and parameters were compared between the high- and low-risk groups by the Wilcoxon rank sum test. In addition, correlations between the CCR-lncRNAs and immune cell fractions were evaluated.

The human normal thyroid cell line (Nthy-ori3–1) was provided by Professor Hongmei Shen (Center for Endemic Disease Control, Chinese Center for Disease Control and Prevention, Harbin Medical University, Harbin, China). The human PTC cell line (TPC-1) was obtained from the School of Public Health, Shandong First Medical University & Shandong Academy of Medical Sciences. The cell lines were certified by short tandem repeat (STR) validation. TPC-1 and Nthy-ori3–1 cells were cultured in RPMI 1640 supplemented with 10% fetal bovine serum (Moregate, Bulimba, Australia) and 1% penicillin/streptomycin solution (Beyotime Biotechnology, Shanghai, China) at 37°C and 5% CO₂ conditions.

Total RNA was extracted from the cell lines using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA), and reverse transcribed into cDNA using the Roche Reverse Transcription Kit (Roche, Basel, Switzerland). QRT-PCR was performed using a real-time PCR instrument (Q5). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) served as a control. The primers for lncRNAs were provided by Sangon Biotech (Shanghai, China) and are listed in Table S1. The relative expression level of each lncRNA was calculated by the 2^-△△CT method.

The R statistical environment (V4.1.2) and GraphPad Prism 8.0 were used for statistical analyses. Analysis of variance, Wilcoxon rank sum test, and Student’s t-tests were applied to evaluate differences between groups. The K–M analysis and log-rank test were adopted for comparing PFIs between groups. Values of p< 0.05 were considered statistically significant difference for all analyses.

We identified 655 DElncRNAs and 266 CCRGs in a comparison between 502 PTC and 58 adjacent normal samples (|log FC| > 1, p.adjust< 0.05). Based on Pearson correlation coefficients between the CCRGs and lncRNAs (|R| > 0.4, p< 0.05), we obtained 446 CCR-lncRNAs.

The 498 PTC samples with complete clinical information were divided into training and test cohorts at a ratio of 1:1. In the training cohort, 31 CCR-lncRNAs associated with PFI were selected by univariate Cox analyses. Subsequently, we applied LASSO-Cox regression analysis and established a prognostic signature consisting of five CCR-lncRNAs (FOXD2-AS1, LOC100507156, BSG-AS1, EGOT, and TMEM105; Figures 1A, B ). The K–M survival curves for the five CCR-lncRNAs relating to the PFI are shown in Supplementary Figures S1A–E . Next, the risk score for each patient in the training cohort was calculated using the following formula: risk score = [(0.84919) × normalized expression value of FOXD2-AS1] – [(3.86735) × normalized expression value of LOC100507156] – [(0.67474) × normalized expression value of BSG-AS1] – [(0.76288) × normalized expression value of EGOT] – [(0.46501) × normalized expression value of TMEM105]. Next, according to the optimal cut-off value of the risk score, patients were grouped into high- and low-risk groups. As the risk score increased, patients’ prognosis worsened ( Figures 2A, B ). The expression levels of TMEM105, EGOT, FOXD2-AS1, and BSG-AS1 in the CCRLSig increased and the expression level of LOC100507156 decreased as the risk score increased ( Figure 2C ). K–M analysis showed that PTC patients in the high-risk group had a shorter PFI than that of patients in the low-risk group ( Figure 2D ). The area under the ROC curve (AUC) values for the risk score for 1-, 3-, and 5-year PFI were 0.784, 0.722, and 0.681, respectively ( Figure 2E ).

To validate the accuracy of the signature, patients in the test cohort and the entire cohort were separated into high- and low-risk groups according to the optimal cut-off value for each dataset. The risk curves, PFI status, and heatmaps of risk lncRNA expression profiles in the test cohort and entire cohort were consistent with those of the training cohort ( Figures 2F–H, K–M ). Similarly, a K–M curve analysis of the test cohort and entire cohort indicated that the high-risk groups had a shorter PFI than the low-risk groups ( Figures 2I, N ). In the test cohort, the AUC values for the risk score for 1-, 3-, and 5-year PFI were 0.723, 0.677, and 0.668, separately ( Figure 2J ). In the entire cohort, the AUC values for the risk score for 1-, 3-, and 5-year PFI were 0.757, 0.693, and 0.673, separately ( Figure 2O ). These results suggest that the five-CCR-lncRNA risk signature has good predictive performance for PFI in PTC.

Univariate and multivariate Cox regression analyses were employed to explore whether the risk score based on the CCRLSig predicts PFI in PTC. Univariate Cox analyses showed that the risk score was notably correlated to the PFI in the training, test, and entire cohorts ( Figures 3A–C ). A multivariate Cox analysis proved that the risk score was an independent predictor for the PFI of PTC in the training, test, and entire cohorts ( Figures 3D–F ).

In correlation analyses, the risk scores based on the CCRLSig were higher for patients older than 55 years than for patients younger than 55 years ( Figure 4B ). The risk scores of tumor size >2 cm were higher than those for tumor size ≤ 2 ( Figure 4C ). The risk scores for N1 were higher than those for N0 ( Figure 4D ). We also found that the risk score tended to increase with T stage and AJCC stage ( Figures 4E, F ), implying the critical role of the signature in the progression of PTC. Additionally, the risk score was significantly higher in patients with the BRAFV600E mutation than in those without the mutation ( Figure 4H ). Nevertheless, there were no significant correlations between the risk score and gender or M stage (p > 0.05, Figures 4A, G ).

In the entire cohort, we built a nomogram for the prediction of the 1-, 3-, and 5-year PFI of PTC based on the signature’s risk score and significant clinicopathological features (P< 0.05) identified in univariate Cox regression analyses, including the pathological T stage, N stage, AJCC stage, tumor size, age, and aggressiveness ( Figure 5A ). The calibration curve indicated that the nomogram-predicted PFI at one, three, and five years was highly consistent with the practically observed PFI ( Figures 5B–D ). Furthermore, the AUCs of the nomogram for evaluation of 1-, 3-, and 5-year PFI were 0.796, 0.711, and 0.681, respectively, and the predictive performances were superior to those of other clinical characteristics (age and N stage; Figures 5E–G ). These results suggest that the nomogram reliably predicts the 1-, 3-, and 5-year PFI in PTC.

To investigate the molecular mechanisms and pathways by which the signature is related to the risk of PTC progression, we carried out GO and KEGG enrichment analyses and GSEA of DEGs between the two risk groups. The GO enrichment analysis demonstrated that the DEGs were enriched in multiple biological processes and molecular functions, including cellular calcium ion homeostasis, positive regulation of MAPK cascade, cell-substrate adhesion, positive regulation of cell-cell adhesion, humoral immune response, T cell mediated immunity, and chemokine activity ( Figure 6A ). The KEGG analysis showed that the DEGs were involved in the PI3K-Akt, MAPK signaling pathway, cytokine-cytokine receptor interaction, cell adhesion molecules, antigen processing and presentation, and IL-17 signaling pathway ( Figure 6B ). The GSEA revealed that the DEGs were enriched in the cell cycle pathway and several immune-related biological processes ( Figure 6C ), including the cell cycle, P53 signaling pathway, cell cycle checkpoints, innate immune system, antigen response, and MHC class II antigen presentation.

To reveal genetic variation in risk score subtypes, we compared the TMB between the high- and low-risk groups. Compared with the low-risk group, the high-risk group had a markedly higher TMB (p< 0.01, Supplementary Figure S2A ). The top 20 mutated genes in two risk group are shown in Supplementary Figures S2B, C . We observed that the mutation rate of BRAF was markedly higher in the high-risk group (84%) than in the low-risk group (56%). Supplementary Figures S2D, E show a complete view of the somatic mutations in the high- and low-risk groups.

The functional enrichment analysis displayed that the CCRLSig may be associated with immunity. Hence, we further analyzed the relationships between the signature and immune cell infiltration. The relative frequencies of infiltrating immune cells in all PTC patients are shown in Figure 7A . The high-risk group exhibited higher immune scores than the low-risk group ( Figure 7B ). The fractions of B-cells, CD4+ memory T cells, class-switched memory B-cells, DC, macrophages, NKT, Th2 cells, and Tregs in the high-risk group were markedly higher than those in the low-risk group (p< 0.05). In contrast, the fractions of CD4+ Tcm, CD8+ T cells, and CD8+ Tcm cells in the high-risk group were lower than those in the low-risk group ( Figure 7C ). Additionally, we explored correlations between the expression levels of the five lncRNAs in the signature and the infiltration of multiple immune cells in PTC ( Figure 7D ). These results implied that the signature is linked to immune cell infiltration and may regulate immune processes in PTC.

We analyzed the differential expression of the five lncRNAs between normal and PTC tissues in TCGA data, as illustrated in Figure 8A , and thereafter verified the results in cell lines. The expression levels of FOXD2-AS1, LOC100507156, BSG-AS1, EGOT, and TMEM105 were notably higher in TPC-1 cells than in Nthy-ori3–1 cells ( Figures 8B–F ), which was in line with the results of the bioinformatics analysis, thereby supporting the accuracy of our analysis.