We extracted the expression data and the matched clinical information of patients with papillary thyroid carcinoma (PTC) from The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/). Then, we classified the RNAs into either protein-coding function or lncRNA using the Ensemble human genome browser (19). The expression data for all the messenger RNAs (mRNAs) and lncRNAs in the cohort were log2 transformed [log2(FPKM+1)] for downstream analyses. A total of 508 patients whose pathological subtype is PTC were included in the study, and those whose survival time was less than 30 days were excluded because their death is likely to be classified under surgery-related death.

PTC cell lines BCPAP, TPC‐1, and K1 were used in the following experiment. TPC-1 and K1 were cultured in Roswell Park Memorial Institute‐1640 medium (Gibco, USA) containing 10% fetal bovine serum (FBS) (Gibco, USA) and 1% antibiotics (P/S) (Gibco, NY, USA), and BCPAP was cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco, USA) containing 10% FBS (Gibco, USA) and 1% antibiotics. All samples were placed in an incubator (Thermo, USA) with 5% CO₂ at 37°C. Cells were dissociated with 0.25% trypsin at a ratio of 1:3 when their density reached 80%.

The small interfering RNAs (siRNA) were utilized for knocking down CRNDE, and scrambled siRNA (si‐NC) was used as the negative control. The sequences were shown in Supplementary Table S1 . All of them were synthesized by RiboBio (Shanghai, China). Cell transfection was performed with Lipofectamine 2000 (Invitrogen, USA) according to the product instructions. After harvesting for 48 h, the efficiency of knocking down was identified using qRT-PCR.

Cell proliferation was measured using the Cell Counting Kit-8 (CCK-8) Assay (MedChemExpress, China) according to the manufacturer’s instructions. Briefly, PTC cells with and without transfection (4 × 10³ cells/well) were cultured in a 96-well plate in the presence or absence of sorafenib (MedChemExpress, China) (2 µM), and the cell viability was observed by measuring absorbance at 450 nm (24, 48, and 72 h) in a microplate reader after incubation with CCK-8 solution (10 µl) for 2 h.

TRIzol (Invitrogen, CA) was used to extract total RNA from PTC cells. Here, 1 ml of TRIzol reagent was added to a 3.5-cm dish and incubated for 10 min. Following this step, 0.2 ml of chloroform was added to the dish. All the contents of the dish were transferred into tubes and shaken vigorously for 16 s. After standing in ice for 10 min, the tubes were centrifuged at 12,000 rpm for 25 min at 4°C. Finally, the RNA was washed with 75% ethyl alcohol and dissolved in diethylpyrocarbonate (DEPC) water. Then, 1 µg extracted RNA retrieved in the above step was reverse-transcribed into cDNA in a 20-µl reaction volume by using the First Strand cDNA Synthesis Kit (Takara, Japan) according to the manufacturers’ instructions. The cDNA was used as a template for real-time PCR by using matched primers. The sequences used in quantitative PCR were shown in Supplementary Table S1 . The length of amplified products was within 300 bp. The real-time PCR was performed with a 20-µl reaction system containing 2 µl cDNA, 2 µl of mix of primer, 4 µl dd H₂O, and 10 µl SYBR Green Real-Time PCR Master Mix (Takara, Japan). The detailed reaction conditions were set by referring to the manufacturers’ instructions. Data were analyzed by using the ΔΔCt method, where the endogenous housekeeping gene β-actin was used as a quantity and quality control.

A total of 232 autophagy genes were obtained from the Human Autophagy Database (HADb; http://www.autophagy.Lu), which provides a comprehensive and up-to-date list of human genes involved in autophagy and is shown in Supplementary Table S2 . The Pearson correlation coefficient was calculated between the expression of autophagy genes and autophagy-related lncRNAs. Autophagy-related lncRNAs were selected based on the absolute value of the coefficient being >0.3 (|R|>0.3), and a P value <0.001 was considered significant. Detailed information on the correlation between lncRNAs and autophagy-related genes is shown in Supplementary Table S3 .

To further explore the interaction between lncRNAs and autophagy genes, we established an lncRNA and mRNA network using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING). The software Cytoscape 3.8.1 was used to visualize the relationship between lncRNAs and mRNAs.

Univariate and multivariate Cox regression analyses were performed to identify potential autophagy-related lncRNAs with prognostic value. The autophagy lncRNAs with a P value <0.01 were submitted to the next step to perform the multivariate Cox regression analysis. The risk stratification based on the expression of autophagy-related lncRNAs was established by using the candidate lncRNAs obtained from the multivariate Cox regression. The linear risk formula was RiskScore = (Coef1 × Expression Gene1) + (Coef2 × Expression Gene2) + (Coef3 × Expression Gene 3) + (Coef4 × Expression Gene4) + (Coef5 × Expression Gene5) +……+ (Coef Gene N × Expression Gene N).

Based on the expression profiles of autophagy-related lncRNAs, we conducted a principal component analysis (PCA) to investigate the difference between the high-risk group and low-risk group.

Gene set enrichment analysis (GSEA) was performed using GSEA software (Version 4.1.0). All operations were performed as previously described (19–22). The results of GSEA were visualized using an enrichment map. Each analysis of gene set permutation was calculated 1,000 times, and pathways enriched in specific phenotypes were sorted by normalized enrichment score (NES) and P value. Relevant gene sets were obtained from the gene set database of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways detailed in Supplementary Figure S2 .

