A Potential Nine-lncRNAs Signature Identification and Nomogram Diagnostic Model Establishment for Papillary Thyroid Cancer

Abstract

The purpose of our current study was to establish a long non-coding RNA(lncRNA) signature and assess its prognostic and diagnostic power in papillary thyroid cancer (PTC). LncRNA expression profiles were obtained from the Cancer Genome Atlas (TCGA). The key module and hub lncRNAs related to PTC were determined by weighted gene co-expression network analysis (WGCNA) and LASSO Cox regression analyses, respectively. Functional enrichment analyses, including Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene set enrichment analysis were implemented to analyze the possible biological processes and signaling pathways of hub lncRNAs. Associations between key lncRNA expressions and tumor-infiltrating immune cells were identified using the TIMER website, and proportions of immune cells in high/low risk score groups were compared. Kaplan-Meier Plotter was used to evaluate the prognostic significance of hub genes in PTC. A diagnostic model was conducted with logistic regression analysis, and its diagnostic performance was assessed by calibration/receiver operating characteristic curves and principal component analysis. A nine-lncRNAs signature (SLC12A5-AS1, LINC02028, KIZ-AS1, LINC02019, LINC01877, LINC01444, LINC01176, LINC01290, and LINC00581) was established in PTC, which has significant diagnostic and prognostic power. Functional enrichment analyses elucidated the regulatory mechanism of the nine-lncRNAs signature in the development of PTC. This signature and expressions of nine hub lncRNAs were correlated with the distributions of tumor infiltrating immune cells. A diagnostic nomogram was also established for PTC. By comparing with the published models with less than or equal to nine lncRNAs, our signature showed a preferable performace for prognosis prediction. In conclusion, our present research established an innovative nine-lncRNAs signature and a six-lncRNAs nomogram that might act as a potential indicator for PTC prognosis and diagnosis, which could be conducive to the PTC treatment.

Introduction

Thyroid cancer (TC) is well-known as a common endocrine malignancy, which induces about 40,000 deaths worldwide annually [1, 2]. Papillary thyroid cancer (PTC) is the most frequent subtype of TC, and approximately 90% of TC patients fall under the category of PTC [3]. Recent therapeutic strategies for PTC patients mainly include surgical resection and radioactive iodine therapy [4, 5]. Genetic mutation and environmental exposure are major risk factors of PTC [6]. Despite more than 90% of PTC patients having satisfactory overall survival, a small proportion (about 10%) of patients with PTC have high recurrence and metastasis rates, which induces aggressive diseases and dismal prognosis [7]. Thus, it is urgently needed to screen cancer-specific factors and establish a prognostic nomogram for PTC patients’ prognosis, thereby developing a potential and accurate risk evaluation mechanism.

Non-coding RNA (ncRNA), as a type of RNA, does not code for protein but holds an enzymatic, structural, or regulatory role. Based on its transcript length, it can be classified into either small or long ncRNA. microRNA (miRNA), as the most acknowledged class of short ncRNA, is involved in the specific regulation of its target messenger RNAs (mRNAs) through inhibiting translation or inducing degradation [8]. Long non-coding RNAs (lncRNAs), as a novel class of ncRNAs, comprise more than 200 nucleotides and encode no protein [9]. In spite of having no recognized protein-coding capacity, lncRNAs are proven to have diversely structural and functional impacts in manifold processes. Data from multiple studies has implied that the influences of lncRNAs and their mechanisms of gene expression and regulation may be more extensive and complicated than those of miRNAs [10–12]. LncRNAs are expressed in many tissues in vivo, and their expression levels are much lower than protein coding genes, but their tissue specificity is much higher than that of protein coding genes, which makes them potentially advantageous as biomarkers with highly specific diagnostic effects [13].

Emerging evidence has suggested that lncRNAs are implicated in strings of physiological or/and pathological processes, including cell viability, apoptosis, transcriptional regulation and tumorigenesis, through regulating the gene expression at the post-transcriptional level [14, 15]. Importantly, lncRNAs are considered as oncogenic factors or tumor-suppressors in the development of PTC. For example, lncRNA TUG1 contributes to cell motility in PTC cells via targeting miR-145 [16]. Low expression of lncRNA EMX2OS can be used to predict the recurrence-free survival of PTC [17]. Wang et al. have indicated that up-regulation of lncRNA SLC26A4-AS1 could repress the epithelial-mesenchymal transition (EMT) process through mediating the MAPK signaling in PTC [18]. LncRNA DGCR5 has been identified as a tumor inhibitor in PTC [19].

The Cancer Genome Atlas (TCGA), an open source database, contains large-scale genome sequences that could be used to map genomic curves of all kinds of human cancers [20]. Development of omics technology has contributed to the diagnosis and treatment of cancers in a systematic approach [21]. Recently, a method that commonly utilized microarray data for the express profiling of lncRNAs has been well-constructed, which might advance the therapeutic measures for human cancers [22, 23]. A unique prognostic nomogram has been established based on the five-lncRNAs signature in clear cell renal cell carcinoma [24]. Moreover, a five-gene signature has been identified with a prognostic nomogram for PTC [25]. This kind of prognostic model has been considered as an important tool for risk evaluation in tumors with the consideration of roles of diverse genes on overall survival of cancer patients. However, the optimal biomarkers for PTC prognosis still have not been well characterized.

