A novel autophagy-related long non-coding RNAs signature predicting progression-free interval and I-131 therapy benefits in papillary thyroid carcinoma

2.1 Data collection and processing

The workflow of this study is presented in Figure 1. The expression data including mRNA and lncRNA of PTC patients were got from TCGA database (https://www.cbioportal.org) [24]. This dataset consisted of 507 subjects, 510 tumor samples, and 58 adjacent normal tissue samples. Only the patients with clinical data and follow-up time/PFI ≥ 30 days were included. Overall, 498 cases were finally included and were randomly split into a training cohort (n = 299) and a validation cohort (n = 199). A total 222 autophagy-related genes (ARGs) were obtained from the Human Autophagy Database (HADb, http://autophagy.lu/clustering/index.html), which collected those genes from the literature. LncRNAs and ARGs mRNAs expression were extracted according to GENCODE annotations (https://www.gencodegenes.org). LncRNAs with zero expression levels in more than 50% of samples were excluded. Clinical variables including age, gender, race, cancer history, thyroid gland disorder history, histological types, TNM stages, T stage, N stage, M stage, tumor location, residual tumor, American Thyroid Association (ATA) risk stratification, the distant metastasis, patient age, completeness of resection, local invasion, and tumor size (MACIS) scores were collected. The methods of assessing the tumors with ATA risk stratification and MACIS scores have been described in the previous report [24] and are introduced in the Supplementary data. We also re-evaluated the tumors based on the methods. Additionally, BRAF mutation and telomerase reverse transcriptase (TERT) promoter mutation status information were extracted, which have been reported to be associated with prognosis in TC [25]. Patients’ progression-free survival (PFS) and overall survival (OS) were also extracted.

Figure 1

Flowchart.

2.2 Screening of the differentially expressed ARGs-related lncRNAs

We used the “Limma” package to identify the differentially expressed lncRNAs (DElncs) between tumor samples and control samples with the thresholds of fold change (FC) >1.5 or <1/1.5 and false discovery rate (FDR) <0.05 [26]. Subsequently, the correlation coefficient (R ²) between ARGs and DElncs was calculated by Pearson correlation analyses. The lncRNAs with R ² > 0.25 and P < 0.001 were defined as autophagy-related lncRNAs.

2.3 Construction of autophagy-related lncRNAs signature

Univariate Cox regression analyses were performed to explore the associations of the autophagy-related DElncs with PFI. The autophagy-related DElncs that were associated with PFI in both the entire cohort and the training cohort (P < 0.1) were used as candidates to build the prognostic signature in the training cohort. The candidate DElncs were imputed into LASSO Cox regression analysis. The core autophagy-related DElncs tightly related to PFI were obtained when the optimal lambda value was achieved. The selected lncRNAs were then subjected to stepwise multivariate Cox regression analysis to build the autophagy-related lncRNAs risk signature. Each patient’s risk score was calculated by a linear combination of multiplying each lncRNA expression and the corresponding Cox regression coefficient. The patients were stratified into high-risk and low-risk groups by the median value. The PFI difference between the two groups was compared by Kaplan–Meier curves. Multivariate Cox regression analysis was performed to evaluate the independent prognostic value of the risk scores by adjusting the potential confounders. For the adjusted model I, the confounders was selected if they changed the effect estimate of the risk scores on PFI by more than 10% or were significantly associated with PFI. For the adjusted model II, the confounders in the adjusted model I and the remaining demographic data were adjusted. In addition, the area under the time-dependent ROC curves (AUC), as well as Harrell’s concordance index (C-index), were utilized to assess the predictive value of the autophagy-related lncRNAs risk signature. Validation was performed in the validation cohort and the entire cohort. The predictive performance of the risk model was also compared with TNM stages, ATA risk stratification, and MACIS scores by comparing their C-indices in the patients without missing values.

2.4 Association of the risk score with I-131 therapy efficacy

The association of the risk score with I-131 therapy efficacy was also explored with PFI as the primary endpoint and PFS and OS as the secondary endpoints. The patients were first divided into high and low risk with the optimal cutoff value, which was obtained based on the minimum P-value of the interaction test in univariate Cox analysis with PFI as the primary endpoint. I-131 therapy efficacy was investigated in patients with high or low risk by Kaplan–Meier curves with log-rank tests.

2.5 Gene set enrichment analysis (GSEA)

GSEA (version 3.0) was executed to assess the significantly different hallmark gene sets between different risk subgroups. The enriched gene sets whose/normalized enrichment score (NES)/> 1, nominal P-value <0.05, and FDR q-value <0.05 were treated as significant.

2.6 Single-cell RNA sequencing data analysis

The scRNA-seq data of 22 fresh surgical samples from six primary PTC tumors, six paired adjacent normal tissues, eight metastatic lymph nodes (LNs, including three recurrent LNs), and two subcutaneous metastatic loci were extracted from the GSE184362 dataset stored in the Gene Expression Omnibus (GEO) database [27]. Further analyses were performed with the R package, “Seurat,” using the standard data analysis pipeline [28]. Briefly, cells with low quality (the proportion of mitochondrial genes counts >10%, UMIs <500 or >5,000) were first removed; then the cell gene expression matrix was normalized and scaled with the default parameters; subsequently, the top 2,000 highly variable genes were identified by the FindVariableFeatures() function for the principal component analysis; fourth, the functions, FindNeighbors() and FindClusters() were utilized for cell clustering at a resolution of 0.4; next, the uniform manifold approximation and projection (UMAP) were used for visualization; cell markers provided in the previous report were used to annotate the cell clusters [27]. The cells with positive expression of lncRNAs of interest were assigned to the positive group and the other cells to the negative group. The differentially expressed genes between the two groups were identified through the FindAllMarkers() function, ranked by the log FC, and subjected to GSEA to explore the potential mechanisms.

2.7 Statistical analysis

R software v3.4.3 was utilized to perform all statistical analyses. The categorical variables were compared by chi-square tests or Fisher’s exact tests while the continuous variables were compared by Wilcox tests or Kruskal–Wallis tests. Kaplan–Meier curves were plotted to assess the prognosis differences between different groups stratified by risk scores and I-131 therapy efficacy. Cox regression analyses were carried out to identify lncRNAs that were associated with PFI. P < 0.05 was regarded as statistically significant unless otherwise stated.

1. Ethics statement: The current study has been approved by the Institutional Review Board of Tianjin Union Medical Center of Nankai University (2021-B34).

3.1 Identification of differentially expressed autophagy-related lncRNAs

The mRNA expression data of PTC tissues and adjacent normal tissues were obtained from the TCGA database. A total of 222 ARGs were generated from the HADb database and their expression data were extracted. About 14,822 lncRNAs were also extracted and 7,425 lncRNAs whose expressions were non-zero in more than 50% of samples were finally included. DElncs were identified between PTC and normal thyroid tissues via the “Limma” package. Based on the cutoff criteria, 262 DElncs with 178 lncRNAs downregulated and 84 lncRNAs upregulated were found. The volcano plot and heatmap for these 262 DElncs in tumor and adjacent tissues are displayed in Figure 2. Pearson correlation analyses were performed between these lncRNAs and 222 ARGs. Finally, 199 autophagy-related DElncs were attained based on the criteria of R ² > 0.25 and P < 0.001.

Figure 2

Comparison of lncRNAs expression in PTC tissues with adjacent tissues. (a) Volcano plot. The top ten DElncs were indicated based on the FDR, which were all downregulated. (b) Heatmap plot. Note: lfc, log2(FC).

3.2 Establishment of autophagy-related lncRNA risk signature

Four hundred ninety-eight PTC patients were included and randomly divided into a training cohort (n = 299) and a validation cohort (n = 199). The clinical characteristics of the two cohorts were similar (Table 1). Univariate Cox regression analyses were performed to explore the associations of the autophagy-related DElncs with PFI in the entire cohort and training cohort. The autophagy-related DElncs that were associated with PFI in both the entire and training cohort (P < 0.1) were imputed into LASSO Cox regression analysis in the training cohort. Thirteen core lncRNAs tightly related to PFI were obtained at optimal lambda value (Figure 3). The selected lncRNAs were then subjected to stepwise multivariate Cox regression analysis. Finally, an autophagy-related lncRNAs risk signature consisting of six lncRNAs was constructed by multiplying each lncRNA expression and the corresponding Cox regression coefficient (Table 2). The six lncRNAs were co-expressed with 31 ARGs (R ² > 0.25; Figure 4).

Figure 3

Screening of the core autophagy-related lncRNAs associated with PFI by LASSO COX regression analysis. (a) LASSO coefficient of the lncRNAs by 10-fold cross-validation. (b) Partial likelihood deviance with corresponding log(λ) values at the minimal deviance.

Figure 4

Correlation of the selected six lncRNAs with ARGs.

3.3 Validation of the autophagy-related lncRNA risk signature

The prognostic performance of the autophagy-related lncRNAs risk model was further assessed by Kaplan–Meier plotting curves and time-dependent ROC curves in the training, validation, and entire cohorts. The patients in these cohorts were split into high-risk and low-risk groups with the median value of the risk scores. In the training cohort, the risk curve and the scatterplot showed that the low-risk group had lower risk scores and progression rates (Figure 5a). The corresponding expression profiles of the six autophagy-related lncRNAs were also visualized by heatmap (Figure 5a). The AUCs of 1-, 3-, and 5-year PFI were 0.78, 0.79, and 0.76, respectively (Figure 5b). Kaplan–Meier curve indicated that patients with low risk exerted more favorable PFI than those with high risk (P < 0.0001, Figure 5c). In the validation cohort, the risk score distribution, progression status, and the corresponding expression profiles of the six lncRNAs were also determined (Figure 6a). The AUCs of 1-, 3-, and 5-year PFI were 0.64, 0.71, and 0.79, respectively (Figure 6b). Kaplan–Meier curve showed that patients with low risk exerted more favorable PFI than those with high risk (P = 0.0033, Figure 6c). Similar results were identified in the entire cohort (Figure 7). The Harrell’s C-indices (95% CIs) of the six lncRNAs risk signature were 0.776 (0.692, 0.861), 0.717 (0.593, 0.841), and 0.756 (0.686, 0.825) in the training, validation, and entire cohorts, respectively. The predictive performance of the lncRNAs risk model was also compared with TNM stages, ATA risk stratification, and MACIS scores. In the patients without missing values (n = 480), Harrell’s C-index of the six lncRNAs signature was 0.759 and bigger than that of TNM stages, ATA risk stratification, and MACIS scores, whose C-indices were 0.631, 0.651, and 0.646, respectively (all P < 0.05; Table S1). Collectively, these results revealed that the risk signature exhibited good performance to predict the PFI of the PTC patients, and was superior to TNM stages, ATA risk stratification, and MACIS scores.

Figure 5

Predictive performance of the six autophagy-related lncRNAs signature in the training cohort. (a) Risk scores, PFI/progression status, and expression heatmap. (b) Time-dependent ROC curves for predicting 1-, 3-, and 5-year PFI based on the risk scores. (c) Kaplan–Meier curves of the patients with high and low-risk scores, which were divided by the median value.

Figure 6

Predictive performance of the six autophagy-related lncRNAs signature in the validation cohort. (a) Risk scores, PFI/progression status, and expression heatmap. (b) Time-dependent ROC curves for predicting 1-, 3-, and 5-year PFI based on the risk scores. (c) Kaplan–Meier curves of the patients with high and low-risk scores, which were divided by the median value of the risk scores.

Figure 7

Predictive performance of the six autophagy-related lncRNAs signature in the entire cohort. (a) Risk scores, PFI/progression status, and expression heatmap. (b) Time-dependent ROC curves for predicting 1-, 3-, and 5-year PFI based on the risk scores. (c) Kaplan–Meier curves of the patients with high and low-risk scores, which were divided by the median value of the risk scores.

Finally, multivariate analysis results suggested that, whether in the training, validation, or entire cohorts, the six-lncRNAs risk scores (serve as a continuous variable) were an independent prognostic factor (Table 3). When the risk scores were equally split into two groups or three groups, high-risk scores were also identified as an independent prognostic factor.

3.4 Associations of the risk scores with clinical features and BRAF/TERT mutation status in PTC patients

The associations of the six-lncRNAs risk scores with patients’ clinical features and BRAF/TERT mutation status were analyzed in the entire cohort. Patients were divided into high and low-risk groups by the median value and the features in different risk groups were compared. The results suggested that the patients with high risk had a higher prevalence of tall cell carcinomas, advanced stages, high ATA risk stratification, and higher TERT promoter mutation rate (Table 4). Comparisons of the original risk scores in patients with different features suggested patients with tall cell carcinomas (Figure 8a), III–IV stages (Figure 8b), high ATA risk stratification (Figure 8c), mutated BRAF (Figure 8d), or mutated TERT promoter (Figure 8e) had higher risk scores. The lncRNAs signature risk scores were also found to be positively correlated with MACIS scores (Figure 8f).

Figure 8

Comparison of the risk scores in patients with different features. (a) Comparison of the risk scores in patients with different histological types. (b) Comparison of the risk scores in patients with different TNM stages. (c) Comparison of the risk scores in patients with different ATA risk stratification. (d) Comparison of the risk scores in patients with or without BRAF mutation. (e) Comparison of the risk scores in patients with or without TERT mutation (f) Correlation of the risk scores and MACIS scores. Note: ns, not significant; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

3.5 Association of the risk scores with I-131 therapy benefits

Of the 390 patients, 215 (60%) with treatment information received postoperative I-131 therapy. No significant association of I-131 therapy with patients’ PFI was identified (P = 0.51; Figure 9a). However, when the patients were equally split into three groups by the LncRNAs signature risk scores, we observed a trend that I-131 therapy was associated with favorable PFI in the high-risk group while being associated with poor PFI in the low-risk group (data not shown), suggesting that the patients with high-risk scores but not those with low-risk scores might get benefit from I-131 therapy. To identify the patients favorable to I-131 therapy, the patients were divided into two groups (unfavorable group and favorable group) by a series of cutoff values, and the optimal cutoff value was selected by minimal P-value for interaction tests in univariate Cox analyses. One hundred ten (28%) and 280 (72%) patients in favorable (with high-risk scores) and unfavorable (with low-risk scores) groups were identified, respectively. As expected, I-131 therapy was associated with poor PFI in patients from the unfavorable group (P = 0.23; Figure 9b) while being associated with improved PFI in patients from the favorable group (P = 0.057; Figure 9c) although the associations did not reach statistical significance. We also used PFS and OS as secondary endpoints to explore the effects of the risk groups on I-131 therapy benefits. Although improved PFS was observed after I-131 therapy in whole 390 patients, the PFS difference was not significant (P = 0.34; Figure 9d). Similarly, a trend of an association between I-131 therapy and poor PFS was found in the unfavorable group (P = 0.15; Figure 9e) and I-131 therapy was significantly associated with improved PFS in the favorable group (P = 0.02; Figure 9f). As for OS (Figure 9g–i), I-131 therapy exerted no effects on OS in patients from the unfavorable group (P = 0.88; Figure 9h). In the favorable group, no death occurred in patients with I-131 therapy (P = 0.012; Figure 9i).

Figure 9

Comparisons of PFI, PFS, and OS between PTC patients with and without postoperative I-131 therapy in whole cohort or subgroups stratified by the risk scores. (a–c) Effects of I-131 therapy on PFI in the whole cohort (a), unfavorable subgroup (b), and favorable subgroup (c). (d–f) Effects of I-131 therapy on PFS in the whole cohort (d), unfavorable subgroup (e), and favorable subgroup (f). (g–i) Effects of I-131 therapy on OS in the whole cohort (g), unfavorable subgroup (h), and favorable subgroup (i).

3.6 GSEA

GSEA was conducted to identify the differentially enriched hallmark gene sets between the high-risk (upper tertile) and low-risk (lower tertile) groups. The gene expression profiles of the high-risk and low-risk groups were compared and subjected to GSEA against hallmark gene sets. The results revealed that 21 gene sets were enriched in the high-risk group while none were enriched in the low-risk group based on the cut criteria. The top ten enriched gene sets in the high-risk group were allograft rejection, interferon-gamma response, inflammatory response, interferon-alpha response, epithelial–mesenchymal transition, coagulation, IL6 JAK STAT3 signaling, Kras signaling up, complement, and IL2 STAT5 signaling (Figure 10).

Figure 10

Top enriched hallmark gene sets between high and low-risk groups determined by GSEA: (a) allograft rejection, (b) interferon-gamma response, (c) inflammatory response, (d) interferon-alpha response, (e) epithelial–mesenchymal transition, (f) coagulation, (g) IL6 JAK STAT3 signaling, (h) Kras signaling up, (i) complement, and (j) IL2 STAT5 signaling.

3.7 Exploration of the lncRNAs in PTC at a single-cell level

scRNA-seq data of 22 samples were extracted from the GEO database including six primary PTC tumors, six paired adjacent normal tissues, five initially treated involved LNs, three recurrent LNs, and two subcutaneous metastases. After quality control, a total of 156,295 cells remained for further analyses. Six main cell populations including B cells, endothelial cells, fibroblasts, myeloid cells, T/natural killer cells, and thyroid cells were identified according to corresponding markers (Figure 11a). The six lncRNAs in our signature were mainly expressed in thyroid cells (Figure 11b), then only thyroid cells were selected for subsequent analyses. The proportion of CTA-398F10.2, RP4-794H19.1, RP11-677M14.3 positive cells in tumor and LN samples was lower than that in adjacent normal samples while the proportion of RP11-60L3.1, RP11-20J15.3, and AE000661.50 positive cells was higher (Figure 11c), similar with the results from TCGA bulk RNA sequencing (Figure S1). The thyroid cells could be further clustered into nine subgroups, which included three developmental hierarchies (State 1–3) based on the cell markers derived from trajectory analysis in the previous report [27]. State 1 indicated the normal thyroid cells, state 2 indicated the premalignant cells, and state 3 indicated the malignant cells. Thyroid cells in adjacent normal samples were all with states 1 and 2 whereas cells with state 3 were enriched in the tumor, LN, and distance metastasis samples (Figure 11d). CTA-398F10.2, RP4-794H19.1, and RP11-677M14.3 were mainly expressed in states 1 and 2 cells while RP11-60L3.1, RP11-20J15.3, and AE000661.50 were mainly expressed in state 3 cells (Figure 11e). Furthermore, we performed hallmark pathway enrichment analyses between the thyroid cells with positive or negative LncRNA expression (Figure 12). A series of hallmark pathways were negatively associated with CTA-398F10.2, RP4-794H19.1, and RP11-677M14.3 expression, including angiogenesis, apoptosis, cholesterol homeostasis, coagulation, complement, epithelial–mesenchymal transition, IL2 STAT5 signaling, inflammatory response, interferon-alpha response, interferon-gamma response, Kras signaling up, P53 pathway, and TNFα signaling via NFκB. Most of those pathways were positively associated with RP11-60L3.1 expression. And only several pathways were significantly associated with RP11-20J15.3 and AE000661.50 expression.

Figure 11

Exploration of the lncRNAs in papillary thyroid carcinoma at a single cell level. Single-cell RNA sequencing data of 22 samples from PTC patients were obtained from the GSE184362 dataset. (a) Cells were clustered into six main types and visualized by UMAP plot. (b) Proportions of the cells with different lncRNAs expression in different cell clusters. (c) Proportions of the thyroid cells with different lncRNAs expression in different sample types. P, paired adjacent normal tissues; T, primary tumors; LN, initially treated involved lymph nodes; rLN, recurrent LN; rSC, recurrent subcutaneous metastases. (d) Proportions of the thyroid cells with different state statuses in different sample types. (e) Proportions of the thyroid cells with different lncRNAs expression in different state cells.

Figure 12

Hallmark pathways enrichment analysis between the thyroid cells with or without the corresponding LncRNA expression based on the single-cell RNA sequencing dataset GSE184362. NES, normalized enrichment score.