Cuproptosis-related lncRNA SNHG16 as a biomarker for the diagnosis and prognosis of head and neck squamous cell carcinoma

Data preprocessing and filtering

RNA sequencing data, clinical data, and gene mutation data of 502 HNSCC and 44 normal head and neck tissue specimens were obtained from the Genomic Data Commons portal (https://portal.gdc.cancer.gov/) (Jensen et al., 2017). A total of 2,876 lncRNAs in the TCGA cohort were screened according to gene annotation by using 19 cuproptosis-related genes retrieved from previous studies (Table S1). The coexpression correlation between lncRNA and cuproptosis-related genes was analyzed with limma software package was used to analyze the at |R| > 0.4 and P < 0.001 to determine potential cuproptosis-related lncRNAs. A total of 197 cuproptosis-related lncRNAs were identified, and 111 cuproptosis-related lncRNAs with prognostic value were selected through univariate Cox regression analysis on overall survival (OS) with P < 0.05 considered as significant.

Construction of the prognosis model of cuproptosis-related lncRNAs

A prognostic model was developed using cuproptosis-related lncRNA expression data from a training set and tested on another set. LncRNAs’ correlation with survival was significant (P < 0.05) (Li, Bai & Gao, 2023). Eleven crucial lncRNAs were identified using Lasso and Cox regressions (Baliakas et al., 2023). Risk scores determined patient categorization into low or high-risk groups, with Kaplan-Meier curves depicting survival (Bentzen & Vogelius, 2023). ROC and C-index curves validated the model’s accuracy (Yang & Smith, 2023), and a nomogram predicted overall survival (Kang et al., 2023). PCA was constructed to explore the distribution of patients with different risk scores (Khamis & Alanazi, 2023). In addition, the main biological properties of cuproptosis-related lncRNAs were comprehensively examined by Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Gao et al., 2023).

Immunoassay and tumor mutation burden analysis and drug screening

Differences in immune-related function between the high- and low-risk group based on the cuproptosis-related lncRNA characteristics were evaluated using limma and GSVA packages, with P < 0.05 considered as statistically significant. The relationship between high- and low-risk groups and tumor mutation burden (TMB) was compared using Maftools package (Huang et al., 2023b). In addition, the tumor immune dysfunction and exclusion (TIDE) of HNSCC was downloaded from http://tide.dfci.harvard.edu/, and the relationship between the two groups and TMB was analyzed using ggpubr and limma packages.

Verification of screened lncRNAs by qRT-PCR

Human bronchial epithelial (HBE) cell line and laryngeal squamous carcinoma cell lines TU212 and TU686 were obtained from Shanghai Zhong Qiao Xin Zhou Biotechnology. All cells were maintained at 37 °C in a 5% CO₂ incubator in RPMI-1640 medium with 10% fetal bovine serum and 1% penicillin-streptomycin. Meanwhile, 12 pairs of fresh tissue samples were obtained from patients who had recently undergone surgical treatment in the Department of Otorhinolaryngology, Head and Neck Surgery, Zhongnan Hospital of Wuhan University from October 2021 to October 2022 (all cases had no other acute and chronic diseases; No prior radiotherapy and/or chemotherapy treatment). The process of collecting fresh tissue samples was as follows: following the basic principles of medical ethics and after obtaining the consent of patients pathologically diagnosed with laryngeal squamous cell carcinoma, the middle and two small tissue samples (about 100 and 50 mg, respectively) were taken from adjacent tissues (5 mm away from the lesion). The blood on the surface was washed with normal saline and transferred to 1.5 ml cryopreservation tubes. The frozen tubes were immediately placed in the liquid nitrogen tank and brought back to the laboratory for storage at 80 °C ultralow temperature refrigerator before being taken out for total RNA extraction. Postoperative pathologically confirmed laryngeal carcinoma tissues and peritumoral tissue specimens were qualified. All the basic and clinical data of the patient were obtained, and T, N, M staging and clinical staging were evaluated according to the patient’s specialist examination, auxiliary examination results, and the “International Union for Anti-Cancer TNM Classification and Staging Standards” (for specific clinical and pathological data, see Table S2). The research protocol was approved by the Ethics Committee of Zhongnan Hospital of Wuhan University (approval number: 2021058), and all patients signed the informed consent. The qPCR primer sequences are listed in Table S3, and β-actin was used as a housekeeping gene.

