Construction and validation of a PANoptosis-related lncRNA signature for predicting prognosis and targeted drug response in thyroid cancer

Thyroid cancer (TC) is the most prevalent malignancy of the endocrine system. PANoptosis, a newly discovered cell death pathway, is of interest in tumor research. However, the relationship between PANoptosis-related lncRNAs (PRlncRNAs) and TC remains unclear. The study aimed to develop a prognostic model based on PRlncRNAs in TC. Gene expression data of PANoptosis-associated genes and clinical information on TC from The Cancer Genome Atlas (TCGA) database were analyzed by Pearson correlation analysis, univariate/multivariate Cox analysis, and Lasso Cox regression analysis. A PRlncRNA signature was constructed and used to develop a nomogram to predict overall survival (OS). We further explored the correlation between the risk score and tumor immune microenvironment, immune checkpoints, and drug sensitivity. Moreover, we verified the expression and biological function of lncRNAs in TC cell lines. Finally, seven PRlncRNAs were used to construct a prognostic model for predicting the OS of TC patients. We found that the risk score was associated with the tumor microenvironment (TME) and the expression of critical immune checkpoints. In addition, we screened for drugs that high- or low-risk TC groups might be sensitive to. Quantitative real-time polymerase chain reaction (qRT-PCR) results showed differential expression of four PRlncRNAs (GAPLINC, IDI2-AS1, LINC02154, and RBPMS-AS1) between tumor and normal tissues. Besides, a GEO database (GSE33630) was used to verify the expression differences of PRLncRNAs in THCA tissues and normal tissues. Finally, RBPMS-AS1 was found to inhibit the proliferation and migration of TC cells. In conclusion, we developed a PANoptosis-related lncRNA prognostic risk model that offers a comprehensive understanding of TME status in patients with TC and establishes a foundation for the choice of sensitive medications and immunotherapy.


INTRODUCTION
Thyroid cancer (TC) is the most prevalent malignancy of the endocrine system, and its incidence has rapidly increased worldwide in recent years (Lim et al., 2017).Owing to

Functional analysis of PANoptosis-related genes
Based on previous literature and databases, we identified 35 PANoptosis-related genes (PRGs) (Table S1).The selected PRGs were characterized by Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment utilizing the "clusterprofiler" package of R software (version 4.0.5;R Core Team, 2021).The results were visualized and plotted with the help of "enrichplot" and "ggplot2" packages.

