A total of 568 thyroid cancer RNA-Seq profile samples in TCGA, including 58 matched normal samples, 502 thyroid cancer samples, and eight metastatic thyroid carcinoma samples, were downloaded using the “TCGAbiolinks” package (21). Whole-transcriptome sequencing data was performed using FPKM expression level in transcripts per million (TPM). The definition and outcome of progression-free survival (PFS) were obtained from the TCGA-Clinical Data Resource (CDR) (22). Mutation status was obtained from Mutation Annotation Format files (derived from MuTect2) from the Genomic Data Commons portal. After excluding two patients with follicular thyroid cancer and two patients who received pretreatment, the expression data, prognostic information, and mutation profiles were included in 487 patients in this study. To exclude genes with high variability across patients, we calculated the median absolute deviation (MAD) of the 487 samples. LncRNAs with MAD > 0.5 were defined as genes with high variability and excluded in the RNA-Seq matrix.

Seventy-nine pairs of thyroid tissues were obtained from the Department of Thyroid Surgery, The First Affiliated Hospital of Wenzhou Medical University. Fresh tissues were immediately snap-frozen in liquid nitrogen and stored at -80°C until further use. All pathological reports were independently confirmed by two experienced pathologists.

The RNA-seq experimental protocol was performed as described in a previous publication (23). Briefly, total RNA was used to construct cDNA libraries through high-throughput RNA sequencing. The RNA expression profile was determined from the sequencing libraries generated from a NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, United States). The clustering of samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina), and the library was sequenced on an Illumina NovaSeq platform. The read counts were normalized to TPM, and the TPM expression values were further log2 transformed.

To identify GIlncs, we examined differentially expressed genes (DEGs) in PTC from TCGA by Wilcoxon rank-sum test in “limma” package (|log₂foldchange| > 0.5, false discovery rate (FDR) < 0.05). To assess the genomic instability, we proposed a mutator hypothesis-derived calculation method: we determined the cumulative number of somatic mutations based on the number of changed sites for each gene on each sample and categorized the patients in descending order. The top 20% of patients were titled with genomic unstable-like (GU-like) group and the last 20% as genomic stable-like (GS-like) group. Then, the GILS for prognosis prediction was developed based on the coefficient of each prognostic GIlncs in the model and their expression levels.

To evaluate the performance of GIlncs on the prognosis, GIlncs were selected by Univariate and multivariate Cox proportional hazard regression analysis. GIlncs with p < 0.05 in univariate Cox were retained, and multivariate Cox was performed in the model group by the “glmnet” package. Based on the coefficients from the multivariate regression analysis and the expression levels of prognostic GIlncs, we constructed a GI-associated lncRNA signature (GILS) for prognostic prediction as follows:

Where N is the number of prognostic GILS, expression (lncRNA) _i is the expression value of prognostic GILS _i , and Coefficient (lncRNA) _i is the estimated multivariable Cox regression coefficient of GILS _i .

Function “cor. test” in R was used to measure the Spearman correlation coefficients between GI-lncRNAs and mRNAs, and the top 10 mRNAs were selected as PTC-specific lncRNA targets. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed using “enrichGo” and “enrichKEGG” functions in the “clusterprofiler” package, respectively.

ESTIMATE is a ssGSEA based algorism for tumor-stroma purity detection, which uses the gene expression profiles of 141 immune and stromal genes. ESTIMATE score is a combination of immune and stromal scores, calculated by the “ESTIMATE” R package. The xCell project (24) used six public cells sorted bulk gene expression data sets to generate gene signatures and score each TCGA sample.