Screening of cuproptosis-related lncRNAs and construction of a prognostic model

Figure 1 depicts the workflow for establishing the prognostic model and subsequent analysis. Our study comprehensively investigated the pivotal roles and prognostic implications of cuproptosis-related lncRNAs in HNSCC. Initially, we collected mRNA and lncRNA expression profiles and clinical information from 502 HNSCC tissue samples and 44 nontumor samples. After analyzing the expression values of 19 cuproptosis genes and 16,876 lncRNAs, we utilized the rigorous screening criteria of |R| > 0.4 and P < 0.001 and identified 783 cuproptosis-related lncRNAs. We then visually represented the coexpression relationships between cuproptosis genes and cuproptosis-related lncRNAs using a Sankey diagram (Fig. 2A). We integrated the expression data of cuproptosis-related lncRNAs with the corresponding clinical information and then randomly divided the data into training and testing sets. Our aim was to construct a prognostic model for patient outcomes using the training set. As an initial step, we conducted univariate Cox regression analysis to pinpoint the lncRNAs significantly associated with survival. We found 35 lncRNAs that were markedly correlated with survival rates (P < 0.05) (Fig. S1A and Table S4). We employed Lasso regression to further refine our model and mitigate the risk of overfitting, which could introduce false positive parameters. We identified the Lasso coefficients for 20 lncRNAs determined by the optimal lambda value as illustrated in Fig. 2B. The best parameter (lambda) was then visualized with vertical lines in Fig. 2C and Table S5. Multivariate Cox regression was utilized to establish the final model. Our research highlighted 11 lncRNAs as independent prognostic factors for HNSCC. These 11 differentially expressed lncRNAs can act as standalone prognostic predictors for HNSCC (Table S6). Finally, Fig. 2D presents a correlation heatmap, offering supplementary evidence of the association between the mutation-related genes and lncRNAs integrated into the prognostic model.

Validation of HNSCC-specific predictive prognostic model

After the prognostic model was established, the patients were divided into high-risk and low-risk groups based on the median risk score. K-M survival analysis was conducted between the training and testing groups and for all groups. The results indicated that OS and PFS were significantly shorter in the high-risk group than in the low-risk group (Figs. 2E and 2F). Moreover, the OS of the high-risk group was lower in the training group (Fig. S1B). To further analyze the survival distribution, we ranked all the patients with HNSCC in the training group according to the risk score (Fig. 2G). A scatter plot illustrated the survival status of patients with different risk scores, with the mortality rate increasing with the risk scores (Fig. 2H). Figure 2I shows the heatmap of the level of 11 lncRNAs in the training group. For example, SNHG16, AC012184.3, LINC01012, AC243965.2, and AP000866.1 are high-risk lncRNAs, and ECE1-AS1, AL358777.3, AL033527.2, AL162458.1, AC090587.1, and CDKN2A-DT are low-risk lncRNAs. We then verified the findings in the testing group, and the results showed that the OS of the high-risk group was significantly shorter than that of the low-risk group (Fig. 3A). In addition, the mortality of the testing group increased with the risk score (Figs. S1C and S1D). The heat map also showed that the expression of cuproptosis-related lncRNAs was associated with the risk score (Fig. 3B). The results for the testing group suggest the reliability of the constructed prognostic model. To determine whether our constructed signal could be used as an independent prognostic factor and is independent of other clinical features, we performed Cox regression analysis on prognostic factors. As shown in Fig. 3C, univariate Cox regression analysis revealed a significant association between OS and age, histological grade, pathological stage, and risk score of patients with HNSCC (P < 0.05). However, upon conducting multiple regression analysis, we observed that only age, N stage, and risk score were the independent prognostic factors significantly associated with OS (P < 0.05) (Fig. 3D). We then employed the ROC curve to assess the predictive accuracy of the risk score, yielding an area under the curve (AUC) of 0.708. The AUC of the risk score was higher than those of age (0.581), gender (0.50), grade (0.544), and stage (0.560) (Fig. 3E). In addition, the AUCs at 1, 3, and 5 years were 0.708, 0.725, and 0.632, respectively (Fig. 3F), signifying that our model was superior in forecasting the prognosis of HNSCC compared with the conventional clinicopathological features.