Here, we attempted to construct a tumor-specific lncRNA signature through analyzing the RNA-seq derived from TCGA database. According to the instruction of Cox regression analysis and WGCNA, a nine-lncRNAs signature with significant prognostic value in PTC was identified. Clinical association of it was confirmed through a Chi-square test. ROC curves, Kaplan-Meier (KM) plotter, and analysis with the help of TIMER website were performed to measure the importance of nine hub lncRNAs in PTC. Through logistic regression, a diagnostic model was established, whose predictive value was determined with calibration/ROC curves and PCA analysis. In short, the present investigation provides a favorable model for predicting PTC prognosis and contributes to elucidating the mechanism of patient poor prognosis.

Materials and Methods

Data Obtained and Analysis

Expression profiles were obtained from TCGA (https://cancergenome.nih.gov/). Patients with primary PTC were selected. Samples with incomplete clinical data and a follow-up time of <30 days were removed and finally the first group of 448 samples was obtained (Table 1 for detailed clinical information), namely the tumor group; and the second group that contained 56 normal samples, was named the normal group. The age and gender of PTC and normal control samples were matched (Supplementary Table S1).

Perl 5.0 (http://www.perl.org/) was applied for background correction and normalization of the lncRNAs expression profiles that we harvested. These downloaded lncRNAs were annotated by using Ensembl (http://www.ensembl.org/index.html), and then the differentially expressed lncRNAs (DElncRNAs) in PTC were identified using the R software “edgeR” package (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) with the criterion of |log fold change (FC)| ≥ 1 and p < 0.05 [26].

WGCNA

WGCNA was employed to build a co-expression network using the 1,601 DElncRNAs that have been screened and 448 samples with complete clinical information (survival time, survival status, age, sex, pathological stage, primary tumor (T) stage, distant metastasis(M) stage, and lymph node metastasis (N) stage, etc.) [27]. After evaluating the quality of the expression matrix, we employed the WGCNA package in R to filter the soft threshold and choose β = 4 for the scale-free network construction. By setting min Module Size as 30 and MEDissThres as 0.25, key modules were determined.

Establishment of lncRNAs Signature in PTC

The LASSO Cox regression model was established to target the genes evidently related with prognosis of PTC patients [28]. DElncRNAs with p < 0.05 were classified into two types, risky lncRNAs (hazard ratio (HR) for death > 1) and protective lncRNAs (HR for death < 1). Due to the linear combination of the expression level of lncRNAs multiplied regression coefficient derived from the univariate Cox regression analysis, the prognostic signature was established by multivariate Cox regression analysis: Risk Score = exp_lncRNA1×beta_lncRNA1 + exp_lncRNA12 × beta_lncRNA2 + exp_lncRNA3 × beta_lncRNA3 + exp_lncRNA4 × beta_lncRNA4 + exp_lncRNA5 × beta_lncRNA5+…+ exp_lncRNAn × beta_lncRNA1n. The “exp” represents the standardized expression of each identified hub lncRNA and the “beta” was determined by multivariate Cox regression analysis.

Assessment of the Nine-lncRNAs Signature

According to the median value of risk score, PTC patients were assigned into high risk (n = 224) and low risk (n = 224) groups. KM method was utilized to predict the prognosis of PTC patients in high- and low-risk groups. The univariate and multivariate Cox regression analyses were performed to detect whether this lncRNAs signature could act as an independent factor for prognostic prediction of PTC patients. A ROC curve was plotted to evaluate the accuracy of the nine-lncRNAs signature in PTC patients. Moreover, a prognostic nomogram of nine-lncRNAs signature was conducted to predict 1-, 2-, and 3-years overall survival rates.

Functional Enrichment Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were implemented to estimate the roles of the co-expressed genes of nine hub lncRNAs. Pearson correlation analyses were implemented for the examination of nine hub lncRNAs’ co-expressed genes based on the TCGA database, and the thresholds were set at the absolute value of correlation coefficient ≥ 0.5 as well as p < 0.05. To clarify the potential mechanisms of the nine-lncRNAs signature, we subsequently performed gene set enrichment analysis (GSEA, http://software.broadinstitute.org/gsea/index.jsp) as the following thresholds: false discovery rate (FDR) < 0.25 & p < 0.05.

Analysis of Gene Expression and Tumor-Infiltrating Immune Cells

CIBERSORT algorithm was utilized to calculate tumor infiltrating immune cells, according to on gene expression predicting the proportions of different immune cells types in the gene expression profiles. The correlation between immune cells and risk score was evaluated by R software “vioplot” package.

TIMER (https://cistrome.shinyapps.io/timer/) was applied to explore the correlation between the hub lncRNAs and tumor infiltrating immune cells (B cell, CD4 T cell, CD8 T cell, Dendritic cell, Macrophage, Neutrophil), relying on 10,897 clinical specimens retrieved from the TCGA database [29].

Construction of Diagnostic Model of hub lncRNAs

The logistic regression analysis was utilized to screen important hub lncRNAs and establish a diagnostic model based on the expression values in normal and PTC samples. ROC and calibration curve (1,000 bootstrap resamples) were utilized to assess the performance of the diagnostic model. PCA was conducted to detect the effectiveness of the diagnostic model with six important lncRNAs.

Statistics

SPSS version 22.0 (IBM Corporation, Armonk, NY, United States) and GraphPad Prism 8.0 (SanDiego, CA, United States) were utilized to analyze the statistical data. Kaplan-Meier method was employed to construct the overall survival curves. The univariate and multivariate Cox regression analyses were implemented to determine the independent prognostic indicators. The relevance of signature and clinical parameters were assessed by Chi-square test. p < 0.05 was considered as a cutoff value in all circumstances.

Results

Clinical Characteristics of Patients

The schematic diagram of our present study was exhibited in Figure 1. Clinical data derived from TCGA database that contains a total of 448 primary PTC patients whose follow-up time over 30 days was collected for the present study. The detailed information of 448 PTC patients was summarized in Table 1. In Table 1, slightly more than half of the patients (56.25%) was aged over 45 years, and the vast majority of PTC patients were female which accounted for 73.88%. The vital status of most patients (96.65%) was alive. Furthermore, 66.96% of PTC patients were at an early stage of cancer as well as T/N/M stage (T1/T2, 62.72%; M0, 98.44%). Additionally, a total of 56 healthy volunteers were enrolled as the normal group, whose age and gender were matched (Supplementary Table S1).

FIGURE 1

Co-Expression Network Construction and key Module Determination

By using the statistical software R and “edgeR” package, the remarkable DElncRNAs were screened by comparing the PTC tumor cases with normal samples. In total, 1,601 DElncRNAs were identified, including 1,000 upregulated DElncRNAs and 601 down-regulated DElncRNAs. To determine the key modules that most related with clinical traits of PTC, WGCNA on the 1,601 DElncRNAs and 448 samples with complete information was performed. By setting soft-thresholding power and cutting height as 0.25, a total of seven key modules were achieved (Figures 2A–C). In the heatmap of module-trait correlations, we observed that the blue and turquoise modules were highly related with clinical traits (Figure 3A). As the turquoise module possessed the largest number of DElncRNAs and a higher correlation with certain clinical traits (such as T, N, and survival time) relative to the blue module, we selected turquoise module for further analysis. The exact associations between genes of turquoise module and clinical features of PTC patients were revealed in Figure 3B.

FIGURE 2

FIGURE 3

Validation of hub lncRNAs

A LASSO Cox regression model was constructed and shown in Figures 4A,B and 9 hub lncRNAs were eventually identified, including six high-risk lncRNAs SLC12A5-AS1, LINC02028, KIZ-AS1, LINC02019, LINC01877, LINC01444 and three low-risk lncRNAs LINC01176, LINC01290 and LINC00581 (Figure 4C). The risk score was calculated for the nine-lncRNAs signature as follows:

FIGURE 4

Risk score=(0.081*the expression level of “SLC12A5-AS1”)+(0.014*the expression level of LINC02028)+(0.141*the expression level of “KIZ-AS1”)+(−0.013*the expression level of “LINC01176”)+(0.033*the expression level of “LINC02019”)+(0.158*the expression level of “LINC01877”)+(−0.455*the expression level of “LINC01290”)+(0.005*the expression level of “LINC01444”) (−1.039*the expression level of “LINC00581”).

Based on the median value of risk score, 448 PTC patients were classified into two groups: high risk group (n = 224) and low risk group (n = 224). Compared with the low risk group, PTC patients of high risk group showed the poorer outcomes (Figure 5A). Death was taken as the end point for overall survival analysis. The ROC analysis was performed to examine the prognostic value of nine-lncRNAs signature. The area under curve (AUC) of nine-lncRNAs signature was 0.971 (Figure 5B). Results showed that patients in the high risk group had lower Disease-specific survival (DSS) than those in the low risk group (p = 0.0112) (Figure 5C). DSS refers specifically to death due to thyroid cancer as the study end point. Progression-free survival (PFS) of PTC patients was not statistically different between high risk and low risk groups. However, as can be seen from Figure 5D, the PFS of the high risk group was also lower than that of the low risk group. PFS took the appearance of disease progression as the end point. The nomogram indicated that the nine-lncRNAs signature may serve as a prognostic candidate for PTC (Figure 5E).

FIGURE 5

Clinical Correlation of Signature in PTC

To better illustrate the clinical significance of this model, we analyzed the correlation between risk score and clinical features, including age, gender, pathological stage, T/N/M stage, and survival status. The results exhibited that age, pathological stage, T stage, and survival status of PTC patients were correlated with the risk score with statistical significance (Table 2). Furthermore, the number of patients with high risk scores in advanced pathological stage and T stage was more than that of patients with low risk scores. These findings confirmed that this signature was highly related with clinical features of PTC patients.

The Nine-lncRNAs Signature is an Independent Prognostic Factor in PTC

To detect whether the prognostic significance of nine-lncRNAs signature is depended on the clinical parameters, univariate and multivariate Cox regression analyses were adopted to analyze the following variables like nine-lncRNAs signature’s risk score, age, gender, pathological stage, and T/N/M stage. In univariate Cox analysis, risk score (HR = 1.001, p < 0.001), age (HR = 1.114, p < 0.001) and pathological stage (HR = 2.346, p < 0.001) were all powerful variables related with prognosis of PTC patients in TCGA-PTC cohort (Figure 6A). Furthermore, results of multivariate Cox analysis demonstrated that risk score (HR = 1.001, p < 0.001) and age (HR = 1.127, p < 0.001) were still proved to be independent factors based on the TCGA-PTC cohort (Figure 6B). All findings elucidated that the nine-lncRNAs signature might function as an independent prognostic factor in PTC.

FIGURE 6

Function Annotation of the hub lncRNAs and Nine-lncRNAs Signature

In order to predict the possible biological processes (BP) and signaling pathways that the nine lncRNAs may participate in, functional enrichment analyses, GO, and KEGG analyses were performed based on the co-expressed genes of nine hub lncRNAs in PTC. The top 10 GO BP terms with highest enrichment were exhibited in Figure7A, including T cell activation, positive regulation of cell adhesion, focal adhesion, cell adhesion molecule binding, etc. (p < 0.05). In addition, the top 16 pathways that were enriched by the nine lncRNAs in KEGG analysis were presented in Figure7B, including Cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, and Th17 cell differentiation (p < 0.05).

FIGURE 7

In addition to GO and KEGG analyses, GSEA was also implemented to inquire into the potential biological processes and signaling pathways of the nine-lncRNAs signature (Figures 7C,D). Eighteen KEGG pathways were enriched. Among these, three pathways related to PTC pathogenesis, namely adipocytokine signaling pathway, glycolysis and gluconeogenesis, and type_II_diabetes mellitus, were enriched in the high risk group, indicating that the nine-lncRNAs signature may be implicated in the progression of PTC (Figure 7C, p < 0.05). Meanwhile, a total of 782 GO terms were identified and the pathways-related with PTC development were as follows: interleukin 1 mediated signaling pathway, negative regulation of endoplasmic reticulum stress induced intrinsic apoptotic signaling pathway, thyroid hormone generation (Figure 7D, p < 0.05).

Comparisons of Tumor-Infiltrating Immune Cells in High/low Risk Group and Correlation of Nine hub lncRNAs Expression With Tumor-Infiltrating Immune Cells

The tumor microenvironment is composed of tumor cells, stromal cells and infiltrating immune cells. Thus, we tested proportions of 22 types of immune cells in the high/low risk group. Memory B cells, CD8 T cells, follicular helper T cells, monocytes, activated mast cells, and eosinophils were higher in risk than the low-risk group, while M0 macrophages, resting dendritic cells, and resting mast cells were lower in risk than the low-risk group (Figure 8). These findings showed different relationships between signatures and different types of infiltrating immune cells, resulting in different survival outcomes.

FIGURE 8

Then, we used TIMER to determine the associations between expressions of nine hub lncRNAs and tumor infiltrating immune cells. As illustrated in Supplementary Figure S1, except for KIZ-AS1, all other lncRNAs were correlated with the content of more than two kinds of immune cells, which was statistically significant.

Survival Analysis of Nine hub lncRNAs in PTC

As previously described, results of the multivariate Cox regression analyses disclosed that SLC12A5-AS1, LINC02028, KIZ-AS1, LINC02019, LINC01877, and LINC01444 were high-risk genes, while LINC01176, LINC01290, and LINC00581 were low-risk genes (Figure 4C; Supplementary Table S2). To further detect the prognostic significances of nine hub lncRNAs in PTC, we accessed the KM Plotter website. The high expressions of LINC01877, LINC02028, and KIZ-AS1 were linked with the unfavorable outcomes of PTC patients, and low expressions of LINC00581, LINC01290 and LINC01176 were related with poor outcomes, with statistical significance (Figure 9, p < 0.05). There was no data of LINC02019 in the KM Plotter website. In summary, these data provide some support for the predictable role of nine hub lncRNAs in PTC.

FIGURE 9

Diagnostic Nomograms Establishment and Performance Evaluation

To build a potential diagnostic model, we first performed the logistic regression analysis for the identification of important hub lncRNAs. Through screening, six important hub lncRNAs were identified, including LINC02028, KIZ-AS1, LINC01877, LINC01444, LINC01176, and LINC00581. On the basis of these six important lncRNAs, a nomogram was built (Figure 10A). The calibration curve showed high consistence between actual probability and predicted survival proportion (Figure 10B). The AUC of the diagnostic model was 0.975, which suggested that the nomogram might have high sensitivity for PTC diagnosis (Figure 10C). PCA further confirmed the effectiveness of the diagnostic model with six important lncRNAs, and the results showed that the normal and tumor samples can be better distinguished based on the expression level of six important lncRNAs (Figure 10D).

FIGURE 10

Comparison of Prognostic Prediction Score Between our Signature and Published Signatures in PTC or TC

By retrieving the published lncRNA signature in others’ studies, we chose three signatures, including six-lncRNA signature [30], seven-lncRNA signature [31] and nine-lncRNA signature [32], and performed ROC evaluation to assess the prognostic prediction performance based on these lncRNAs. The results exhibited in Figure 11 showed that AUC values of these three signatures were 0.861, 0.776, and 0.936, all lower than AUC value of our model (0.971). The data indicates a favorable potential of our signature for predicting the prognosis of PTC.

FIGURE 11

Discussion

PTC, as the most common malignancy in the endocrine system, is a great challenge for public health [33]. Thus, it is of great urgency to identify specific cancer-related biomarkers for risk evaluation, which could be employed to predict the prognosis of PTC patients and promote the development of effective therapies to PTC. In the investigation, a systematic bioinformatics analysis identified a nine-lncRNAs signature that significantly correlated with the pathogenesis of PTC.

WGCNA is a systematic biological algorithm that could distinguish highly synergistically altered genes and unearth relevance between genes and clinical characteristics [34]. Therefore, it has been widely utilized to identify the candidate biomarkers of diverse diseases, such as Alzheimer’s disease, NSCLC, colon cancer and PTC [35–38]. For example, a total of 11 genes have been reported to exert important roles in PTC relapse, which shed novel insights in the lymph node metastasis of PTC [39]. Han et al. determined that nine genes containing CDH5, KDR, CD34, FLT4, EMCN, FLT1, ROBO4, PTPRB, and CD93 could mediate the I¹³¹ radiotherapy in TC patients [40]. Three genes (NGF, FOS, and GRIA1) are found to be involved in regulation of the features of TC stem cells using WGCNA [41]. Additionally, for a long time in the past, genomic libraries of many non-protein transcripts, including lncRNAs, were considered as insignificant transcription “junk”. However, the implementation of TCGA program has focused a tremendous amount of attention on the roles of lncRNAs in tumorigenesis [42, 43]. Therefore, in the present study, to investigate possible prognostic lncRNAs for PTC, we downloaded the lncRNAs expression profiles from TCGA database and mined the underlying microarrays to perform a comprehensive analysis of lncRNA expression matrix files. In the TCGA-PTC cohort, a total of 1,601 DElncRNAs, 1,000 up-regulated lncRNAs and 601 down-regulated lncRNAs, were identified. Co-expression network construction and key modules identification through WGCNA identified seven modules, of which, the genes of turquoise module were all linked with age, survival status, survival time, stage, stage of tumor, and stage of metastasis.

After analyzing the genes of turquoise module by LASSO analysis, we finally obtained nine hub lncRNAs, including SLC12A5-AS1, LINC02028, KIZ-AS1, LINC02019, LINC01877, LINC01444, LINC01176, LINC01290, and LINC00581. The overall survival curve was plotted by Kaplan-Meier and revealed that PTC patients in the high risk group had unfavorable outcomes compared with the patients in the low risk group. The ROC analysis elucidated that the AUC of the identified nine-lncRNAs signature for overall survival rates was 0.971, indicating a satisfactory predictive value for PTC prognosis. The nomogram for 1-, 2- and 3-years overall survival suggested that the nine-lncRNAs signature was suitable for predicting the probability of overall survival of PTC patients. Furthermore, univariate and multivariate Cox regression analyses were adopted to identify the associations between following factors like risk score, age, gender, stage and T/N/M stage and PTC prognosis, which demonstrated that the nine-lncRNAs signature may be an independent candidate for predicting the prognosis of PTC patients.

Notably, up to date, there are fe reports regarding the biological roles of nine hub lncRNAs. Therefore, the results of the functional enrichment analyses could further promote the understanding of potential functions of nine hub lncRNAs. In GO and KEGG analyses, 2036 mRNAs that co-expressed with the nine hub lncRNAs were enriched in 228 GO terms and 16 pathways, including T cell activation, positive regulation of cell adhesion, focal adhesion, cell adhesion molecule binding, Cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, and Th17 cell differentiation. As previously described, the activation of T cell intracellular antigen can suppress the development of anaplastic TC [44]. Song et al. indicated that the variants of genes-related with cell adhesion could be utilized to predict PTC aggressiveness [45]. NF-kappa B signaling pathway has been identified as a putative target for the treatment of advanced TC [46]. Besides, the underlying mechanism of nine-lncRNAs signature was detected using GSEA and 18 pathways were obtained, such as adipocytokine signaling pathway, glycolysis and gluconeogenesis, type II diabetes mellitus as well as 782 GO terms (for instance, interleukin 1 mediated signaling pathway, negative regulation of endoplasmic reticulum stress induced intrinsic apoptotic signaling pathway, thyroid hormone generation). It is well known that adipokines could regulate body mass and affect glucose homeostasis, cell proliferation, and other important cell procedures, which are all involved in tumor growth. Moreover, the deregulation of adipokines is observed in the head and neck cancer [47]. The energy provision of tumor cells relies nearly entirely on glycolysis and gluconeogenesis [48]. A previous report suggests that early type II diabetes mellitus might have a protective effect against TC [49]. Interleukin 1 is involved in the progression of medullary TC via mediating other essential pathways [50]. These findings of functional enrichment analyses determined that the nine hub lncRNAs may be required for PTC by modulating various signaling pathways and biological processes.

Tumor microenvironment contains tumor cells and tumor infiltrating immune cells. In order to explore the relationship of this signature and immune infiltration, we assessed the differences of 22 types of tumor infiltrating immune cells in the high/low risk score group. The results displayed that memory B cells, CD8 T cells, follicular helper T cells, monocytes, activated mast cells, eosinophils were higher in high-risk than low-risk group, while M0 macrophages, resting dendritic cells, and resting mast cells were lower in risk than the low-risk group. In the tumor environment, mast cells stimulate angiogenesis resulting in the growth and development of tumors [51]. Monocytes can indirectly promote extracellular matrix degradation and tumor angiogenesis, thereby promoting tumor progression and metastasis [52]. Studies have found that low CD4+/CD8 + ratio is related to a poor survival of tumor patients [53]. To further investigate the possible functions of this signature, nine hub lncRNAs, we referred to the TIMER website and found that except for KIZ-AS1, all other lncRNAs were correlated with the content of more than two kinds of immune cells with significant differences. These suggested that the functions of this signature composed of nine hub lncRNAs related to immunological regulation of the tumor microenvironment.

Additionally, the diagnostic model was then built using logistic regression, containing LINC02028, KIZ-AS1, LINC01877, LINC01444, LINC01176, and LINC00581. Calibration plot and ROC curves unearthed the high accuracy and sensitivity of this diagnostic model, which was enhanced by PCA analysis.

Nevertheless, a few limitations of this present study needed to be noted. Firstly, because all of the clinical samples were derived from TCGA database, we might overlook other potential lncRNAs candidates. A larger sample size should be used to validate the results of this study. Secondly, the data of LINC02019 in KM Plotter website were missing, therefore PTC samples should be harvested in the future for further validation of nine hub genes’ potential role in PTC. Thirdly, the in vitro and vivo functional experiments are required to verify their accurate biological effects. Fourth, the tumor microenvironment of PTC patients has been proven to be affected by concomitant Hashimoto Thyroiditis(HT). The coexistence of HT is related with the inf immunophenotype, whereas PTC patients without HT are mainly characterized by immune desert (ID) or immune-excluded (IE) immunophenotype [54]. Thus, further research should be performed to explore the relationship between PTC and HT as well as the accurate phenotype of immune cells.

Conclusion

In conclusion, a nine-lncRNAs signature was constructed in our present study, which might be an independent prognostic variable for patients with PTC and possess potential value for prognosis and risk assessment. A predictive model with high accuracy and sensitivity was also built for PTC diagnosis. These findings would advance our understanding of PTC pathogenesis and optimize the strategies for improving prognosis of PTC.

References

• 1.

  Bray F Ferlay J Soerjomataram I Siegel RL Torre LA Jemal A . Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: a Cancer J clinicians (2018) 68:394–424. 10.3322/caac.21492

  • CrossRef
  • Google Scholar

• 2.

  Sung H Ferlay J Siegel RL Laversanne M Soerjomataram I Jemal A et al Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A Cancer J Clin (2021) 71:209–49. 10.3322/caac.21660

  • CrossRef
  • Google Scholar

• 3.

  Fagin JA Wells SA . Biologic and Clinical Perspectives on Thyroid Cancer. N Engl J Med (2016) 375:1054–67. 10.1056/nejmra1501993

  • CrossRef
  • Google Scholar

• 4.

  Singer PA Cooper DS Daniels GH Ladenson PW Greenspan FS Levy EG et al Treatment Guidelines for Patients with Thyroid Nodules and Well-Differentiated Thyroid Cancer. American Thyroid Association. Am Thyroid Assoc Arch Intern Med (1996) 156:2165–72. 10.1001/archinte.156.19.2165

  • CrossRef
  • Google Scholar

• 5.

  Mirian C Grønhøj C Jensen DH Jakobsen KK Karnov K Jensen JS et al Trends in Thyroid Cancer: Retrospective Analysis of Incidence and Survival in Denmark 1980-2014. Cancer Epidemiol (2018) 55:81–7. 10.1016/j.canep.2018.05.009

  • CrossRef
  • Google Scholar

• 6.

  Clarke CA Reynolds P Oakley-Girvan I Lee E Lu Y Yang J et al Indicators of Microbial-Rich Environments and the Development of Papillary Thyroid Cancer in the California Teachers Study. Cancer Epidemiol (2015) 39:548–53. 10.1016/j.canep.2015.04.014

  • CrossRef
  • Google Scholar

• 7.

  Toniato A Boschin I Casara D Mazzarotto R Rubello D Pelizzo M . Papillary Thyroid Carcinoma: Factors Influencing Recurrence and Survival. Ann Surg Oncol (2008) 15:1518–22. 10.1245/s10434-008-9859-4

  • CrossRef
  • Google Scholar

• 8.

  Garzon R Calin GA Croce CM . MicroRNAs in Cancer. Annu Rev Med (2009) 60:167–79. 10.1146/annurev.med.59.053006.104707

  • CrossRef
  • Google Scholar

• 9.

  Jaé N Heumüller AW Fouani Y Dimmeler S . Long Non-coding RNAs in Vascular Biology and Disease. Vasc Pharmacol (2019) 114:13–22. 10.1016/j.vph.2018.03.003

  • CrossRef
  • Google Scholar

• 10.

  Rinn JL Chang HY . Genome Regulation by Long Noncoding RNAs. Annu Rev Biochem (2012) 81:145–66. 10.1146/annurev-biochem-051410-092902

  • CrossRef
  • Google Scholar

• 11.

  Mercer TR Dinger ME Mattick JS . Long Non-coding RNAs: Insights into Functions. Nat Rev Genet (2009) 10:155–9. 10.1038/nrg2521

  • CrossRef
  • Google Scholar

• 12.

  Wilusz JE Sunwoo H Spector DL . Long Noncoding RNAs: Functional Surprises from the RNA World. Genes Dev (2009) 23:1494–504. 10.1101/gad.1800909

  • CrossRef
  • Google Scholar

• 13.

  Spizzo R Almeida MI Colombatti A Calin GA . Long Non-coding RNAs and Cancer: a New Frontier of Translational Research. Oncogene (2012) 31:4577–87. 10.1038/onc.2011.621

  • CrossRef
  • Google Scholar

• 14.

  Fatica A Bozzoni I . Long Non-coding RNAs: New Players in Cell Differentiation and Development. Nat Rev Genet (2014) 15:7–21. 10.1038/nrg3606

  • CrossRef
  • Google Scholar

• 15.

  Kornienko AE Guenzl PM Barlow DP Pauler FM . Gene Regulation by the Act of Long Non-coding RNA Transcription. BMC Biol (2013) 11:59. 10.1186/1741-7007-11-59

  • CrossRef
  • Google Scholar

• 16.

  Lei H Gao Y Xu X . LncRNA TUG1 Influences Papillary Thyroid Cancer Cell Proliferation, Migration and EMT Formation through Targeting miR-145. Acta Biochim Biophys Sinica (2017) 49:588–97. 10.1093/abbs/gmx047

  • CrossRef
  • Google Scholar

• 17.

  Gu Y Feng C Liu T Zhang B Yang L The Downregulation of lncRNA EMX2OS Might Independently Predict Shorter Recurrence-free Survival of Classical Papillary Thyroid Cancer. PLoS ONE (2018). 13:e0209338. 10.1371/journal.pone.0209338

  • CrossRef
  • Google Scholar

• 18.

  Wang Z Cheng Y Zhao D Pliss A Liu J Luan P . Synergic Treatment of Alzheimer's Disease with Brain Targeted Nanoparticles Incorporating NgR-siRNA and Brain Derived Neurotrophic Factor. Smart Mater Med (2020) 1:125–30. 10.1016/j.smaim.2020.08.001

  • CrossRef
  • Google Scholar

• 19.

  Chen F Yin S Zhu J Liu P Yang C Feng Z et al lncRNA DGCR5 Acts as a Tumor Suppressor in Papillary Thyroid Carcinoma via Sequestering miR-2861. Exp Ther Med (2019) 17:895–900. 10.3892/etm.2018.7012

  • CrossRef
  • Google Scholar

• 20.

  Liu K Kang M Liao X Wang R . Genome-wide Investigation of the Clinical Significance and Prospective Molecular Mechanism of Minichromosome Maintenance Protein Family Genes in Patients with Lung Adenocarcinoma. PLoS ONE (2019) 14:e0219467. 10.1371/journal.pone.0219467

  • CrossRef
  • Google Scholar

• 21.

  Subbannayya Y Di Fiore R . The Role of Omics Approaches to Characterize Molecular Mechanisms of Rare Ovarian Cancers. Recent Adv Future Perspect (2021) 9. 10.3390/biomedicines9101481

  • CrossRef
  • Google Scholar

• 22.

  Tang J Cui Q Zhang D Kong D Liao X Ren J et al A Prognostic eight‐lncRNA Expression Signature in Predicting Recurrence of ER‐positive Breast Cancer Receiving Endocrine Therapy. J Cel Physiol (2020) 235:4746–55. 10.1002/jcp.29352

  • CrossRef
  • Google Scholar

• 23.

  Tang J Ren J Cui Q Zhang D Kong D Liao X et al A Prognostic 10-lncRNA Expression Signature for Predicting the Risk of Tumour Recurrence in Breast Cancer Patient. J Cel Mol Med (2019) 23:6775–84. 10.1111/jcmm.14556

  • CrossRef
  • Google Scholar

• 24.

  Wang S Chai K Chen J . A Novel Prognostic Nomogram Based on 5 Long Non-coding RNAs in clear Cell Renal Cell Carcinoma. Oncol Lett (2019) 18:6605–13. 10.3892/ol.2019.11009

  • CrossRef
  • Google Scholar

• 25.

  Wu M Yuan H Li X Liao Q Liu Z . Identification of a Five-Gene Signature and Establishment of a Prognostic Nomogram to Predict Progression-free Interval of Papillary Thyroid Carcinoma. Front Endocrinol (2019) 10:790. 10.3389/fendo.2019.00790

  • CrossRef
  • Google Scholar

• 26.

  Yang G Zhang Y Yang J . A Five-microRNA Signature as Prognostic Biomarker in Colorectal Cancer by Bioinformatics Analysis. Front Oncol (2019) 9:1207. 10.3389/fonc.2019.01207

  • CrossRef
  • Google Scholar

• 27.

  Niemira M Collin F Szalkowska A Bielska A Chwialkowska K Reszec J et al Molecular Signature of Subtypes of Non-small-cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-expression Network Analysis (WGCNA). Cancers (Basel) (2019) 12:12. 10.3390/cancers12010037

  • CrossRef
  • Google Scholar

• 28.

  Wu S Zheng J Li Y Wu Z Shi S Huang M et al Development and Validation of an MRI-Based Radiomics Signature for the Preoperative Prediction of Lymph Node Metastasis in Bladder Cancer. EBioMedicine (2018) 34:76–84. 10.1016/j.ebiom.2018.07.029

  • CrossRef
  • Google Scholar

• 29.

  Weinstein JN Collisson EA Collisson EA Mills GB Shaw KRM Ozenberger BA et al The Cancer Genome Atlas Pan-Cancer Analysis Project. Nat Genet (2013) 45:1113–20. 10.1038/ng.2764

  • CrossRef
  • Google Scholar

• 30.

  Li Z Lin W Zheng J Hong W Zou J Zhang T et al Identification of Immune-Related lncRNAs to Improve the Prognosis Prediction for Patients with Papillary Thyroid Cancer. Biosci Rep (2021) 41. 10.1042/bsr20204086

  • CrossRef
  • Google Scholar

• 31.

  Huang Z Xiao C Zhang F Zhou Z Yu L Ye C et al A Novel Framework to Predict Breast Cancer Prognosis Using Immune-Associated LncRNAs. Front Genet (2020) 11:634195. 10.3389/fgene.2020.634195

  • CrossRef
  • Google Scholar

• 32.

  Shan Y He R Yang X Zang S Yao S Gao M et al An Autophagy-Related lncRNA Prognostic Risk Model for Thyroid Cancer. Eur Arch Otorhinolaryngol (2021):1–14. 10.1007/s00405-021-07134-4

  • CrossRef
  • Google Scholar

• 33.

  Schlumberger M Leboulleux S . Current Practice in Patients with Differentiated Thyroid Cancer. Nat Rev Endocrinol (2021) 17:176–88. 10.1038/s41574-020-00448-z

  • CrossRef
  • Google Scholar

• 34.

  Luo Y Shen D Chen L Wang G Liu X Qian K et al Identification of 9 Key Genes and Small Molecule Drugs in clear Cell Renal Cell Carcinoma. Aging (2019) 11:6029–52. 10.18632/aging.102161

  • CrossRef
  • Google Scholar

• 35.

  Sun J Zhang Z Bao S Yan C Hou P Wu N et al Identification of Tumor Immune Infiltration-Associated lncRNAs for Improving Prognosis and Immunotherapy Response of Patients with Non-small Cell Lung Cancer. J Immunother Cancer (2020) 8:e000110. 10.1136/jitc-2019-000110

  • CrossRef
  • Google Scholar

• 36.

  Zhai X Xue Q Liu Q Guo Y Chen Z . Colon Cancer Recurrence-Associated Genes Revealed by WGCNA Co-expression Network Analysis. Mol Med Rep (2017) 16:6499–505. 10.3892/mmr.2017.7412

  • CrossRef
  • Google Scholar

• 37.

  Wu Y Liang S Zhu H Zhu Y . Analysis of Immune-Related Key Genes in Alzheimer's Disease. Bioengineered (2021). 10.1080/21655979.2021.1999553

  • CrossRef
  • Google Scholar

• 38.

  Lin R Fogarty CE Ma B Li H Ni G Liu X et al Identification of Ferroptosis Genes in Immune Infiltration and Prognosis in Thyroid Papillary Carcinoma Using Network Analysis. BMC genomics (2021) 22:576. 10.1186/s12864-021-07895-6

  • CrossRef
  • Google Scholar

• 39.

  Zhai T Muhanhali D Jia X Wu Z Cai Z Ling Y . Identification of Gene Co-expression Modules and Hub Genes Associated with Lymph Node Metastasis of Papillary Thyroid Cancer. Endocrine (2019) 66:573–84. 10.1007/s12020-019-02021-9

  • CrossRef
  • Google Scholar

• 40.

  Shuwen H Xi Y Miao D Jiamin X Jing Z Weili G . Nine Genes Mediate the Therapeutic Effects of Iodine-131 Radiotherapy in Thyroid Carcinoma Patients. Dis Markers (2020) 2020:9369341. 10.1155/2020/9369341

  • CrossRef
  • Google Scholar

• 41.

  Zhang T Yan C Ye Z Yin X Jiang TA . The Identification of Three Key Genes Related to Stemness in Thyroid Carcinoma through Comprehensive Analysis. Comb Chem High Throughput Screen (2021) 24:423–32. 10.2174/1386207323666200806164003

  • CrossRef
  • Google Scholar

• 42.

  Li J Han L Roebuck P Diao L Liu L Yuan Y et al TANRIC: An Interactive Open Platform to Explore the Function of lncRNAs in Cancer. Cancer Res (2015) 75:3728–37. 10.1158/0008-5472.can-15-0273

  • CrossRef
  • Google Scholar

• 43.

  Song P Jiang B Liu Z Ding J Liu S Guan W . A Three-lncRNA Expression Signature Associated with the Prognosis of Gastric Cancer Patients. Cancer Med (2017) 6:1154–64. 10.1002/cam4.1047

  • CrossRef
  • Google Scholar

• 44.

  Bi JW Zou YL Qian JT Chen WB . MiR-599 Serves a Suppressive Role in Anaplastic Thyroid Cancer by Activating the T-Cell Intracellular Antigen. Exp Ther Med (2019) 18:2413–20. 10.3892/etm.2019.7864

  • CrossRef
  • Google Scholar

• 45.

  Song J Wu S Xia X Wang Y Fan Y Yang Z . Cell Adhesion-Related Gene Somatic Mutations Are Enriched in Aggressive Papillary Thyroid Microcarcinomas. J Transl Med (2018) 16:269. 10.1186/s12967-018-1642-0

  • CrossRef
  • Google Scholar

• 46.

  Pozdeyev N Berlinberg A Zhou Q Wuensch K Shibata H Wood WM et al Targeting the NF-Κb Pathway as a Combination Therapy for Advanced Thyroid Cancer. PloS one (2015) 10:e0134901. 10.1371/journal.pone.0134901

  • CrossRef
  • Google Scholar

• 47.

  Tzanavari T Tasoulas J Vakaki C Mihailidou C Tsourouflis G Theocharis S . The Role of Adipokines in the Establishment and Progression of Head and Neck Neoplasms. Cmc (2019) 26:4726–48. 10.2174/0929867325666180713154505

  • CrossRef
  • Google Scholar

• 48.

  Halestrap AP Wilson MC . The Monocarboxylate Transporter Family-Role and Regulation. IUBMB life (2012) 64:109–19. 10.1002/iub.572

  • CrossRef
  • Google Scholar

• 49.

  Seo YG Choi HC . The Association between Type 2 Diabetes Mellitus and Thyroid Cancer (2017) 2017:5850879. 10.1155/2017/5850879

  • CrossRef
  • Google Scholar

• 50.

  Park J-I Strock CJ Ball DW Nelkin BD . Interleukin-1β Can Mediate Growth Arrest and Differentiation via the Leukemia Inhibitory factor/JAK/STAT Pathway in Medullary Thyroid Carcinoma Cells. Cytokine (2005) 29:125–34. 10.1016/j.cyto.2004.10.005

  • CrossRef
  • Google Scholar

• 51.

  Melillo RM Guarino V Avilla E Galdiero MR Liotti F Prevete N et al Mast Cells Have a Protumorigenic Role in Human Thyroid Cancer. Oncogene (2010) 29:6203–15. 10.1038/onc.2010.348

  • CrossRef
  • Google Scholar

• 52.

  Matsubara T Kanto T Kuroda S Yoshio S Higashitani K Kakita N et al TIE2-expressing Monocytes as a Diagnostic Marker for Hepatocellular Carcinoma Correlates with Angiogenesis. Hepatology (2013) 57:1416–25. 10.1002/hep.25965

  • CrossRef
  • Google Scholar

• 53.

  Brambilla E Le Teuff G Marguet S Lantuejoul S Dunant A Graziano S et al Prognostic Effect of Tumor Lymphocytic Infiltration in Resectable Non-small-cell Lung Cancer. Jco (2016) 34:1223–30. 10.1200/jco.2015.63.0970

  • CrossRef
  • Google Scholar

• 54.

  Sulaieva O Chernenko O Selesnov O Nechay O Maievskyi O Falalyeyeva T et al Mechanisms of the Impact of Hashimoto Thyroiditis on Papillary Thyroid Carcinoma Progression: Relationship with the Tumor Immune Microenvironment. Endocrinol Metab (2020) 35:443–55. 10.3803/enm.2020.35.2.443

  • CrossRef
  • Google Scholar