Genomic Instability-Related LncRNA Signature Predicts the Prognosis and Highlights LINC01614 Is a Tumor Microenvironment-Related Oncogenic lncRNA of Papillary Thyroid Carcinoma

Background Genomic instability (GI) is among the top ten characteristics of malignancy. Long non-coding RNAs (lncRNAs) are promising cancer biomarkers that are reportedly involved in GI. So far, the clinical value of GI-related lncRNAs (GIlncs) in papillary thyroid cancer (PTC) has not been clarified. Methods Integrative analysis of lncRNA expression and somatic mutation profiles was performed to identify GIlncs. Analysis of differentially expressed lncRNAs in the group with high- and low- cumulative number of somatic mutations revealed significant GIlncs in PTC. Univariate and multivariate Cox proportional hazard regression analyses were performed to identify hub-GIlncs. Results A computational model based on four lncRNAs (FOXD2-AS1, LINC01614, AC073257.2, and AC005082.1) was identified as a quantitative index using an in-silicon discovery cohort. GILS score was significantly associated with poor prognosis, as validated in the TCGA dataset and further tested in our local RNA-Seq cohort. Moreover, a combination of clinical characteristics and the composite GILS-clinical prognostic nomogram demonstrates satisfactory discrimination and calibration. Furthermore, the GILS score and FOXD2-AS1, LINC01614, AC073257.2, and AC005082.1 were also associated with driver mutations and multiple clinical-pathological variables, respectively. Moreover, RNA-Seq confirmed the expression patterns of FOXD2-AS1, LINC01614, AC073257.2, and AC005082.1 in PTC and normal thyroid tissues. Biological experiments demonstrated that downregulated or overexpressed LINC01614 affect PTC cell proliferation, migration, and invasion in vitro. Activation of the stromal and immune cell infiltration was also observed in the high LINC01614 group in the PTC microenvironment. Conclusion In summary, we identified a signature for clinical outcome prediction in PTC comprising four lncRNAs associated with GI. A better understanding of the GI providing an alternative evaluation of the progression risk of PTC. Our study also demonstrated LINC01614 as a novel oncogenic lncRNA and verified its phenotype in PTC.


INTRODUCTION
Papillary thyroid cancer (PTC) is the most common type of thyroid cancer, with heterogeneous biological behavior and a favorable prognosis (1)(2)(3). In a retrospective study of PTC with 27 years of median follow-up times, the recurrence rate and PTC-specific death rate were 28% and 9%, respectively (4). Although the overall prognosis of most PTC patients was satisfactory, there was still a small part of PTC patients having aggressive characteristics, and even after standard surgical treatment, they are still prone to recurrence and metastasis (5).
Genomic instability (GI) endowed tumors with an inherent survival advantage. As a hallmark of cancer, GI triggers selfsufficiency in evasion of programmed cell death, limitless replicative potential, sustained angiogenesis, and tissue invasion (6,7). PTC exhibits frequent alterations in oncogenes (such as BRAF and RAS), DNA repair defects, and GI (8). Emerging studies demonstrated that GI was highly regulated through DNA damage checkpoints, DNA repair mechanisms, and mitotic checkpoints. Besides, aberrant transcriptional regulation and epigenetic modification are also implicated in genome instability. The exploration of a novel biomarker utilizing the GI signature of PTC might be helpful for risk stratification and prognostic evaluation of PTC patients.
Increasing evidence has shown that many epigenetic regulators in tumors reside in non-coding regions, mostly transcribed into long non-coding RNAs (lncRNAs), which transcripts of more than 200 nucleotides, involved in the survival, proliferation, migration, and genomic stability of cells (9,10). LncRNA NORAD, CUPID1, CUPID2, and DDSR1 promote genomic stability by regulating the expression of DNA repair-associated genes or by interacting with damagerelated proteins (11)(12)(13). A recent study revealed that BGL3 lncRNA mediates BRCA1/BARD1 retention at double-strand breaks (14). lncRNAs act as vital epigenetic regulators of multiple biological processes by binding to different DNA, mRNA, or proteins and are requisite in regulating and maintaining genomic stability and tumorigenesis (15). These mechanisms by which lncRNAs regulate DNA repair-associated genes facilitate our understanding of the link between lncRNA and GI. Moreover, several lncRNA signatures have been established for predicting the prognosis of cancer patients (16)(17)(18)(19)(20), whereas the potential biological and clinical significance of GI-associated lncRNAs in PTC remains largely unknown.
In the following research, we constructed a lncRNA signature from the genomic and transcriptional levels to predict the clinical outcome of PTC. Besides, we also investigated the molecular changes related to PTC associated with the GILS score. GILS predicted the function of LINC01614 in mediating cell proliferation and migration, which was further validated by overexpression and knockdown experiments in PTC cell lines. In summary, our findings showed that the GILS robustly predict patient prognosis and revealed the oncogenic functions of LINC01614, which have great potentials in the future development of PTC biomarkers.