Construction of prediction nomogram and principal component analysis

A C-index curve was constructed to compare the consistency index of the risk score with other clinical features such as age, gender, and stage. The results showed that the C-index value of the risk score was higher than those of the other clinical features (Fig. 3G). A nomogram was also developed using age, sex, stage, T, risk score, and N in the signature and accurately predicted the 1-, 3-, and 5-year survival of the patients (Figs. 3H and 3I). The patients in the low-risk group with stages I–II and III–IV exhibited better OS than those in the high-risk group (Figs. 4A and 4B) over time. PCA indicated a high degree of differentiation between the high- and low-risk groups for all genes (Fig. 4C), cuproptosis-related genes (Fig. 4D), cuproptosis-related lncRNAs (Fig. 4E), and cuproptosis-related lncRNA prognostic markers (Fig. 4F).

Functional enrichment analysis and immune-related function analysis

Enrichment analysis revealed significant GO and KEGG pathways for cuproptosis-related lncRNAs. BP analysis showed enrichment in immunoglobulin production, B cell activation regulation, and humoral immune response. MF analysis showed enrichment in antigen binding, immunoglobulin receptor binding, and immunoglobulin binding. Cellular component analysis showed enrichment in immunoglobulin complex, T cell receptor complex, and blood microparticles. KEGG analysis showed enrichment in cytokine–cytokine receptor interaction, cell adhesion molecules, B cell receptor signaling pathway, and NF-kappa B signaling pathway. Furthermore, immune cell correlation analysis showed significant differences in immune cell infiltration between high-risk and low-risk groups. Figures 4G–4I indicate that the cuproptosis-related lncRNAs are associated with immune cell infiltration in HNSCC.

Tumor mutation burden and drug sensitivity analysis

Figures 5A and 5B shows that the mutation frequency was generally higher in the high-risk group than in the low-risk group (e.g., TP53: high-risk, 76%; low-risk, 57%. FAT1: high-risk, 24%; low-risk, 16%. CSMD3: high-risk, 22%; low-risk, 16%). We also found that the TMB was significantly higher in the high-risk group than in the low-risk group (Fig. 5C). TIDE algorithm revealed higher immunotherapy resistance in the high-risk group than in the low-risk group (Fig. 5D). Figure 5E shows that high TMB was associated with low OS in patients. A combined survival analysis of TMB and risk score further demonstrated their significant impact on the OS in patients with HNSCC (Fig. 5F). Finally, potential antitumor drugs were screened using the pRRophetic software package. Analysis of the relationship between risk score and drug IC50 in HNSCC treatment revealed eight drugs with significant sensitivity differences between the high- and low-risk groups. Figures 6A–6C indicate that three drugs (PD-0332991, cytarabine, and WH-4-023) had low IC50s in the high-risk group, implying high sensitivity to drug treatment. These drugs may have a potential role in the future treatment of HNSCC. However, the high IC50 values of other drugs in the high-risk group suggested that the low-risk group was more sensitive to these drugs (Figs. 6D–6H).

Validation of lncRNA expression features

After obtaining 11 lncRNAs from the prognostic model, we evaluated their prognostic value. K-M analysis in Fig. 7 indicated that ECE1-AS1, AL033527.2, AL162458.1, and AC090587.1 were positively correlated with OS, suggesting a good prognosis. Meanwhile, SNHG16 and AC012184.3 were negatively correlated with OS, indicating a poor prognosis. Overall, the results of K-M analysis were consistent with those of univariate Cox analysis, indicating that most lncRNAs had strong predictive ability when injected into specific prognostic models.

For experimental validation, we ultimately screened three prognostic lncRNAs with consistent expression. First, RT-qPCR was used to verify their expression in cell lines. Figures 8A–8C show that compared with those in the normal HBE line, AC090587.1 was downregulated in the HNSCC cell line TU686. AC012184.3 was upregulated in TU686, and SNHG16 expression was higher in TU686 and TU212 cell lines than in normal cell lines. Among the 12 pairs of clinically collected cancer and adjacent tissues, only SNHG16 showed differential expression (Figs. 8D–8F). These results suggest that SNHG16 can serve as a biomarker for the diagnosis and prognosis of HNSCC.