Liuqing Yang
Jinling Yu- Department of Prosthodontics, Hospital of Stomatology, Jilin University, Changchun, China
Background: Cuproptosis is a new type of cell death that induces protein toxic stress and eventually leads to cell death. Hence, regulating cuproptosis in tumor cells is a new therapeutic approach. However, studies on cuproptosis-related long noncoding RNA (lncRNA) in head and neck squamous cell carcinoma (HNSC) have not been found. This study aimed to explore the cuproptosis-related lncRNAs prognostic marker and their relationship to immune microenvironment in HNSC by using bioinformatics methods.
Methods: RNA sequencing, genomic mutations, and clinical data of TCGA_HNSC were downloaded from The Cancer Genome Atlas. HNSC patients were randomly assigned to either a training group or a validation cohort. The least absolute shrinkage and selection operator Cox regression and multivariate Cox regression models were used to determine the prognostic model in the training cohort, and its independent prognostic effect was further confirmed in the validation and entire cohorts.
Results: Based on previous literature, we collected 19 genes associated with cuproptosis. Afterward, 783 cuproptosis-related lncRNAs were obtained through coexpression. Cox model revealed and constructed eight cuproptosis-related lncRNAs prognostic marker (AL132800.1, AC090587.1, AC079160.1, AC011462.4, AL157888.1, GRHL3-AS1, SNHG16, and AC021148.2). Patients were divided into high- and low-risk groups based on the median risk score. The Kaplan–Meier survival curve revealed that the overall survival between the high- and low-risk groups was statistically significant. The receiver operating characteristic curve and principal component analysis demonstrated the accurate prognostic ability of the model. Univariate and multivariate Cox regression analysis showed that risk score was an independent prognostic factor. In addition, we used multivariate Cox regression to establish a nomogram of the predictive power of prognostic markers. The tumor mutation burden showed significant differences between the high- and low-risk groups. HNSC patients in the high-risk group responded better to immunotherapy than those in the low-risk group. We also found that risk scores were significantly associated with drug sensitivity in HNSC.
Conclusion: In summary, our study identified eight cuprotosis-related lncRNAs signature of HNSC as the prognostic predictor, which may be promising biomarkers for predicting the benefit of HNSC immunotherapy as well as drug sensitivity.
Introduction
Head and neck squamous cell carcinoma (HNSC) is a malignant tumor that affects different tissues and organs of the head and neck and poses a serious threat to human health (Siegel et al., 2018). HNSC is a complex disease. Typical risk factors are smoking, excessive alcohol consumption, and human papillomavirus (Jamal et al., 2016; Rooper et al., 2020). Although the current treatment methods for HNSC include surgery and chemoradiotherapy, the recurrence rate and metastasis risk of HNSC are still very high (Camisasca et al., 2011). In addition, the 5-year survival rate is less than 50% (Amit et al., 2013). An emerging approach to the treatment of HNSC is urgently needed. Hence, it is of great clinical significance to identify reliable biomarkers for predicting treatment response and prognosis and to develop effective treatment strategies for HNSC patients.
Cuproptosis is a new type of cell death that is different from the known cell death mechanism such as apoptosis, autophagy, and ferroptosis. When the known cell death mechanism is blocked, copper ions can still induce cell death. Cuproptosis occurs through the direct binding of copper ions to lipoacylated components of the tricarboxylic acid cycle in mitochondrial respiration, resulting in the aggregation of lipoacylated proteins. In addition, copper ions can reduce the protein level of the Fe–S cluster. They both induce protein toxic stress response and eventually lead to death (Tsvetkov et al., 2022). Based on this novel approach to cell death, we are developing new therapies for HNSC patients. Therefore, identifying the key regulators of cuproptosis is an important step toward further understanding.
Long noncoding RNA (lncRNA) are single-stranded RNAs with over 200 nucleotides in length, most of which do not have protein-coding capabilities (Gao et al., 2020). LncRNA regulates a variety of physiological and biochemical cellular processes by mediating chromosomal modification, transcriptional activation, and interference (Statello et al., 2021). Studies have shown that lncRNA is abnormally expressed and regulated in a variety of tumors (Castro-Oropeza et al., 2018). It has been reported that abnormal lncRNAs can be used as prognostic indicators of various cancers (Ai et al., 2020; Gai et al., 2020; Jiang et al., 2021). At present, there are few studies on cuproptosis-related lncRNAs and their association with the prognosis of HNSC patients. Therefore, this study aims to explore prognostic cuproptosis-related lncRNAs markers, to improve current strategies for diagnosis, treatment, follow-up, and prevention of HNSC.
In this study, we obtained HNSC RNA sequencing (RNA-seq) data downloaded from The Cancer Genome Atlas (TCGA) database and randomly assigned patients to a training and test datasets. We identified the cuproptosis-related lncRNAs prognostic marker (CRLPM) and developed an lncRNA signature prognostic model, which might represent potential therapeutic targets and provide valuable clinical utility for prognostic prediction of patients with HNSC. Last, we verified the predictive capacity of the model in order to provide a basis for the development of appropriate clinical strategies and revealed its potential to predict immunotherapy and drug sensitivity of HNSC.
Materials and Methods
Download and Processing of Transcriptomic Data, Mutation Data, and Clinical Information
The RNA-seq transcriptome profiling dataset comprised 44 normal tissues and 504 HNSC samples, which were downloaded from TCGA (https://portal.gdc.cancer.gov/) database on April 20, 2022. The tumor somatic mutation data and clinical information including survival time, survival status, age, gender, grade, stage, and tumor-node-metastasis classification were also obtained from TCGA. The annotations for lncRNAs were obtained from the GENCODE website (https://www.gencodegenes.org/). Furthermore, cuproptosis-related genes were obtained based on previous literature (Tsvetkov et al., 2022).
Generation and Assessment of the Cuproptosis-Related Long Noncoding RNA
The coexpression analysis between cuproptosis-related genes and lncRNAs was performed by the “limma” package in R to obtain the cuproptosis-related lncRNAs. Meeting the |Cor|>0.4 and p < 0.001 criteria indicated an association. According to the results of the coexpression analysis, we used R “ggplot2,” “ggalluvial,” and “dply” packages to generate the Sankey plot.
Prognostic Model Construction
The samples were randomly divided into training and validation groups through the R package “caret.” Univariate Cox proportional risk regression was performed for each cuproptosis-related lncRNAs with survival data using the survival R package. We performed the least absolute shrinkage and selection operator (Lasso)-penalized Cox regression by using the “glmnet” package in R software to avoid overfitting. The optimal and minimum criteria for the penalty (λ) using 10 times cross-validation were selected. Next, multiple stepwise COX regression analyses were performed to identify the CRLPM. Afterward, the formula of the risk scoring model was established as follows:
Validation of Risk Models
Patients in the training and validation groups were categorized into the high- and low-risk groups based on the median risk score and the corresponding coefficient of the training group. The Kaplan–Meier method was then conducted to display the prognostic performance of the risk score model in both the training and validation groups. In addition, the receiver operating characteristic (ROC) curve and the area under the curve (AUC) were used to evaluate the accuracy and diagnostic value of the CRLPM through the use of the survival ROC and time ROC packages in R. The principal component analysis (PCA) was also conducted to validate risk models, and the results were visualized using “scatterplot3D” packages in R software. The progression-free survival (PFS) was performed through “survival” and “survminer” packages in R. We used the C-index to predict the accuracy of risk models by using the R package “rms,” “dplyr,” “survival,” and “pec.” The validation and entire cohorts were performed to validate this model.
Establish and Evaluate a Nomogram
We used univariate and multivariate Cox regression to investigate the independent prognostic role of the risk model. Based on the results of univariate and multivariate COX regression, we developed a nomogram by employing the R package “rms,” “regplot,” and “survival.” The accuracy of the nomogram was evaluated using a calibration curve.
Exploration of the Relationship Between the Prognostic Risk Score and Clinical Stage
To verify whether the model is suitable for patients with different clinical stages, we explore the relationships between risk score and clinical stage to reveal their possible roles in HNSC using univariate and multivariate Cox regression analyses.
Pathway Enrichment Analysis and Gene set Enrichment Analysis
The differentially expressed genes (DEGs) between the high- and low-risk groups were identified using the R package “limma,” with the limited condition set to log2 |fold change| >1 and false discovery rate < 0.05. Based on the R package “clusterProfiler,” “org.Hs.eg.db,” and “enrichplot,” we explored the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database pathways to clarify the molecular functions and key signaling pathways.
Estimation of Intratumoural Immune Cell Infiltration and Immunotherapy
To investigate the relationship between the CRLPM risk score and immune cell infiltration, the single-sample gene set enrichment analysis (ssGSEA) algorithms function package in the R software genome variation analysis package was used to evaluate the infiltration and function of tumor-infiltrating immune cells. The related heat map was utilized to drawn. Next, based on the simulation of the tumor immune escape mechanism, the tumor immune dysfunction and exclusion (TIDE) algorithm was applied to predict the response to immunotherapy (http://tide.dfci.harvard.edu) (Jiang et al., 2018). Therefore, we observed the effect of immunotherapy in the high- and low-risk groups based on the TIDE algorithm.
Evaluation of Drug Sensitivity
IC50 represented the semiinhibitory concentration of the measured antagonist. To evaluate CRLPM in the clinic for HNSC treatment, we calculated the IC50 of the chemotherapeutic drugs through the “pRRophetic” R package and its dependencies including “car, ridge preprocessCore, genefilter and sva.” A total of 138 drugs were included such as midostaurin, temsirolimus, tipifarnib, and imatinib. The Wilcoxon sign rank test was used to compare IC50 differences between common antineoplastic agents in the high- and low-risk groups. The boxplot was presented using the R package “ggplot2.”
Calculation of Tumor Mutation Burden Scores
Tumor mutational burden (TMB) reflects the number of mutations in cancer mutation. The mutation data of HNSC samples downloaded from TCGA were analyzed using the R package “maftools.” The waterfall diagram showed the relationship between risk scores and TMB in HNSC patients.
Statistical Analysis
All statistical analyses were processed by the R programming language (Version 4.0.3) on R studio. RNA-seq transcriptome data and somatic mutation data downloaded from TCGA were combined using the “limma” package in R. Pearson correlation test was used to analyze the correlations between cuproptosis-related genes and cuproptosis-related lncRNAs. CRLPM were screened for differential genes using the “limma” R package. Cox regression and survival analysis were performed through “survival” and “survminer” packages in R. Cox proportional risk regression model was used to calculate the hazard ratios of univariate and multivariate analyses. The GO terms and KEGG pathways were analyzed by using “clusterProfiler” in the R package. The “Pheatmap” R package was used to draw heat maps in cluster analysis. We applied the Wilcoxon rank-sum test to compare the difference between two groups of quantitative data. The overall survival (OS) time of the different groups was evaluated using the Kaplan–Meier analysis with a log-rank test. The chi-square test was used to compare categorical data between different groups. A p value of <0.05 was considered statistically significant.
Results
Data Processing
We removed the genes encoding proteins and identified 16876 lncRNAs in the TCGA_HNSC dataset through the “GENCODE” database. In total, we collected 19 cuproptosis-related genes. Based on Pearson analysis, 783 cuprotosis-related lncRNAs were obtained. The Sankey plot showed the association between cuproptosis-related genes and cuproptosis-related lncRNAs (Figure 1A). Thenceforth, univariate COX regression analysis was applied to explore cuproptosis-related lncRNAs (p < 0.05). The 501 patients were divided into the training group (n = 251) and the validation group (n = 250), and the clinical information on HNSC was presented in Table 1. The results showed that there was no difference between the training and validation groups in all clinical traits.
FIGURE 1. Sankey diagram and heat map. (A) Sankey diagram of coexpression between 19 cuproptosis-related genes and 783 cuproptosis-related long noncoding RNA (lncRNAs). (B) correlation 19 cuproptosis-related genes and 8 prognostic cuproptosis-related lncRNAs. *p < 0.05, **p < 0.01, and ***p < 0.001.
TABLE 1. Characteristic of head and neck squamous cell carcinoma patients.
Construction and Validation of the Cuproptosis-Related Long Noncoding RNAs Prognostic Marker
The univariate COX analysis of 21 cuproptosis-related lncRNAs was shown in Figure 2A. We further screened 17 lncRNAs using Lasso–Cox regression. We identified trajectory changes in regression coefficients of lncRNAs and cross-validation results of model construction (Figure 2B, C). Afterward, multiple stepwise Cox regression analysis was performed, and we screened out eight CRLPM with survival to establish the risk score models. In conclusion, we generated a total of eight CRLPMs to participate in the construction of a prognostic model to predict the OS of patients with HNSC. Afterward, we obtained a heat map of the correlation between cuproptosis-related genes and CRLPM (Figure 1B). Risk score = (0.351771323*ExpressionAL132800.1) + (−0.378321346*ExpressionAC090587.1) + (0.404882025*ExpressionAC079160.1) + (−0.314303555*ExpressionAC011462.4) + (0.716547372* ExpressionAL157888.1) + (−0.593212656*ExpressionGRHL3-AS1) + (0.390726744*ExpressionSNHG16) + (−0.892732753*ExpressionAC021148.2). Based on the median risk score, we divided the patients in the training group into the high- and low-risk groups for survival analysis. The KM method was used to analyze the OS of patients in the two groups, and the results showed that the OS of patients in the high-risk group was significantly poorer than that in the low-risk group (p < 0.05; Figure 2D). The distribution of risk scores and the survival status of patients were shown in Figure 2E. The expression level of eight cuproptosis-related lncRNAs involved in the high- and low-risk groups was shown in a heatmap (Figure 2F). It can be observed that with the increase in risk score, the survival time was shortened and the number of deaths increased. We also observed statistically significant differences in OS between the high- and low-risk groups in the validation and entire cohorts (p < 0.05; Figure 2G–I, Figure 3A–C). The PCA revealed a high degree of differentiation between the high- and low-risk groups. Based on the risk model of cuproptosis-related lncRNAs, we intuitively observed that HNSC patients were effectively divided into two clusters (Figure 3D–G).
FIGURE 2. Construction of the prognostic cuproptosis-related long noncoding RNA (lncRNAs) risk model in head and neck squamous cell carcinoma (HNSC). (A) univariate Cox regression analysis for identifying the prognostic cuproptosis-related lncRNAs. (B–C) Lasso–Cox regression analysis was performed to construct prognostic prediction models. (D) Kaplan–Meier curves for survival analysis in the high- and low-risk groups. (E) risk score distribution and survival status in patients with HNSC. (F) heatmap of the prognostic markers and overall survival. (G) Kaplan–Meier curves for survival analysis in the validation cohort. (H) risk score distribution and survival status in the validation cohort. (I) heatmap of the prognostic markers and overall survival in the validation cohort.
FIGURE 3. Validation of the risk model in the entire cohort and principal component analysis. (A) Kaplan–Meier curves for survival analysis in the entire cohort. (B) risk score distribution and survival status in the entire cohort. (C) heatmap of the prognostic markers and overall survival in the entire cohort. PCA between the high- and low-risk groups based on the (D) all genes, (E) cuproptosis-related genes, (F) cuproptosis-related long noncoding RNA (lncRNAs), and (G) cuproptosis-related lncRNAs prognostic marker.
Independence of the Cuproptosis-related lncRNAs Prognostic Marker in Predicting Overall Survival
Univariate and multivariate Cox regression analyses were performed to assess the predictive value of the prognostic model. In univariate Cox analysis, there were statistically significant differences among age, stage, and risk score (Figure 4A). In multivariate Cox regression analysis, they remained prognostic value for OS (Figure 4B). The PFS indicated significant differences in progression-free survival between the high- and low-risk groups (p < 0.05, Figure 4C). The ROC curves demonstrated the accuracy and diagnostic value of the cuproptosis-related lncRNAs for OS, and the AUC reached 0.690 at 1 year, 0.701 at 2 years, and 0.668 at 3 years (Figure 4D). Both C-index and ROC curve indicated the predictive accuracy of the prognostic model was superior to other clinical including age, gender, grade, and stage (Figures 4E, F). A nomogram plot is a predictive tool for quantitative analysis of clinical outcomes in patients with HNSC. Thus, we initiated a prognostic nomogram based on the risk score and other clinical characteristics (Figure 5A). The calibration plots showed good conformity with the prediction of this nomogram (Figure 5B).
FIGURE 4. Independent prognostic analysis of head and neck squamous cell carcinoma (HNSC) overall survival (OS). (A) univariate Cox analysis. Age, stage, and risk score were statistically significant. (B) multivariate Cox analysis. Age, stage, and risk score were statistically significant. (C) Kaplan–Meier curves of progression-free survival (PFS). (D) TimeROC curve predicted 1, 3, and 5 years of OS for HNSC patients. (E) ROC demonstrated the predictive accuracy of the risk model was superior to other clinical parameters. (F) C-index showed the predictive accuracy of the risk model was superior to other clinical parameters.
FIGURE 5. Construction and evaluation of a nomogram based on CRLPM. (A) nomogram used to predict prognosis was constructed based on CRLPM. (B) calibration curves are used to predict 1-, 3-, and 5-year overall survival. (C) Kaplan–Meier curves of patients with stage I-II. (D) Kaplan–Meier curves of patients with stage III–IV.
Relationship Between the Marker and the Clinical Features in Head and Neck Squamous Cell Carcinoma
Next, to investigate the clinical utility of the CRLPM, we explored the relationship of the CRLPM with clinical features. The results indicated that there were significant differences between the distribution of risk scores and clinical stages. In specific, stages I–II and III–IV were statistically significant (p < 0.05, Figure 5C, D).
Pathway Enrichment Analysis and Gene Set Enrichment Analysis
To explore the biological functions and pathway analysis of DEGs between the high- and low-risk groups, we further performed the GO and KEGG enrichment analyses. A total of 359 DEGs were identified. In the biological process category, the genes were primarily concentrated in response to humoral immune response, immune response-activating cell surface receptor signaling pathway, and adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains. In the cellular component category, it was mainly enriched in the immunoglobulin complex, external side of plasma membrane and apical part of cell. In the molecule function category, it was antigen binding, immunoglobulin receptor binding, and receptor ligand activity (Figure 6A, C, E). Genes in the KEGG category were enriched in the IL-17 signaling pathway, hematopoietic cell lineage, and amoebiasis (Figure 6B, D, F).
FIGURE 6. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. (A) barplot of the top 10 GO enrichment terms. (B) barplot of the top 30 KEGG enrichment terms. (C) bubble chart of the top 10 GO enrichment terms. (D) bubble chart of the top 30 KEGG enrichment terms. (E) circle diagram of GO enrichment analysis. (F) circle diagram of KEGG enrichment analysis. Biological process, cellular component, and molecular function.
Estimation of Intratumoural Immune Cell Infiltration and Immunotherapy
Figure 7A shows the heatmap of immune response based on the ssGSEA algorithm. Based on ssGSEA of TCGA-HNSC data, correlation analysis between immune cell populations and related functions revealed that T cell functions including regulation of inflammation, HLA, checkpoint (inhibition), and costimulation and coinhibition were significantly different between the high- and low-risk groups. These results indicated that GRLPM was associated with immune cell infiltration in HNSC. Based on the TIDE algorithm, we predicted the effect of patients receiving immunotherapy. Figure 7B shows significant differences in TIDE scores between the high- and low-risk groups and the lower TIDE scores in the high-risk group. This further proves that patients in the high-risk group have a low potential for immune escape and may receive better results from immunotherapy.
FIGURE 7. Immunological landscape in head and neck squamous cell carcinoma (HNSC) patients and relationship between tumor mutation burden (TMB) and risk score. (A) heatmap of the tumor-infiltrating lymphocytes based on single-sample gene set enrichment analysis algorithms among the high- and low-risk groups in HNSC. *p < 0.05, **p < 0.01, and ***p < 0.001. (B) comparison of TIDE prediction score between the high- and low-risk groups. (C) Analysis of TMB differences between the high- and low-risk groups in HNSC. (D) survival analysis curves of the high- and low-TMB groups. (E) TMB risk combined with survival curve in HNSC. (F) waterfall plot of top 15 mutant genes in the high-risk group in HNSC. (G) waterfall plot of top 15 mutant genes in the low-risk group in HNSC.
Tumor Mutational Burden of the Cuproptosis-related lncRNAs Prognostic Marker in Head and Neck Squamous Cell Carcinoma Samples
We collected data on somatic mutation in HNSC and calculated corresponding TMB scores in order to investigate the potential role of tumor mutation load in HNSC. As shown in Figure 7C, the high-risk group had a higher mutation load than the low-risk group in HNSC. We divided patients into “Hight-TMB” and “Low-TMB” by median cutoff points and performed survival analyses. The results showed that the high-risk group had a lower survival rate than the low-risk group in HNSC (Figure 7D). A combined survival analysis of tumor mutation load and risk scores can obtain the combined survival curve. It revealed that the TMB and risk scores had significant effects on the OS of HNSC patients (Figure 7E). The mutation landscapes in the CRLPM high- and low-risk groups were compared. Waterfall plots visualized the 15 genes with the highest mutation frequency in the high- and low-risk groups. The results showed that more mutation events occurred in the high-risk group (Figures 7F, G). TP53 was the gene with the highest mutation frequency.
Drug Sensitivity
To explore the possible use of CRLPM in the individualized treatment of HNSC, we investigated the relationship between risk scores and IC50 of drugs in HNSC treatment. To this end, we compared the sensitivity of 30 common anticancer drugs between the high- and low-risk groups. As shown in Figure 8, the sensitivity of 12 of the 30 anticancer drugs was significantly different in the high- and low-risk groups (p <0.05). Meanwhile, 11 of the drugs had lower IC50s in the high-risk group, further proving that the high-risk group was more sensitive to drug treatment. This means that these drugs have a potential role in the treatment of HNSC in the future. However, the IC50s of temsirolimus were higher in the high-risk group, which suggests that the low-risk group had a high sensitivity to this drug. Results showed that in HNSC patients, except for temsirolimus, the risk score was inversely associated with IC50 (Supplementary Figure S1).
FIGURE 8. Drug sensitivity (IC50) correlated with high- and low-risk patients in head and neck squamous cell carcinoma.
Discussion
Despite advances in surgery and chemotherapy in recent years, the prognosis for patients with advanced and metastatic HNSC remains poor (Miyauchi et al., 2019; Johnson et al., 2020). Cuproptosis overcomes the resistance of malignant cells to chemotherapy and helps remove defective cells. Therefore, cuproptosis may be an effective way to treat many types of cancer in the future. In addition, lncRNA affects the development and treatment of cancer through biological means (Gao et al., 2021; Tan et al., 2021). LncRNAs have been found to play an important role in the prognosis of HNSC and may be a potential effective molecular target for the treatment of HNSC (Ban et al., 2020; Tang et al., 2021; Zhu et al., 2021). However, the regulatory mechanisms of cuproptosis remain largely unknown, especially in the field of lncRNA. Therefore, we should focus on the potential interaction between lncRNA and cuproptosis to uncover potential prognostic markers.
In the current study, we identified 21 prognostic cuproptosis-related lncRNAs, eight of which were selected for prognostic construction to predict OS in patients with HNSC. First, 19 cuproptosis-related genes and 783 cuproptosis-related lncRNAs were obtained. Then, we used Lasso regression and COX regression to identify the prognostic cuproptosis-related lncRNAs. In addition, we further explored the CRLPM and common clinical variables, upstream regulatory mechanisms, immune cell infiltration and immunotherapy, and drug sensitivity of HNSC.
Based on Cox, Lasso, and multivariate Cox regression analyses, we identified eight lncRNAs associated with prognosis, including AL132800.1, AC090587.1, AC079160.1, AC011462.4, AL157888.1, GRHL3-AS1, SNHG16, and AC021148.2. Guo et al. (2021) found that AC079160.1 is a prognostic biomarker for gastric cancer, and AC079160.1 was found to be overexpressed. In general, high expression of AC079160.1 was associated with better survival in gastric cancer. AC011462.4 has been reported to play an important oncogenic role in tumors. Li et al. (2021) revealed that the high expression level of AC011462.4 was associated with longer OS. The expression level of AC011462.4 increased with the increase of risk score in colon cancer. GRHL3-AS1 expression was upregulated in patients with primary HNSC, and its expression was associated with the survival of patients with HNSC (Feng et al., 2021). Small nucleolar RNA host gene 16 (SNHG16) is considered to be a cancer-associated lncRNA that promotes tumor development primarily by acting as a competing endogenous RNA (ceRNA) (Grüll and Massé, 2019). There is evidence that SNHG16 acts as a ceRNA in various cancer by sponging corresponding miRNA to regulate mRNA (Zhao et al., 2018). In addition to the ceRNA mechanism, SNHG16 plays a role in promoting cancer through other mechanisms. SNHG16 can promote the proliferation and inhibit apoptosis of bladder cancer by inhibiting the expression of P21 (Cao et al., 2018). SNHG16 was found to be significantly upregulated in a variety of tumor tissues and cell lines, such as hepatocellular carcinoma, lung cancer, colorectal cancer, glioma, and other tumor types (Christensen et al., 2016; Yang et al., 2018; Chen et al., 2019; Su et al., 2019). In addition, high SNHG16 expression was associated with a poor prognosis. One study showed that in oral squamous cell carcinoma, the expression of SNHG16 was upregulated by c-Myc (Li et al., 2019). However, there are few reports on AL132800.1, AC090587.1, AL157888.1, and AC021148.2. Thus, it is necessary to further determine their mechanisms during cuproptosis through experiments in our future studies.
Afterward, we verify the accuracy of the risk model. The Kaplan-Meier method showed that the OS of the high-risk group was lower than that of the low-risk group. The ROC curves showed that CRLPM had high accuracy in predicting 1-, 3-, and 5-year survival, and AUC were all greater than 0.65. PCA intuitively showed differences between the high- and low-risk groups. PFS, C-index combining CRLPM with clinical information, a new nomogram was created to predict prognosis, lymph node metastasis, and distant metastasis in HNSC patients.
In addition, functional enrichment analyses revealed the potential biological mechanism of the involved CRLPM. We explored the key signaling pathways of eight cuproptosis-related lncRNAs. GO and KEGG analysis indicated that this differentially expressed CRLPM was mainly enriched in the IL-17 signaling pathway, hematopoietic cell lineage, and amoebiasis.
Our results found that the TMB was statistically higher in the high-risk group than in the low-risk group, suggesting that patients at high risk of HNSC have a better response to immunotherapy. Among the first 15 mutated genes, TP53 was mutated more frequently in HNSC patients. Moreover, the drug sensitivity of these CRPLM was analyzed to guide clinical treatment. There were significant differences in IC50 between high-risk and low-risk patients for all 12 drugs.
Cuproptosis is a new form of cell death that may play an important role in future cancer treatments. On the other hand, some lncRNAs influence cancer progression and treatment in a variety of biological ways. However, there is still a lot of unexplored territory between cuproptosis and lncRNA. Overall, this study provides new insights into the tumorigenesis and progression of HNSC from the perspective of cuproptosis. Biomarkers of cuproptosis that can be used for the prognosis of HNSC were explored, which could inform the treatment of the disease. We inevitably used only the TCGA validation and entire cohorts, and additional patients could improve the reliability of the model; thus, further validation is needed through preclinical studies. In the meantime, our study needs further validation in vivo and in vitro soon.
Conclusion
To sum up, we identified a risk model based on seven cuproptosis-related lncRNAs that accurately predicted the prognosis of HNSC. Hence, our study provided a new therapeutic strategy for individualized therapy and immunotherapy response in HNSC patients. These eight cuproptosis-related lncRNAs may be therapeutic targets for HNSC.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author Contributions
LY analyzed the data regarding the oral tongue squamous cell carcinomas and was a major contributor to writing the manuscript. JY modified the grammar of the article. HH, YG, LT, ZL, and JY collected data. All authors read and approved the final manuscript.
Funding
This research were funded by Jilin Provincial Department of Finance Science and Technology Project jcsz grant number 2021893-9 and Jilin Provincial Development and Reform Commission 2021 Industrial Technology Research and Development Project grant number 2021C043-2.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as potential conflicts of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.947551/full#supplementary-material
Supplementary Figure S1 | Linear correlation between IC50 and risk score.
References
Ai, W., Li, F., Yu, H. H., Liang, Z. H., and Zhao, H. P. (2020). Up‐regulation of Long Noncoding RNA LINC00858 Is Associated with Poor Prognosis in Gastric Cancer. J. Gene Med. 22 (7), e3179. doi:10.1002/jgm.3179
Amit, M., Yen, T.-C., Liao, C.-T., Chaturvedi, P., Agarwal, J. P., Kowalski, L. P., et al. (2013). Improvement in Survival of Patients with Oral Cavity Squamous Cell Carcinoma: An International Collaborative Study. Cancer 119 (24), 4242–4248. doi:10.1002/cncr.28357
Ban, Y., Tan, P., Cai, J., Li, J., Hu, M., Zhou, Y., et al. (2020). LNCAROD Is Stabilized by m6A Methylation and Promotes Cancer Progression via Forming a Ternary Complex with HSPA1A and YBX1 in Head and Neck Squamous Cell Carcinoma. Mol. Oncol. 14 (6), 1282–1296. doi:10.1002/1878-0261.12676
Camisasca, D. R., Silami, M. A. N. C., Honorato, J., Dias, F. L., Faria, P. A. S. d., and Lourenço, S. d. Q. C. (2011). Oral Squamous Cell Carcinoma: Clinicopathological Features in Patients with and without Recurrence. Orl 73 (3), 170–176. doi:10.1159/000328340
Cao, X., Xu, J., and Yue, D. (2018). LncRNA-SNHG16 Predicts Poor Prognosis and Promotes Tumor Proliferation through Epigenetically Silencing P21 in Bladder Cancer. Cancer Gene Ther. 25 (1-2), 10–17. doi:10.1038/s41417-017-0006-x
Castro-Oropeza, R., Melendez-Zajgla, J., Maldonado, V., and Vazquez-Santillan, K. (2018). The Emerging Role of lncRNAs in the Regulation of Cancer Stem Cells. Cell Oncol. 41 (6), 585–603. doi:10.1007/s13402-018-0406-4
Chen, H., Li, M., and Huang, P. (2019). LncRNA SNHG16 Promotes Hepatocellular Carcinoma Proliferation, Migration and Invasion by Regulating miR-186 Expression. J. Cancer 10 (15), 3571–3581. Published 2019 Jun 9. doi:10.7150/jca.28428
Christensen, L. L., True, K., Hamilton, M. P., Nielsen, M. M., Damas, N. D., Damgaard, C. K., et al. (2016). SNHG16 Is Regulated by the Wnt Pathway in Colorectal Cancer and Affects Genes Involved in Lipid Metabolism. Mol. Oncol. 10 (8), 1266–1282. doi:10.1016/j.molonc.2016.06.003
Feng, Z.-Y., Gao, H.-Y., and Feng, T.-D. (2021). Immune Infiltrates of m6A RNA Methylation-Related lncRNAs and Identification of PD-L1 in Patients with Primary Head and Neck Squamous Cell Carcinoma. Front. Cell Dev. Biol. 9, 672248. Published 2021 Jun 4. doi:10.3389/fcell.2021.672248
Gai, C., Liu, C., Wu, X., Yu, M., Zheng, J., Zhang, W., et al. (2020). MT1DP Loaded by Folate-Modified Liposomes Sensitizes Erastin-Induced Ferroptosis via Regulating miR-365a-3p/NRF2 axis in Non-small Cell Lung Cancer Cells. Cell Death Dis. 11 (9), 751. Published 2020 Sep 14. doi:10.1038/s41419-020-02939-3
Gao, N., Li, Y., Li, J., Gao, Z., Yang, Z., Li, Y., et al. (2020). Long Non-coding RNAs: The Regulatory Mechanisms, Research Strategies, and Future Directions in Cancers. Front. Oncol. 10, 598817. Published 2020 Dec 18. doi:10.3389/fonc.2020.598817
Gao, Y., Shang, S., Guo, S., Li, X., Zhou, H., Liu, H., et al. (2021). Lnc2Cancer 3.0: an Updated Resource for Experimentally Supported lncRNA/circRNA Cancer Associations and Web Tools Based on RNA-Seq and scRNA-Seq Data. Nucleic Acids Res. 49 (D1), D1251–D1258. doi:10.1093/nar/gkaa1006
Grüll, M. P., and Massé, E. (2019). Mimicry, Deception and Competition: The Life of Competing Endogenous RNAs. WIREs RNA 10 (3), e1525. doi:10.1002/wrna.1525
Guo, Z., Liang, E., Zhang, T., Xu, M., Jiang, X., and Zhi, F. (2021). Identification and Validation of a Potent Multi-lncRNA Molecular Model for Predicting Gastric Cancer Prognosis. Front. Genet. 12, 607748. Published 2021 Dec 20. doi:10.3389/fgene.2021.607748
Jamal, A., King, B. A., Neff, L. J., Whitmill, J., Babb, S. D., and Graffunder, C. M. (2016). Current Cigarette Smoking Among Adults - United States, 2005-2015. MMWR Morb. Mortal. Wkly. Rep. 65 (44), 1205–1211. Published 2016 Nov 11. doi:10.15585/mmwr.mm6544a2
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1
Jiang, W., Song, Y., Zhong, Z., Gao, J., and Meng, X. (2021). Ferroptosis-Related Long Non-coding RNA Signature Contributes to the Prediction of Prognosis Outcomes in Head and Neck Squamous Cell Carcinomas. Front. Genet. 12, 785839. Published 2021 Dec 17. doi:10.3389/fgene.2021.785839
Johnson, D. E., Burtness, B., Leemans, C. R., Lui, V. W. Y., Bauman, J. E., and Grandis, J. R. (2020). Head and Neck Squamous Cell Carcinoma. Nat. Rev. Dis. Prim. 6 (1), 92. Published 2020 Nov 26. doi:10.1038/s41572-020-00224-3
Li, S. S., Zhang, J., and Chen, J. (2019). c-Myc Induced Upregulation of Long Non-coding RNA SNHG16 Enhances Progression and Carcinogenesis in Oral Squamous Cell Carcinoma. Cancer Gene Ther. 26 (11-12), 400–410. doi:10.1038/s41417-018-0072-8
Li, X., Yu, H., Wei, Z., Gou, X., Liang, S., and Liu, F. (2021). A Novel Prognostic Model Based on Autophagy-Related Long Non-coding RNAs for Clear Cell Renal Cell Carcinoma. Front. Oncol. 11, 711736. Published 2021 Aug 3. doi:10.3389/fonc.2021.711736
Miyauchi, S., Kim, S. S., Pang, J., Gold, K. A., Gutkind, J. S., Califano, J. A., et al. (2019). Immune Modulation of Head and Neck Squamous Cell Carcinoma and the Tumor Microenvironment by Conventional Therapeutics. Clin. Cancer Res. 25 (14), 4211–4223. doi:10.1158/1078-0432.CCR-18-0871
Rooper, L. M., Windon, M. J., Hernandez, T., Miles, B., Ha, P. K., Ryan, W. R., et al. (2020). HPV-Positive Squamous Cell Carcinoma of the Larynx, Oral Cavity, and Hypopharynx. Am. J. Surg. Pathol. 44 (5), 691–702. doi:10.1097/PAS.0000000000001433
Siegel, R. L., Miller, K. D., and Jemal, A. (2018). Cancer Statistics, 2018. CA A Cancer J. Clin. 68 (1), 7–30. doi:10.3322/caac.21442
Statello, L., Guo, C.-J., Chen, L.-L., and Huarte, M. (2021). Gene Regulation by Long Non-coding RNAs and its Biological Functions. Nat. Rev. Mol. Cell Biol. 22 (2), 96–118. doi:10.1038/s41580-020-00315-9
Su, P., Mu, S., and Wang, Z. (2019). Long Noncoding RNA SNHG16 Promotes Osteosarcoma Cells Migration and Invasion via Sponging miRNA-340. DNA Cell Biol. 38 (2), 170–175. doi:10.1089/dna.2018.4424
Tan, Y. T., Lin, J. F., Li, T., Li, J. J., Xu, R. H., and Ju, H. Q. (2021). LncRNA‐mediated Posttranslational Modifications and Reprogramming of Energy Metabolism in Cancer. Cancer Commun. 41 (2), 109–120. doi:10.1002/cac2.12108
Tang, Y., Li, C., Zhang, Y.-J., and Wu, Z.-H. (2021). Ferroptosis-Related Long Non-coding RNA Signature Predicts the Prognosis of Head and Neck Squamous Cell Carcinoma. Int. J. Biol. Sci. 17 (3), 702–711. Published 2021 Jan 31. doi:10.7150/ijbs.55552
Tsvetkov, P., Coy, S., Petrova, B., Dreishpoon, M., Verma, A., Abdusamad, M., et al. (2022). Copper Induces Cell Death by Targeting Lipoylated TCA Cycle Proteins. Science 375 (6586), 1254–1261. doi:10.1126/science.abf0529
Yang, B. Y., Meng, Q., Sun, Y., Gao, L., and Yang, J. X. (2018). Long Non-coding RNA SNHG16 Contributes to Glioma Malignancy by Competitively Binding miR-20a-5p with E2F1. J. Biol. Regul. Homeost. Agents 32 (2), 251–261.
Zhao, W., Fu, H., Zhang, S., Sun, S., and Liu, Y. (2018). LncRNA SNHG16 Drives Proliferation, Migration, and Invasion of Hemangioma Endothelial Cell through Modulation of miR‐520d‐3p/STAT3 axis. Cancer Med. 7 (7), 3311–3320. doi:10.1002/cam4.1562
Keywords: head and neck squamous cell carcinoma, immune microenvironment, tumor mutational burden, cuproptosis, lncRNA, prognostic marker, drug sensitivity
Citation: Yang L, Yu J, Tao L, Huang H, Gao Y, Yao J and Liu Z (2022) Cuproptosis-Related lncRNAs are Biomarkers of Prognosis and Immune Microenvironment in Head and Neck Squamous Cell Carcinoma. Front. Genet. 13:947551. doi: 10.3389/fgene.2022.947551
Received: 18 May 2022; Accepted: 20 June 2022;
Published: 22 July 2022.
Edited by:
Hifzur R. Siddique, Aligarh Muslim University, IndiaReviewed by:
Samapika Routray, All India Institute of Medical Sciences Bhubaneswar, IndiaFred Bunz, SOURCE, Johns Hopkins University, United States
Copyright © 2022 Yang, Yu, Tao, Huang, Gao, Yao and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhihui Liu, liu_zh@jlu.edu.cn