Public Dataset Source and Preprocessing
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). Wholetranscriptome 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.

Clinical Specimens and RNA Sequencing
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.

Establishment of GILS
To identify GIlncs, we examined differentially expressed genes (DEGs) in PTC from TCGA by Wilcoxon rank-sum test in "limma" package (|log 2 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 GIassociated 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 .

Functional Enrichment Analysis
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.

TME Analysis
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.

Cellular Proliferation, Colony Formation, Migration, and Invasion Assay
After siRNA transfection, cells were trypsinized, resuspended, seeded in a 96-well plate with a density of 1.5×10 3 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 3 cells/well) and cultured at 37°C in 5% CO 2 . 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.0mM 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.

Statistical Analysis
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 timedependent 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).

Identification of GIlncs in Papillary Thyroid Cancer Patients
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 coexpressed 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 celladhesion 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.

Development of the GILS for Outcome Prediction
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).
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

Independent Validation of GILS in the PTC Dataset
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 highrisk 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).

Construction and Performance Evaluation of the Nomogram
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.

Further Testing of the GILS in the WMU-PTC Cohort
To an outward promotion of our constructed GILS and GILSnomogram, 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).

Independent Prognostic Significance of GILS and Other Clinical Factors
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).

Association Between GILS and Molecular Features of PTC
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 =  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.

GILS Associated With Clinical Parameters of PTC Patients
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.

Effect of LINC01614 on PTC Cell Lines In Vitro
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 LINC01614plasmid ( 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.

Effect of LINC01614 on PTC Microenvironment
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 nonimmune 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).