Total RNA was extracted using TRIzol Reagent (Invitrogen). Complementary DNA (cDNA) was synthesized using ReverTra Ace qPCR RT Kit (Toyobo, Japan). Real-time qPCR was performed using SYBR Select Master Mix on 7500 Fast Instrument. Each sample was tested in triplicate. Experiments were repeated three times. Human primer sequences for qPCR are as below: LINC01614, forward 5’-TCAACCAAGAGCGAAGCCAA-3’, reverse- 5'-TTGGACACAGACCCTAGCAC-3'; GAPDH, forward 5’-GTCTCCTCTGACTTCAACAGCG-3’, reverse- 5'-ACCACCCTGTTGCTGTAGCCAA-3’.

Full-length LINC01614 cDNA was synthesized and inserted into pCDH-GFP+PURO-3xFlag and pCDH-GFP+PURO vectors (Genepharma, Shanghai, China). The resulting vector or empty vector was transfected into PTC cells using Lipofectamine 2000 Transfection Reagent (Life Technologies, Carlsbad, CA) according to the manufacturer’s protocol. Infected cells were selected with puromycin (Invivogen) at 1 μg/ml. For LINC01614 knockdown, siRNA transfections were carried out using Lipofectamine RNAiMAX Reagent (Invitrogen). The siRNAs containing the following two individual siRNAs were used: siLINC01614#1: forward 5’-GCCCACCTCAAATCCTGAA-3’; siLINC01614#2: forward 5’-GCUGGAAGCAUUUCGUAAU-3’.