A Cox regression model was employed to construct the risk stratification. The hazard ratio (HR) and associated 95% confidence interval (CI) were calculated. From the analysis of two subgroups based on the median value of risk score, the overall survival (OS) times of the high-risk and low-risk subgroups were calculated and analyzed via Kaplan–Meier test and compared using the log-rank test. The time-dependent receiver operating characteristic (ROC) curve was used to evaluate the predictive accuracy of the autophagy-related lncRNAs.

By constructing the coexpression network of the 232 autophagy genes, a total of 1,283 autophagy-related lncRNAs were obtained ( Figures 1A–D ). Then, we performed Cox proportional hazards analysis to obtain 144 autophagy lncRNAs with significant prognostic value by analyzing the expression profile ( Figure 1 ). There were 32 lncRNAs with low risk (HR <1) and 112 lncRNAs with high risk (HR >1) ( Supplementary Table S4 ). Following this, multivariate Cox analysis was performed to further screen 10 lncRNAs from 144 autophagy lncRNAs with prognostic value, and they were AC008063.1, AC011297.1, FAM201A, AC092279.1, LINC00900, AL162231.2, CRNDE, TONSL-AS1, LINC02454, and AC004918.3 ( Figures 4A–J ) ( Supplementary Table S5 ).

Then, these optional lncRNAs were submitted to the next step to construct the prognostic risk stratification of autophagy-related lncRNAs. We established a coexpression network of autophagy-related lncRNAs with prognostic value and mRNAs, as shown in ( Figures 1A, C ). And the formula is as follows:

Risk score = (AC008063.1*-1.740157523) + (AC011297.1*0.171231624) + (FAM201A*0.895715416) + (AC092279.1*-3.730086976) + (LINC00900*-2.38496601) + (AL162231.2*0.671649156) + (CRNDE*0.450612252) + (TONSL-AS1*-0.79310585) + (LINC02454*0.588168097) + (AC004918.3*2.412626944).

Based on the risk formula and median risk score, the PTC patients were divided into a high-risk group and a low-risk group ( Figure 2A ). Kaplan–Meier survival analysis showed that patients in the high-risk group had shorter OS than that in the low-risk group ( Figure 2B ), which suggested that the risk stratification had good prognostic prediction ability.

A scatter plot of survival status and risk curve were made in order to display the risk score and the corresponding survival status of PTC patients, the results of which showed that the higher risk score corresponded to higher mortality. At the same time, the heatmap was drawn based on these 10 autophagy-related lncRNAs in the PTC samples ( Figure 2A ). The survival curve ( Figures 2B–D ). The heatmaps displayed the expression of these elements in the high-risk group and low-risk group. Hence, these results indicate that all 10 autophagy lncRNAs have good prognostic value for PTC ( Figures 2 , 3 ).

Univariate and multivariate Cox regression analyses were employed to determine whether the risk stratification based on autophagy-related lncRNAs was an independent prognostic factor of PTC. The HR and 95% CI of the risk score based on the risk stratification were 1.003014104–1.00616136988199 (P < 0.05) in the univariate Cox regression analysis and 1.004196324–1.00559955694039 (P < 0.05) in the multivariate Cox regression analysis ( Figures 3D, E ). The results showed that the risk stratification based on autophagy-related lncRNAs was a powerful significant prognostic factor of PTC, independent of clinicopathological characteristics such as tumor size, metastasis of lymph node and distant metastasis, sex, and TNM stage. The ROC curve was plotted to evaluate and estimate the predictive specificity and sensitivity of the risk stratification based on the autophagy-related lncRNAs ( Figure 8 ). The time-dependent area under the curve (AUC) for 1, 3, and 5 years was 0.981, 0.901, and 0.963, respectively ( Figure 3B ). At the same time, the AUC of the other clinical characteristic including T, N, M, age, and TNM (stage) was calculated, and the AUC value of the risk score exceeded most of the others ( Figure 3C ). These results indicate that the risk stratification based on autophagy-related lncRNAs was an effective independent prognostic factor for PTC patients. Subgroup analysis was performed based on the risk stratification ( Figures 5C–N ). Besides, those patients both with high tumor mutation burden (TMB) levels and high risk exhibited a poorer survival state ( Figures 5A, B ). PCA results show that the model has an excellent ability to distinguish high-risk and low-risk patients ( Figures 6A–D ). DNA stemness score (DNAss) and RNA stemness score (RNAss) presented significant differences among the two subgroups ( Figures 6E, F ). Different risk types possessed different TIME statuses, Immunogenomics features ( Figure 7 ).

Interestingly, most of the risk genes involved in the stratification system with higher expression in PTC samples matched better prognosis ( Figures 4F , 8R ). After further validation, we found that those patients with higher CRNDE expression matched worse survival state.

We tried to evaluate the association between the risk stratification and the efficacy of chemotherapeutics based on the THCA cohort of TCGA project. Interestingly, the results showed that higher risk score based on the risk stratification was associated with higher half inhibitory concentration (IC50) of sorafenib (P < 0.01) ( Figure 8K ).

Next, we chose lncRNA CRNDE from the 10 autophagy lncRNAs and investigated if it could mediate cancer cell survival under sorafenib challenge. By directly knocking down CRNDE in the PTC cell line (Bacap-1, TPC-1, and K1), CCK-8 assay showed that knocking down CRNDE partially increased the effects of sorafenib-induced cell death in PTC cells with si-CRNDE transfection ( Figures 8I–III ). All those results suggested that CRNDE plays a positive role in tumorigenesis and may regulate the sorafenib sensitivity of PTC cells ( Figures 8IV–VI ).