DISCUSSION
Several clinicopathological features have been reported to be significantly associated with poor prognosis in PTC patients, including larger primary tumor diameter, older age at first diagnosis, extrathyroidal invasion, lymph node metastasis, and distant metastasis (26). TNM staging and other scoring systems have been used to predict the clinical outcome of PTC patients (27,28). However, the current acknowledged mortality riskstratification system for PTC, which integrated numerous clinicopathological parameters, seemed uninformative for the foresight of recurrence risk (29)(30)(31). Besides, patients with similar clinical phenotypes do not have identical prognoses, suggesting that the present PTC prognostic evaluation system may not provide an accurate clinical prognosis for every patient. Therefore, the PTC prognostic evaluation system needs improvement, and the accuracy of PTC prognoses could be significantly improved using molecular biomarkers. Over the past decade, accumulating researches have demonstrated that lncRNAs played essential roles in multiple biological events, such as maintaining genomic stability (32). Cancer patients with a higher GI level have a worse clinical prognosis, so GI could be used to evaluate the clinical prognosis of tumor patients (33). In recent years, lncRNAs have also been shown to be involved in genome stability, and no relevant studies have investigated the lncRNA signatures of GI in PTC. For example, aberrant expression of lncRNAs could affect the development and progression of malignant tumors (34). Therefore, we constructed a novel signature of GIlncs in PTC and explored their significance in predicting the patient's prognosis. Numerous genomic classifiers have been used to predict the clinical outcome of PTC. For example, the ThyroSeq v3 classifier integrated the mutational information of 112 genes, a collection of gene fusions, and the expression profile of 19 DEGs (35). As for non-coding genes, miRNAs were proposed as the regulators of the essential phenotypes of PTC by fine-tuning gene expression (signaling, differentiation, and invasion process) (1, 36). By comparison, few classifiers shed light on the significance of lncRNAs in the prognosis of PTC. Unlike miRNAs that acted on major post-transcriptional mechanisms, lncRNAs were able to combine with homologous nucleotide sequences (37). According to our results of functional enrichment analysis, it could be inferred that lncRNA signatures were involved in the cancer hallmarks as described above, which provided sufficient support for the predictive role of GI. Additionally, the complexity of spatial folding configurations endowed lncRNAs with the binding sites, which could combine certain macromolecular proteins. lncRNAs with aberrant expression may dampen the cancer genome's stability by exerting coupling effects to disrupt the regulation of lncRNA-related PCGs. Significantly differentially expressed lncRNAs are thus, profiled as predictors of genome instability.
Recently, an increasing number of studies have revealed the potential links between lncRNA and GI. It was noteworthy that over the past decades, increasing studies of lncRNA in the prognosis of human cancers were published through the organic combination with GI (14,(38)(39)(40)(41). For instance, analogous analysis in breast cancer illustrated a novel panel of lncRNAs as GI-associated tumor biomarkers (16). The GILS model could stratify risk subgroups based on the degree of GI and serve as a robust prognostic index independently. With the development of the function of lncRNA binding sites and their role in genomic variation, increasing researchers have paid attention to the mechanisms of lncRNA in cancer. Up to now, it has been proved that certain lncRNAs were involved in tumorigenesis by promoting DNA damage, altering drug metabolism, regulating cell apoptosis, affecting the EMT process (6,14,(42)(43)(44).
Currently, there is no effective biomarker to predict the prognosis of PTC. We constructed a GILS model as an independent prognostic predictor based on multivariate Cox regression analysis and exhibited the prediction model using a visual nomogram. In our nomogram, ETE and tumor stage were utilized as binary variables, while the GILS score was utilized as a continuous variable. Other variables significantly related to the clinical prognosis of PTC will also be included in the future prediction model. To assess the clinical significance of this competing-risk nomogram, we performed the net benefits analysis using DCA. As a new analysis for evaluating the clinical prediction model, DCA could integrate the preferences of patients or decision-makers into the analysis (45). By calculating the net benefits at different threshold probabilities, we found that the GILS-nomogram had a satisfactory net benefit, which potentially implied clinical practice.
The current study investigated the role of lncRNAs in GI and identified four lncRNAs as independent prognostic biomarkers for PTC. First, we selected 558 lncRNAs differentially expressed in GS-like and GU-like phenotypes, which were enriched in GIrelated pathways based on the functional enrichment analysis with the lncRNA-derived mRNAs. Next, the prognostic significance of GI-derived lncRNA biomarkers was identified after the univariate and multivariate Cox regression analyses.
Further internal validation emphasized the general applicability of these biomarkers. Many studies have demonstrated that some GILS-derived lncRNAs were associated with specific cancers. In other cases, FOXD2-AS1 had been identified as an oncogene promoting the progression of different cancers, such as PTC and bladder cancer (46)(47)(48). FOXD2-AS1 also enhanced chemotherapeutic resistance in esophageal and laryngeal squamous cell carcinoma (49,50). FOXD2-AS1 was recognized as a molecular biomarker and participated in chemo-resistance, malignant proliferation, invasion, and immunosuppression (46,(51)(52)(53)(54)(55). Recently, a study demonstrated that overexpression of FOXD2-AS1 was associated with poor clinical outcomes in PTC (E, F) Colony formation and CCK-8 assays were performed to test the survival of PTC cells after silencing or overexpressing of LINC01614. (G) Transwell assays were carried out to test the migration and invasive activity of LINC01614 altered PTC cells. The results above were summarized as bar graph. All data are the means ± SD of three experiments. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
patients, and the knockdown of FOXD2-AS1 could suppress tumor growth by reversing its sponge effect (56). However, AC073257.2 and AC005082.1 have not been reported in cancer yet. The mechanisms of their function in cancer require further research. Furthermore, we observed that GILS could identify the BRAF mutation status and share commonalities with clinical prognostic parameters. In addition, lncRNA, mRNA & miRNA likely work together to form a competing endogenous RNA (ceRNA) regulatory network, which could be used as prognostic markers in cancers (57)(58)(59). GI-associated ceRNA regulatory network should be evaluated in the future studies. A pan-cancer analysis only briefly demonstrated that the LINC01614 was upregulated in different cancers, but no further in-depth analysis was performed in PTC (60). Another research found that LINC01614 could serve as an emulative sponge to combine with miR-383 and result in aggressive behavior in glioma (61). In esophageal squamous cell carcinoma and breast cancer, LINC01614 was also suggested as a survival predictor and had an impact on tumor invasion (62,63). However, the functions of LINC01614 are not clear in PTC. We found that LINC01614 is significantly overexpressed in PTC compared with the non-tumor thyroid tissues in both TCGA and local cohorts. In experiments, silencing of LINC01614 significantly inhibited survival of PTC cells while over-expressing of LINC01614 reverse it. In addition, Transwell analysis indicated that LINC01614 promotes PTC cells migration and invasion abilities. Both immune cell and stromal cell exerts huge to regulate tumor progression, and TME activity is closely associated with genomic instability (64,65). ESTIMATE and xCell analysis suggested LINC01614 closely associated with infiltration level of stromal and epithelial cells. TIP analysis demonstrated LINC01614 closely associated with cancer-immunity cycle in PTC.
This study has several limitations. First, the number of patients in our local testing cohort is still small, and the median follow-up period wasn't long enough. Thus, additional larger cohorts are required to validate these findings. Second, although we investigated the potential value of GILS using bioinformatics analysis and using local testing cohort, further testing is still lacking in external datasets. Third, the potential ceRNA-regulated mechanism of LINC01614 in PTC progression needs further investigation. Lastly, we only verified that LINC01614 is an oncogene, but its oncogenic function whether related to genomic instability should be explored.
In conclusion, we screened out GIlncs and constructed a corresponding lncRNA signature. Four GIlncs were identified as independent prognosis factors. Based on this, we established a composite GILS-nomogram for PTC patients to predict the clinical outcome, which was verified in internal validation and our local cohort. This GILS model could identify high-risk PTC patients and precisely formulate treatment strategies, which could potentially benefit the management of PTC patients. These findings not only suggested that the GILS might serve as the biomarker for pathological classification system and prognosis prediction but demonstrated LINC01614 as a novel oncogenic TME-related lncRNA in PTC. Comparison of the activity of the cancer immunity cycle between the LINC01614 high and LINC01614 low groups. ns, no-significance, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because the local RNA-sequencing data used in this paper belongs to a larger unpublished clinical research project. According to the terms of the project contract, the complete sequencing data will be uploaded to the GEO when the main articles are published. This paper is a part of the project: Major Science and Technology Projects of Zhejiang Province, 2015C03052. Requests to access the datasets should be directed to the corresponding author(s).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by All research protocols have been approved and implemented through the ethical standards of the institutional review board of the First Affiliated Hospital of Wenzhou Medical University (Approval No. 2012-57). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
XD contributed to study design, bioinformatic analysis, and manuscript editing. XD, YC, and Z-qY contributed to the manuscript draft. DC and CJ collected and analyzed clinical samples for the study. XH, XZ, WZ, and D-nG discussed the results and participated in the critical review of the manuscript.