After siRNA transfection, cells were trypsinized, resuspended, seeded in a 96-well plate with a density of 1.5×10³ cells/well, and incubated at 37°C. At each indicated time-point, 10 µl of CCK-8 was added and incubation was continued for 2 h. The plates were agitated and the absorbance was measured at 450 nm under an absorption spectrophotometer. For colony formation assay, cells were seeded in 6-well plates (1×10³ cells/well) and cultured at 37°C in 5% CO₂. After 10 days, the cells were washed with PBS and stained with crystal violet. For migration assay, Cells were plated at a known density in the upper chamber of 8.0μM membrane transwells (Cat. #3422, Corning, Tewksbury, MA) in media containing 10% FBS. Transwells were placed in wells with media containing 10% FBS. Cells were fixed onto the transwell membrane in 10% formalin. Migrated or invaded cells were imaged in a 10× magnification microscope in 5-10 random fields for each well and quantified by ImageJ software.

Mann–Whitney test or Wilcoxon signed-rank test were used in the two-group analysis. Comparisons between different groups were conducted by Kruskal–Wallis or one-way ANOVA. Survival analysis was conducted using the Kaplan–Meier method and the log-rank test. Univariate and multivariate analyses with Cox proportional hazards regression for PFS were performed on the individual variables by calculating the hazard ratios (HR) and 95% confidence intervals (CI). The time-dependent ROC curve was calculated with the nearest neighbor estimation method. In all experiments, three biological replicates were performed for each group. All statistical analyses were performed using the R software (3.5.2).

A flowchart is shown in Figure 1 . The cumulative number of somatic mutations in each patient was calculated and ranked in descending order. In a mutation cumulation-based method, the top 20% (n = 104) and the last 20% (n = 109) patients were assigned into the GU-like group and GS-like group, respectively. We identified 558 GIlncs were differentially expressed between the two groups ( Supplementary Table S1 ), among which 533 upregulated and 25 downregulated lncRNAs in the GU-like group ( Figure 2A ).

To investigate the potential functions of the GIlncs, we explored the protein-coding genes (PCGs), which were co-expressed with lncRNA. Based on the consistency of expression, the GI-related PCGs profile was finally constitutive of the 2598 PCGs related to the lncRNA profile ( Supplementary Table S2 ). The lncRNA-mRNA co-expression network was displayed ( Figure 2B ). To figure out the coupling effect between lncRNA and GI, we concentrated on the underlying functions and related pathways of these lncRNAs. Functional enrichment analysis of the GI-related PCGs was performed. KEGG pathway ( Figure 2C and Supplementary Table S3 ) analysis revealed that significant pathways played a role in transcriptional misregulation in cancer, deregulation of cellular energy within glycolipid and choline metabolism, and endocrine resistance. Besides, enrichment of canonical pathways, such as PI3K-AKT, MAPK, and JAK-STAT pathway, provided sustained proliferative signaling. EGFR tyrosine kinase inhibitor resistance caused apoptosis inhibition GO analysis listed ( Figure 2D and Supplementary Table S4 ) covered three aspects: biological process (BP), cellular component (CC), and molecular function (MF) in detail (25). In BP terms, the tumor microenvironment components, including leukocyte migration, T cell activation, and lymphocyte differentiation, indicated their involvement in immuno-inflammatory responses. In CC terms, migration and invasion abilities of the tumor cells were also implicated, which was shown to correlate with the extracellular matrix and cell-adhesion membrane. MF analysis stressed the significance of metabolic components such as nucleotides and polysaccharides. The above results elaborated that GIlncs lighted a novel perspective on the mechanism of PTC progression.

The whole TCGA cohort was randomly divided into two datasets for discovery and validation. The clinical factors were not significantly different in both groups (all p > 0.05, Supplementary Table S6 ). The discovery set (n = 249) was used to screen clinically significant lncRNA and construct a corresponding lncRNA signature. The validation set (n = 248) was used to determine the accuracy of the signature. Univariate Cox proportional hazards regression analysis revealed that 30 lncRNAs were derived from differential lncRNAs (p < 0.05; Figure 3A ). Specifically, three lncRNAs were indicated unfavorable outcomes with the HR > 1, whereas 27 lncRNAs with the HR < 1. Next, multivariate Cox proportional hazards regression analysis revealed four lncRNAs (FOXD2-AS1, LINC01614, AC073257.2, and AC005082.1) as independent prognosis biomarkers ( Supplementary Table S7 ).

Based on the coefficients of multivariate analysis and lncRNAs expression, GILS was defined as follows: GILS score = 0.304 * FOXD2-AS1 + 0.036 * LINC01614 + (-0.769) * AC073257.2 + (-0.328) * AC005082.1. Of the GILS, the coefficient of lncRNA FOXD2-AS1 and LINC01614 was positive, suggesting that they are risk factors, while the other lncRNA AC073257.2 and AC005082.1 tended to be protective factors. A comparison of individual lncRNA signature with the degree of GI demonstrated that high somatic mutation was accompanied by upregulation of FOXD2-AS1 and LINC01614 and downregulation of AC073257.2 and AC005082.1 ( Figure 3B ).

Patients in the discovery set were labeled as high- or low-risk based on the calculated GILS score, which utilized the median score as the threshold. The overall mutation rate was significantly higher in the high-risk group compared with the low-risk group (p = 8.2e−07; Figure 3C ). Besides, Kaplan–Meier analysis of the discovery set showed that there was a significant difference in PFS between high and low-risk groups (p < 0.0001, log-rank test; Figure 3D ). As shown in Figure 3E , the GILS showed an accurate estimation performance with an area under the curve (AUC) of 0.798, 0.772, and 0.771 at 1, 3, and 5 years, respectively.

To verify the performance of the GILS for prognosis prediction, we calculated the GILS scores of the validation set and the entire TCGA set. Patients in the low-risk group showed a more prolonged survival than patients in the high-risk group with an AUC value of 0.764 in the validation cohort (p < 0.0001; Figures 4A, B ). Similar results were also observed in the entire TCGA set, where the AUC of the ROC curves for GILS was 0.751 (p < 0.0001; Figures 4E, F ). These results indicate that GILS has a good survival prediction efficacy. In both validation cohorts, FOXD2-AS1 and LINC01614 were upregulated, while AC073257.2 and AC005082.1 were down-regulated in the high-risk group ( Figures 4C, G ). Comparison analysis showed that there were significant differences in the number of somatic mutations between the high-risk and low-risk groups in both the validation (p < 0.0001; Figure 4D ) and the entire TCGA cohort (p < 0.0001; Figure 4H ).

In the discovery cohort, univariate analysis revealed the pathologic stage, ETE, and GILS score were significantly associated with PFS. To improve the clinical practicability of GILS, we established a statistical GILS-nomogram model in the discovery set by integrating GILS score, ETE, and pathologic stage ( Figure 5A ). The calibration curve showed good agreement among the estimations with the GILS-nomogram and actual observations in both discovery and validation cohorts ( Figures 5B, E ). The GILS-nomogram also yielded supreme concordance index (C-index) when compared to pathologic stage and ETE ( Figures 5C, F ). Decision curve analysis (DCA) revealed that the net benefit of the GILS-nomogram was improved when compared with other prognostic factors alone ( Figures 5D, G ). Therefore, these findings indicate improved prediction performance of the GILS-nomogram.

To an outward promotion of our constructed GILS and GILS-nomogram, we carried out the same procedures for external testing in our cohort (n = 79). Only FOXD-AS1 and LINC01614 achieved satisfactory discrimination in both the high- and low-risk cohorts (p = 0.009 and 0.052, respectively; Figures 6A, B ), while AC073257.2 and AC005082.1 had no significant difference (p = 0.185 and 0.554, respectively; Figures 6C, D ). We calculate the GILS score for each patient in the WMU-PTC cohort ( Supplementary Table S5 , 9 ), and we found high GILS score was associated with shorter PFS (p = 0.085, Figure 6E ). As expected, the net benefit of the GILS-nomogram was highest when compared with other prognostic factors ( Figure 6F ).

To verify the independent predictive effects of GILS, we combined GILS with clinical factors for univariate and multivariate Cox analysis. The results revealed that ETE, disease stage, GILS and was independent prognostic factors ( SupplementaryTable S8 ). Besides, to assess the risk clustering ability of GILS in different risk stratifications, we analyzed prognostic significance of GILS in whole TCGA cohort after adjusted by clinical characteristics. There were significant difference in PFS between the young-patient (p < 0.0001; Supplementary Figure S1A ) and old-patient groups (p < 0.0001; Supplementary Figure S1B ); female-patient (p < 0.0001; Supplementary Figure S1C ) and male-patient groups (p = 0.018; Supplementary Figure S1D ); T1+T2 (p < 0.0001; Figure S1E ) and T3+T4 groups (p = 0.007; Supplementary Figure S1F ); Without-LNM (p = 0.003; Supplementary Figure S1G ) and with-LNM groups (p < 0.0001; Supplementary Figure S1H ); Early-stage (p < 0.0001; Supplementary Figure S1I ) and advanced stage groups (p = 0.009; Figure S1J ); Without-RT (p = 0.043; Supplementary Figure S1K ) and with-RT groups (p < 0.0001; Supplementary Figure S1L ). In summary, these results suggested that the prognostic significance of GILS and might widely applied in different risk stratifi;cations.

We divided the whole TCGA-PTC cohort into high and low-risk groups based on stratification by GILS score. The distinctive mutation distribution profiles of the top mutated genes in the two groups are shown ( Figures 7A, B ). The MAPK and PI3K-AKT pathways play an important role in PTC malignancy. The MAPK and PI3K-AKT pathway-related mutations, BRAF and RAS, were frequently observed in both high and low-risk subgroups. We compared the mutated percent between the high-risk and low-risk groups in train, validation, and whole TCGA datasets using the chi-square test. The results showed that patients in the high-risk group displayed a significantly higher proportion of BRAF mutations than those in the low-risk group among the three datasets (Whole set: p = 0.002; discovery set: p = 0.016; validation set: p = 0.078; Figure 7C ). However, we didn’t find a significant difference of RAS mutations between the two groups ( Figure 7D ). We further explored whether four GIlncs are associated with PTC molecular characteristics. The results revealed that four GIlncs were all upregulated in BVL compared to the RL group (p < 0.05; Figure 7E ), especially for FOXD2-AS1 and LINC01614. We also observed that the FOXD2-AS1 and LINC01614 were significantly upregulated in BRAF ^V600E mutation groups (p < 0.001; Figure 7F ), while FOXD2-AS1 and LINC01614 were significantly downregulated in RAS mutation groups (p < 0.01; Figure 7G ). Interestingly, we found AC073257.2 was upregulated in RET fusion group (p < 0.001; Figure 7H ). As described above, we found both GIlncs and GILS scores were correlated with critical molecular characteristics of PTC.

To investigate four GIlncs and GILS in PTC progression, we evaluated four GIlncs levels and GILS scores in different clinical risk stratifications. We found FOXD2-AS1 and LINC01614 levels and GILS score were positively associated with T stage, N stage, disease stage, the degree of ETE, and histological subtypes. In addition, AC073257.2 was differentially expressed in different histological subtypes, while AC005082.1 was negatively associated with T stage, disease stage, and ETE ( Figures 8A-E ). These results indicated both GIlncs and GILS scores were significantly associated with the clinical outcome of PTC patients.

We compared the expression profiles of four GIlncs in both tumor and normal tissue in TCGA, GTEx, and WMU-PTC cohorts. In the pan-cancer analysis, we found that FOXD2-AS1 and LINC01614 were highly expressed in most cancer types ( Supplementary Figure S2 ). Besides, FOXD2-AS1, LINC01614, and AC073257.2 were upregulated in PTC compared to normal tissue, while AC005082.1 was downregulated in the TCGA cohort ( Figure 9A ). In the local cohort, we validated LINC01614 was upregulated in PTC compared to normal tissue, whereas AC005082.1 was downregulated ( Figure 9B ). To explore the potential effect of LINC01614 in PTC tumorigenesis, we silenced LINC01614 with specific siRNA ( Figure 9C ) and overexpressed LINC01614 with LINC01614-plasmid ( Figure 9D ). LINC01614 knockdown and overexpression significantly decreased and increased colony formation, respectively ( Figure 9E ). CCK-8 assay showed that LINC01614 knockdown and overexpression decreased and increased PTC cell proliferation, respectively ( Figure 9F ). Furthermore, results of transwell assay ( Figure 9G ) showed that LINC01614 depletion or overexpression in PTC cells decreased or increased cellular migration and invasion, respectively. Our results indicated that altered LINC01614 may affect cell growth and migration in the PTC cells.

Furthermore, we aimed to explore the underlying correlation between the LINC01614 mRNA level and the PTC microenvironment. As shown in Figure 10A , LINC01614 was found to be positively correlated with the immune score (r = 0.4, p < 0.001), stromal score (r = 0.58, p < 0.001), and TME score (r = 0.5, p < 0.001), whereas negatively correlated with tumor purity (r = -0.5, p < 0.001). By comparing different types of non-immune TME signatures between LINC01614 ^high and LINC01614 ^low subgroups ( Figure 10B ), we found that the LINC01614 ^high subgroup was marked with higher infiltration level of adipocytes, chondrocytes, mesangial cells, astrocytes, epithelial cells, keratinocytes, and sebocytes (All p-value < 0.001), whereas LINC01614 ^low subgroup was marked with higher infiltration level of endothelial cells, MSC, osteoblast, smooth muscle cells, hepatocytes, and neurons (All p-value < 0.001).

The activities of the cancer immunity cycle are a direct comprehensive performance of the functions of the chemokine system and other immunomodulators. In the LINC01614 ^high group, activities of the majority of the steps in the cycle were found to be upregulated, including priming and activation (Step 3), and trafficking of immune cells to tumors (Step 4) (CD8 T cell recruiting, Macrophage recruiting, Th1 cell recruiting, NK cell recruiting, DC recruiting, and TH17 recruiting) and infiltration of immune cells into tumors (Step 5),. Interestingly, the activity of recognition of cancer cells by T cells (Step 6) was downregulated in the LINC01614 ^high group (All p-value < 0.001, Figure 10C ).