Identification of PANoptosis-related lncRNAs
Gene expression data of 571 TC samples (512 tumor and 59 normal) and matching clinical data of 507 TC patients (Patients with incomplete tracking data were not included) were acquired from The Cancer Genome Atlas (TCGA) databank (https://portal.gdc.cancer.gov/).The expression data is available at https://www.jianguoyun.com/p/DeqNJTAQ9pbeCxjIxY0FIAA, and clinical data are listed in Table S2.To identify PRlncRNAs, Pearson correlation was used to determine the correlation between lncRNAs and PRGs in the present study.Correlation coefficient |R 2 | > 0.4 and p < 0.001 were used as inclusion criteria for PRlncRNAs.The differentially expressed LncRNAs in TC were obtained using the "limma" package.

Establishment of the prognostic model
We randomly distributed all TC patients in the TCGA dataset into the training set (n = 256) and the test set (n = 256) in a 1:1 ratio, and the training set was used for model construction.Data from the test set and entire set were utilized to estimate the performance of the prognostic model.The "Survival" package was applied to conduct univariate/multivariate Cox analysis for PRlncRNAs with p-value ≤ 0.05.Subsequently, we determine the optimal prediction model by means of the least absolute choice operator (LASSO) (Simon et al., 2011).The risk scores of individual patients were computed by incorporating the expression level of lncRNAs in the model and the regression coefficient in the multivariate analysis.The formula is as follows: Σcoefficient of (PRlncRNAi) × expression of (PRlncRNAi).

Validation of model
The samples were categorized as high-or low-risk according to their average risk scores.Kaplan-Meier method was used for survival analysis of high-risk and low-risk groups.The "pheatmap" package was used to map the risk profile, survival status, and lncRNA expression heat map of the groups.The "survivalROC" and "timeROC" packages were used to plot the Receiver Operating Characteristic (ROC) curve and the Area Under the Curve (AUC) to determine the level of accuracy of the signature.Univariate and multivariate Cox regression analyses were performed to verify whether the predictive signature was an independent prognostic variable for Overall Survival (OS).Using the "rms" R package (Iasonos et al., 2008), we combined the OS of TC patients with age, gender, and tumor stage to create a nomogram.A calibration curve was used to test the agreement between the actual results and model predictions according to the Hosmer-Lemeshow test.To compare the accuracy, we calculated the consistency index (cindex) and restricted mean survival (RMS) of the published models (Guo et al., 2022;Qin et al., 2022;Ding et al., 2022;Liu et al., 2022;Shan et al., 2022) and the model we developed.

Estimation of tumor-infiltrating immune cells
THCA immune cell infiltration assessment data were were firstly estimated by the Tumor Immune Estimation Resource (http://timer.cistrome.org/).Based on the TCGA database expression data, we applied the TIMER (Li et al., 2017), CIBERSORT (Newman et al., 2015), CIBERSORT-ABS, QUANTISEQ, MCPcounter (Dienstmann et al., 2019), xCELL (Aran, Hu & Butte, 2017), and EPIC (Racle & Gfeller, 2020) algorithms to obtain the cell composition for the TCGA-THCA dataset samples by R software (Table S3).This enabled comparison of immune cell subpopulations in the groups.ssGSEA was also performed using the "GSVA" packages in R to compare the immune-cell infiltration in the groups.ImmuneScore, StromalScore, and ESTIMATEScore (ImmuneScore + StromalScore) were quantified for each patient using the "limma" and "estimate" packages in R. Finally, the "ggpubr" package in R was used to map immune checkpoints for risk groups (Charoentong et al., 2017).

Clusters based on prognostic PRlncRNAs
Given the PRlncRNA expression associated with the risk signature, we obtained two subtypes using the "ConsensusClusterPlus" package in R. We then evaluated the differences in immune cell infiltration and immune scores between the two subtypes.

Drug sensitivity prediction
The accuracy of the signature in forecasting the response to common antitumor drugs on TC was also evaluated.To determine the therapeutic response, the half inhibitory concentration (IC50) of each sample was estimated with the help of the "pRRophetic" package (Geeleher, Cox & Huang, 2014).Lower IC50 values indicate higher drug sensitivity.The Wilcoxon test was used to calculate the p-values of differences in drug sensitivity.

Quantitative PCR
Ten paired tumor and para-tumor TC tissue samples were acquired from Qilu Hospital of Shandong University.Our study was approved by the Research Ethics Committee of Qilu Hospital of Shandong University (Approval Number: KYLL-2017(KS)-248), and all patients signed a written informed consent form.RNA was isolated from tissues with the total RNA extraction reagent TRIzol (DP424; TIANGEN Biotech, Beijing, China).The obtained RNA was then reverse-transcribed to cDNA using HiScript II Reverse Transcriptase (Vazyme, Nanjing, China).Quantitative PCR assays were performed using a Bio-Rad C1000 Thermal Cycler with SYBR Green Master Mix (Q711; Vazyme Biotech, Nanjing, China), and expression levels were calculated using the 2 −ΔΔCt method.β-actin was used as a standard control.The primers designed for qRT-PCR were obtained from Sangon Biotech (Shanghai, China).And the sequences are presented in Table S4.

Cell culture and transfection
The human thyroid cell lines TPC-1 (RRID: CVCL_6298) was obtained from the Stem Cell Bank, Chinese Academy of Sciences (Shanghai, China).The PTC cell line BHP10-3 (RRID: CVCL_6278) were supplied by Qilu Hospital of Shandong University (Zhang et al., 2019).TPC-1 and BHP10-3 cells were cultured in Roswell Park Memorial Institute (RPMI 1640) complete medium (Shanghai Institute of Biological Sciences, Shanghai, China) at 37 C in an incubator with a constant temperature and 5% CO 2 .The cells were dispersed in six Well Cell Culture Plates and transfected the following day.The overexpression plasmid was obtained from MiaoLing Plasmid Platform (Wuhan, China).When the cell density reached approximately 70%, 2,000 ng of negative control or plasmid was transfected into cells using Lipofectamine TM 2000 Transfection Reagent (11668019; Invitrogen, Waltham, MA, USA).After 6 h, the medium was replaced with 2 mL of complete medium.After 24 h of transfection, cells were collected for subsequent experiments.

Proliferation and migration assay
After 24 h of transfection, the cells were evenly distributed in 96-well plates, five parallel replicate wells were set for each group, and 10 mL of CCK-8 reagent was added on days 1, 2, and 3.The absorbance was measured at 450 and 630 nm, which represents the cell growth rate.Thyroid cancer cells (5 × 10 4 ) transfected with RBPMS-AS1 overexpression plasmid or negative control were introduced to the upper chamber of Transwell plates in 200 mL serum-free culture medium, and 800 mL medium comprising 10% fetal bovine serum was introduced into the lower chamber as an inducer.After 24 h, the migrated cells were washed with PBS, fixed in methanol, stained with 1% crystal violet, and non-migrating cells in the upper chamber were removed with a cotton swab.At least four randomly selected fields were observed and counted under a microscope (Olympus, Tokyo, Japan).

Wound healing assay
The cells were incubated in 6-well plates.When the cells were cultured to 80% confluence, a 200 mL pipette tip was used to create a scratch.The cells were then washed with PBS and placed in serum-deficient medium for an additional 24 h.The wound healing area at 0 and 24 h was recorded using a microscope (Olympus, Tokyo, Japan).

Functional enrichment analysis of PANoptosis-related genes
From the previous literature, we obtained 35 PRGs.Based on GO analysis, PRGs primarily participate in biological processes (BP), such as the positive regulation of cytokine production, positive regulation of response to external stimulus, and response to virus.GO analysis of cellular components (CC) revealed that PRGs were mainly positioned in the membrane raft, membrane microdomain, and inflammasome complex.According to molecular function (MF), PRGs mostly carried out the cytokine receptor binding function (Figs.1A, 1B).KEGG enrichment analysis showed that the selected PRGs were mainly involved in signaling pathways including the NOD-like receptor signaling pathway, Influenza A-mediated pathways, and necroptosis (Figs.1C, 1D).

Evaluating the predictive capacity of the risk signature
To examine the efficacy of the model, the patients in the public database were assigned to training and validation cohorts.The basic information of the patients in both sets is listed in Table 1.Based on the risk score, the patients in each set were classified into high-risk and low-risk groups.Patients in the high-risk group had poor survival according to the Kaplan-Meier curve analysis (p < 0.001) (Fig. 3A).Figures 4B and 4C show the survival status, distribution of patients, and the mortality rate increased with risk scores.Both testing and the entire set suggested a poor prognosis in high-risk patients (Figs.3F-3H, 3K-3M).Furthermore, ROC was used to validate the risk signature.The AUC values of this risk signature were 1.000, 0.986, and 0.882 for 1, 3, and 5 years respectively, suggesting  good predictive power (Fig. 3D).In addition, the ROC analysis results had high accuracy for both the testing set and the entire set (Figs. 3I, 3M).Figures 3E, 3J and 3O depict the heat maps of the expression levels of the seven PRlncRNAs in all sets.

Association of clinicopathological factors with PRlncRNAs risk signature
To verify the predictive effect on patients with different clinicopathologic variables, patients with TC were stratified in terms of age, gender, and clinical stage.Patients in the high-risk group had notably shorter OS in all subgroups of age (≤65 years), gender (female or male), TNM stage (T1-2 or T3-4, N0 or N1, M0), and clinical stage (stage-III-IV) (Fig. 4).Moreover, the model demonstrated the highest sensitivity and specificity compared to other clinical variables such as age, gender, grade, and stage in all sets (Figs. 5A-5C).To better assess the independent forecasting capacity of the signature, the Univariate/Multivariate Cox regression analyses were conducted on risk scores and patient clinical data.Univariate Cox regression analysis revealed that age, stage, and risk scores acted as independent prognostic indicators (p < 0.001; Fig. 5D), while multivariate independent prognostic analysis showed that only age and risk scores acted as independent prognostic factors for TC patients (p < 0.001; Fig. 6E).We then incorporated risk scores with clinical factors and plotted a predictive nomogram for patients with TC (Fig. 5F).
Combined with the calibration plots, it was observed that the predicted OS was consistent with the actual observations (Fig. 5G).In addition, a comparison of the c-index and the RMS of seven previously reported risk models shows that the developed model has a significant advantage (Figs.5H, 5I).Therefore, the risk signature showed high performance in predicting OS in patients with TC in the TCGA dataset.

Discovery of the association between tumor-infiltrating immune cells and risk signature
We explored the underlying differences in signaling pathways between the risk groups by Gene Set Enrichment Analysis (GSEA).GSEA of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways showed that the enriched pathways in the high-risk group involved ABC transporters, ECM-receptor interaction pathway, hedgehog signaling athway, primary bile acid biosynthesis, and tryptophan metabolism (Fig. 6A).Conversely, the accumulated pathways in the group at low risk included DNA replication, Nucleotide excision repair, Ribosome, and Spliceosome (Fig. 6B).
To elucidate the relationship between the tumor immune microenvironment (TME) and the risk signature, we performed an analysis of RNA microarray data to explore the correlation between risk characteristics and immune cells.The bar graph shows the association between risk score and tumor-infiltrating immune cells (Fig. 6C).The XCELL algorithm results showed that the risk score was significantly positively correlated with B cell plasma, plasmacytoid dendritic cells, and T cell regulatory cells (Tregs), while it was positively correlated with T cell CD4 + Th1, CD8 + , and mast cells.The results of the QUANTISEQ algorithm showed that the risk score was significantly positively correlated with regulatory T cells (Tregs), while negatively correlated with macrophage M1, neutrophils, etc.The EPIC algorithm showed that the risk score was positively correlated with cancer-associated fibroblasts and negatively correlated with macrophages and NK cells.To better explore the immune status between the two groups, we evaluated the stromal score, the immune score, and the estimated score (Fig. 6D).Subsequently, we investigated the variations between the risk points and tumor infiltrating immune cells using ssGSEA.The scores for B cells, macrophages, and T helper cells were markedly higher in the group at high risk (Fig. 6E).Furthermore, we observed that there was a considerable difference (p < 0.001) between the groups in immune checkpoints such as IDO2, CD44, CD40, HHLA2, ADOARA2A, and CD80 (Fig. 6F).Cluster analysis according to prognostic PRlncRNAs In an effort to evaluate the immune microenvironment and response of different subtypes, cluster analysis was performed to obtain clustered subtypes.On the basis of the seven PRlncRNAs, we divided the patients into two subtypes by applying the ConsensusClusterPlus package (Fig. 7A).As shown in the Sankey diagram (Fig. 7B), most of the low-risk patients were grouped in cluster 1, whereas the majority of the high-risk patients were assigned to cluster 2. Survival analysis showed that patients in subgroup 1 had a superior OS (p < 0.001) (Fig. 7C).Based on the cluster analysis, the relationship between the different clusters and TME were examined.Box plots show that cluster 1 had higher immune, stromal, and ESTIMATE scores than cluster 2 (Fig. 7D).The heat map shows the variation in immune cell infiltration across different tumor clusters (Fig. 7E).
Drug sensitivity in the PANoptosis-related LncRNA signature We further explored the differences in drug resistance between the risk groups.
We investigated the IC50 values of 80 chemotherapy drugs and inhibitors in both groups.
The values of six typical drugs are shown in Fig. 8. Low-risk patients had substantially lower IC50 values for Imatinib, Sunitinib, Pazopanib, Mitomycin C, 5-Fluorouracil, and tipifarnib, suggesting that the drugs may be eligible for the treatment of patients in the low-risk group.In contrast, drugs such as Bortezomib and Phenformin showed better effects in the high-risk group.This means that candidate immunotherapeutic candidates for patients can be identified with TC based on their risk scores.

Assessing the expression and biological function of PRLncRNAs
Among the PRLncRNAs, we verified the gene levels of four lncRNAs in pairs of samples obtained from patients with TC in our center.Comparable expression patterns were observed in the clinical samples (Figs.9A-9D).GAPLINC showed increased expression in tumor tissues (T) compared to in normal tissues (N), whereas RBPMS-AS1, IDI2-AS1, and LINC02154 showed lower expression levels.In addition, we also used a GEO database (GSE33630) to verify the expression differences of PRLncRNAs between thyroid cancer and normal tissues.However, we were able to obtain only the expression of RBPMS-AS1, IDI2-AS1, and GAPLINC because of limited data from the sequencing platform (Fig. 9E).
As shown in Fig. S2, the RNA expression of RBPMS-AS1 was successfully upregulated in both cell lines after transfection.Functionally, its role in regulating TC cell proliferation was tested using the Cell Counting Kit-8 (CCK-8) assay (Fig. 9F).Overexpression of RBPMS-AS1 inhibited the proliferation of BHP10-3 cells significantly but did not affect TPC-1 cell growth (Fig. 9E).By wound healing assay and transwell experiments, the overexpression of RBPMS-AS1 decreased the migration ability of the two TC cell lines significantly (Figs.9G, 9H).

DISCUSSION
Thyroid cancer is the most prevalent malignancy of the endocrine system and is highly heterogeneous in terms of morphological features and prognosis (Bray et al., 2018).
In clinical practice, the diagnosis and treatment of TC has shown promising efficacy.However, owing to the high rate of postoperative recurrence, patients still do not have a satisfactory prognosis.Thus, there is a great need to explore new biomarkers for screening high-risk groups and individual treatments.PANoptosis is a newly recognized mode of inflammatory cell death characterized by the interaction of apoptosis, pyroptosis, and necroptosis and has been shown to be present in a variety of diseases (Wang & Kanneganti, 2021).In 2020, Karki et al. (2020) found that activation of IRF1-dependent "PANoptosis" in a mouse model of colorectal cancer could prevent AOM/DSS-induced colorectal carcinogenesis.In addition, accumulating evidence suggests that a broad range of programmed cell death processes, including apoptosis, pyroptosis, and necroptosis, play an important role in TC (Holm et al., 2022).Most of the research on PANoptosis has focused on molecular mechanisms.However, few researchers have investigated the prognostic value of PRGs in tumors.Therefore, there is an urgent need to construct prognostic models based on the crucial lncRNAs of TC to guide prognostic treatment and management.
In this study, RNA-seq data of PRLncRNAs and related clinical data were obtained from TCGA database.First, we performed KEGG and GO enrichment analyses using the 35 PRGs that were mainly involved in apoptosis, pyroptosis, and necroptosis pathways.By Cox regression analyses, we identified seven lncRNAs for building a PANoptosisrelated lncRNA risk model.Based on the risk scores, we stratified the TC patients into high-and low-risk groups.Survival analysis verified a statistically meaningful prognostic difference between groups (p < 0.001).Univariate and multivariate Cox regression analyses showed that the model could independently predict OS in THCA (p < 0.001).According to the ROC curves, the signature was more accurate and reliable than other clinical characteristics, and the AUC value of the test set at 1 year was 1.000.Compared with previously reported RNA or lncRNA risk models for TC, our signature is more powerful according to the c-index and RMS.
In recent years, with advancements in genomic technologies, many noncoding genes have been identified to have key roles in tumorigenesis and evolution.Among the seven PRlncRNAs, five lncRNAs (IDI2-AS1, LINC02154, C5orf34-AS1, LINC01705, and GAPLINC) were poor prognostic factors for TC, while the other two lncRNAs (DPP4-DT and RBPMS-AS1) were protective factors.Yue et al. (2022) demonstrated that LINC02154 has an imperative effect on the progression of hepatocellular carcinoma by increasing the activity of the SPC24 promoter and triggering the PI3K-AKT signaling pathway (Yue et al., 2022).According to Du et al. (2020), LINC01705 is a crucial oncogene that prevents apoptosis and accelerates proliferation of breast cancer cells (Du et al., 2020).Hu et al. (2021) revealed that DPP4 gene knockdown inhibits proliferation and epithelialmesenchymal transition in PTC cells through suppressing the MAPK pathway (Hu et al., 2021).RBPMS-AS1 enhances radiosensitivity and apoptosis, which facilitate regulation of GBM cell proliferation (Li et al., 2022).Xu et al. (2022) reported that IDI2-AS1 is a potential prognostic marker for uveal melanoma.The results are in line with our own findings and further demonstrate their reliability.
To further explore the role of the seven PRlncRNAs in TC, we compared tumor-infiltrating immune cells between the groups.The results showed that patients in the group at high risk had a significantly higher level of immune cells, such as B cells, macrophages, and T helper cells.A previous study showed that macrophages promote tumor growth and metastasis by secreting various cytokines (Vlaicu et al., 2013).B-cell infiltration has diagnostic and prognostic value for PTC lymph node metastases (Yang et al., 2021).Meanwhile, we observed a notable decrease in the immune, stromal, and ESTIMATE scores in the high-risk group.Therefore, these results suggest that the signature has predictive capacity to reflect the immune infiltration of TC and may provide a reference for immunotherapy.
Moreover, we verified the expression of four lncRNAs (GAPLINC, IDI2-AS1, LINC02154, and RBPMS-AS1) in tumor and adjacent normal tissues, which was consistent with the expression of PRlncRNAs in the database.Of note, RBPMS-AS1, as a protective factor, was expressed less in TC.We then overexpressed RBPMS-AS1 in two types of TC cell lines and found that their migration and proliferation abilities were largely inhibited.The expression and biological function of PRlncRNAs were in line with our previous analysis.
In summary, this is the first study to use bioinformatics to construct a prognostic risk model for TC using PRlncRNAs.Our work may provide predictive tools for immune and targeted therapy for patients with TC.Although our PRlncRNA risk model achieved promising results in THCA, it still has some limitations.Our signature was established using the TCGA database only, lacking further verification of external datasets.The predictive effect of the model is yet to be validated in a large cohort of thyroid cancer patients.In addition, bioinformatics analysis is the main basis of our study, and the mechanism by which PRlncRNAs affect TC progression is uncertain.Therefore, our study provides new insights into the possible designs of PANoptosis-related drugs to treat thyroid cancer.The molecular mechanisms underlying the role of PRlncRNAs in thyroid cancer remain to be further elucidated in the future.

CONCLUSION
In conclusion, this analysis was the first to identify seven PANoptosis-related lncRNAs (IDI2-AS1, LINC02154, DPP4-DT, RBPMS-AS1, C5orf34-AS1, LINC01705, and GAPLINC) in TC.Furthermore, the constructed model can accurately predict OS of TC and guide targeted drugs, which holds promise for predicting the immune infiltration of tumors and guiding targeted drug therapy.

Figure 5
Figure 5 Independent prognostic factor assessment and nomogram construction.(A-C) Area under the curve (AUC) of receiver operating characteristic (ROC) curves comparing the predictive accuracy of risk score and other prognostic factors in training, validation, and overall groups.(D, E) Univariate and multifactorial Cox analyses of clinicopathologic factors and risk score with overall survival (OS).(F) The nomogram for the prediction of 1-, 3-, and 5-year OS. (G) The calibration plots of the 1-, 3-and 5-year nomograms.(H) The comparison of C-index among risk models.(I) The comparison of restricted mean survival (RMS) among risk models.Full-size  DOI: 10.7717/peerj.15884/fig-5

Figure 6
Figure 6 Differences in tumor immune microenvironment between low and high-risk groups.(A, B) Gene set enrichment analysis of KEGG in high-risk and low-risk groups.(C) The immune cell bubble of risk groups.(D) The box plots comparing StromalScore, ImmuneScore, and ESTIMATEScore between low and high-risk groups.(E) The single-sample gene set enrichment analysis (ssGSEA) scores of immune cells.(F) The difference of common immune checkpoint expression in the risk groups.Ã p < 0.05, ÃÃ p < 0.01, ÃÃÃ p < 0.001.Full-size  DOI: 10.7717/peerj.15884/fig-6

Table 1
The basic information of TCGA-THCA.