Identification and verification of aging-related lncRNAs for prognosis prediction and immune microenvironment in patients with head and neck squamous carcinoma

Aging is highly associated with tumor formation and progression. However, little research has explored the association of aging-related lncRNAs (ARLs) with the prognosis and tumor immune microenvironment (TIME) of head and neck squamous cell carcinoma (HNSCC). RNA sequences and clinicopathological data of HNSCC patients and normal subjects were downloaded from The Cancer Genome Atlas. In the training group, we used Pearson correlation, univariate Cox regression, least absolute shrinkage/selection operator regression analyses, and multivariate Cox regression to build a prognostic model. In the test group, we evaluated the model. Multivariate Cox regression was done to screen out independent prognostic factors, with which we constructed a nomogram. Afterward, we demonstrated the predictive value of the risk scores based on the model and the nomogram using time-dependent receiver operating characteristics. Gene set enrichment analysis, immune correlation analysis, and half-maximal inhibitory concentration were also performed to reveal the different landscapes of TIME between risk groups and to predict immuno- and chemo-therapeutic responses. The most important LINC00861 in the model was examined in HNE1, CNE1, and CNE2 nasopharyngeal carcinoma cell lines and transfected into the cell lines CNE1 and CNE2 using the LINC00861-pcDNA3.1 construct plasmid. In addition, CCK-8, Edu, and SA-β-gal staining assays were conducted to test the biofunction of LINC00861 in the CNE1 and CNE2 cells. The signature based on nine ARLs has a good predictive value in survival time, immune infiltration, immune checkpoint expression, and sensitivity to multiple drugs. LINC00861 expression in CNE2 was significantly lower than in the HNE1 and CNE1 cells, and LINC00861 overexpression significantly inhibited the proliferation and increased the senescence of nasopharyngeal carcinoma cell lines. This work built and verified a new prognostic model for HNSCC based on ARLs and mapped the immune landscape in HNSCC. LINC00861 is a protective factor for the development of HNSCC.


Introduction
Head and neck squamous cell carcinoma (HNSCC) is a type of malignancy originating from the mucosal epithelium of the lip, oral cavity, larynx, or naso-, oro-, hypo-pharynx. HNSCC has a mortality rate of nearly 50% in the first five years after diagnosis, with over 850,000 new cases and more than 400,000 deaths globally in 2020 [1][2][3]. Though diverse treatment modalities such as surgery, targeted therapy, chemo-, radio-and immuno-therapy are used in HNSCC patients, their clinical benefits are far from satisfactory [4].
Aging is an irreversible growth arrest condition closely associated with the development of many chronic diseases and cancers [5]. Cellular senescence can be induced by stimuli such as telomere shortening due to extensive replication, oxidative stress, DNA damage, and oncogene overexpression, which are regulated by senescence-related genes [6]. After senescent cells lose proliferative capacity, hindering tumor progression, they remain viable and metabolically active, secreting various cytokines that promote the growth, migration, and invasion of tumor cells [7,8]. The potential prognostic value of aging-related genes in HNSCC and the correlations with inflammation and tumor immunity have recently been explored and confirmed [9].
Long non-coding RNAs (lncRNAs) are a class of nonprotein-coding RNAs that are more than 200 nucleotides in length and have specific functions (e.g., splicing, transcriptional and post-transcriptional modulation of mRNA), which are largely linked to tumor development, metastasis, and tumor immunity [10]. The action mechanisms of lncRNAs in aging have not yet been elucidated, and the prognostic impact of lncRNAs associated with aging-related genes in HNSCC remains unknown.
Here, we use bioinformatics to probe into the mechanisms of aging-related lncRNAs (ARLs) in HNSCC, identify some new molecular biomarkers linked with prognosis, and build an effective prognostic model to forecast the survival of HNSCC patients and function as their new treatment target.

Data collection
We downloaded the RNA sequencing data and clinical data of HNSCC and normal tissues from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). Aging-related genes were acquired from the Human Ageing Genomic Resources (HAGR ; Table S1). Expression data were standardized to fragment per kilobase million. Based on Strawberry Perl, we integrated and processed the data, distinguished between lncRNA and mRNA, extracted complete clinical information, and removed data missing survival status and unknown survival time or <30 days.

Identification of ARLs in HNSCC
We performed co-expression analysis of lncRNAs and AGs in tumor samples (correlation coefficient > 0.5 and p < 0.001) with packages Strawberry Perl and limma R and mapped a network with igraph R. Differential expression of ARLs in cancer versus normal tissue samples was analyzed (Log2 fold change (FC) > 1, false discovery rate (FDR) < 0.05, p < 0.05) [11], and a heatmap was plotted with the R package pheatmap.

Construction of the Prognostic Signature
Based on the clinical data of HNSCC, lncRNAs with significant effects on OS were screened from ARLs (p < 0.05) using univariate Cox analysis. Samples were randomly and evenly separated into a training risk group and a test risk group on the R package caret. In the training group, we conducted the least absolute shrinkage and selection operator (LASSO) Cox regression on the R package glmnet and multivariate Cox analysis to recognize the ARLs for risk models and computed the risk score as follows: Risk score ¼ Expression of the 1 st lncRNA Ã coefficient þ Expression of the 2 nd lncRNA Ã coefficient þ Expression of the nth lncRNA Ã coefficient where the coefficient denoted the regression coefficient of the corresponding lncRNA [12]. The overall group was classified by the median risk score into low-and high-risk groups.

Verification of risk signature and establishment of nomogram
We performed the Kaplan-Meier survival curve and log-rank test on the R package survminer to clarify whether the lowrisk group has longer OS than the high-risk group. The 1, 3and 5-year survival receiver's operating characteristic (ROC) curves and the areas under the curves (AUCs) were conducted on the R package timeROC to assess the prediction efficacy in comparison with other clinical data [13]. Based on the R package rms, a nomogram for forecasting survival was built on basis of the risk score and clinical features, which were detected by both uni-/multivariate Cox regression. We assessed the accuracy of the nomogram via calibration curve and ROC analysis.
Gene set enrichment analyses (GSEA) and assessment of the immune landscape We performed GSEA (version 4.2.3) using the curated gene set (kegg.v7.5.1 symbols.gmt), screened at p < 0.05 and FDR < 0.25, to determine the different biofunctions and pathways between the two groups [14]. Immune cell infiltration files for HNSCC in the TCGA were obtained from UCSC Xena (https://xena.ucsc.edu/). Pearson correlation analyses between risk score and immune cell infiltration scores were conducted by the R packages "limma", "scales", "ggplot2", "reshape2", "tidyverse", "ggpubr" and "ggtext" and we screened out and visualized the results at p < 0.05. The TME scores of HNSCC were calculated using the ESTIMATE algorithm with estimate R package. Then we assessed immune cell infiltration and immune cell function by singlesample GSEA (ssGSEA), and we accessed immune checkpoint activity by comparing differential expression of immune checkpoint genes.
Prediction of drug therapy response We adopted the Genomics of Cancer Drug Sensitivity in Cancer (https://www.cancerrxgene.org) and R package pRRophetic to forecast the sensitivity to drug therapies as per the half-maximal inhibitory concentration (IC50) in the two groups. All the code files mentioned above can be downloaded and used at the following website: https:// github.com/yexian123/ARLs-HNSCC-analysis.git.
Quantitative real-time PCR (qRT-PCR) Total RNA from each cell line was extracted with an RNA-easy isolation reagent, and reversed with HiScript Ò III RT SuperMix for qPCR (+gDNA wiper) (both Vazyme Biotech, China) to synthesize complementary DNA (cDNA). The PCR system was made from SYBR Green Ò Premix Ex Taq TM (Vazyme Biotech, China). Results of individual lncRNAs were standardized to the GAPDH expression. The primer sequences for LINC00861 and GAPDH are listed in Table S2.
Cell counting kit-8 (CCK-8) assay Cell proliferation was detected using CCK-8 assay (Beyotime, Shanghai, China) according to the instructions. Cells (1000/ well) were planted into 96-well plates and grown in an RPMI-1640 medium containing 10% FBS. At the same time each day, 10 μl of a CCK-8 solution was added to each well, and further incubated for 2 h at 37°C. Absorbance at 450 nm was detected using a microplate spectrophotometer (Thermo, USA) and used to estimate the proliferative capacity of CNE1 and CNE2 cells.

EdU assay
We assessed cell proliferation using an EdU assay kit (RiboBio, China). CNE1 and CNE2 cells (5 × 10 4 cells/well) were inoculated in a confocal laser cuvette, then cultured in a medium with 50 μM EdU (C10310-1, RiboBio) for 2 h and treated as per the instructions. Images were obtained under a fluorescent microscope. The mean proportion of EdUpositive cells in three random fields of view was analyzed.
Senescence-associated β-galactosidase (SA-β-gal) staining Cellular SA-β-gal activity were assayed according to the manufacturer's instructions (C0602, Beyotime, Nanjing, China). CNE1 and CNE2 cells were fixed and stained with the fixative solution and staining solution mixture provided in the kit. Observe and photograph the cells under an ordinary light microscope. The number of SA-β-gal-positive cells in 3 to 5 regions of the 6-well plate was randomly counted as a percentage of the total number of cells, and the mean value was calculated.

Statistical analysis
All statistical analysis and graphical visualization were conducted on R software 4.1.2. Categorical data between groups were compared via the chi-square test. Continuous variables in normal distribution (containing risk score and TME scores) were compared between two or more groups via Student's t-test or one-way ANOVA. Differences in expression levels between immune checkpoint genes and ARLs were investigated via Wilcox test. The significance level was probability less than 0.05. Fig. 1 illustrates the flowchart of the study. Firstly, we downloaded transcript data of 44 normal and 504 cancer samples from TCGA on 2022-04-15. Then 307 aging-related genes (ARGs) were obtained from HAGR (Table S1). After that, we got 156 AGs and 565 ALRs by correlation analysis (correlation coefficient > 0.5 and p < 0.001).

Identification of aging-related LncRNAs in HNSCC patients
The regulatory network of these lncRNAs with AGs is shown in Fig. 2A. Based on differential expression analysis of the ALRs in cancer versus normal tissue samples (|Log 2 FC| > 1 and p < 0.05), we acquired 276 significantly different ALRs, of which 244 were up-regulated and 32 were down-regulated in cancer tissues (Table S3, Fig. 2B).

Establishment and verification of the prognostic model
The 501 eligible tumor samples with complete information were equally divided into a train (N = 251) group and a test group (N = 250). In the training group, univariate Cox regression analysis identified 38 ALRs significantly related with overall survival (OS) (p < 0.05) (Figs. 3A and 3B). To avoid over-fitting of the prognostic features, we used LASSO regression of these lncRNAs and isolated sixteen lncRNAs associated with AGs in HNSCC when the first-order value of Log(λ) was the smallest likelihood of bias (Figs. 3C and 3D), where Log(λ) was −3.96. After that, we performed multi-COX regression analysis and finalized nine ARLs used to construct the model. The Sankey diagram (Fig. 3E) With the risk score formula, we divided the train, test, and entire groups into low-and high-risk groups. We compared the distribution of risk scores, survival status, lncRNA expression, and survival time among the risk groups. All results showed significant differences between the low-and high-risk groups (Figs. 4A-4D). Timedependent receiver operating characteristics (ROC) were used to assess the sensitivity and specificity of the model for prognosis. We also illustrated the ROC results in terms of the area under the ROC curve (AUC); the 1-, 3-and 5-year  Independent prognostic ability of the signature In univariate Cox regression analyses for all tumor samples, the hazard ratio (HR) for the risk score was 1.324, with a 95% confidence interval (CI) of 1.193-1.468 (p < 0.001). In addition, age, and stage were also statistically different (Fig. 5A). In the multivariate Cox regression, risk score (HR: 1.379, 95% CI: 1.225-1.553, p < 0.001), age (HR: 1.023, 95% CI: 1.007-1.038; p = 0.003) and stage (HR: 1.610, 95% CI: 1.322-1.960; p < 0.001) were independent prognostic factors of HNSCC (Fig. 5B).
Then based on risk score, age, and stage, we constructed a nomogram to predict the incidence rates of 1-, 3-, and 5-year OS in HNSCC patients (Fig. 5C). The calibration curves demonstrated a high concordance between the actual and nomogram-predicted survival rates of HNSCC patients (Fig. 5D). The sensitivity and specificity of the model for prognosis were evaluated via time-dependent ROC. The risk score (AUC = 0.684) and nomogram score (AUC = 0.715) had better predictive ability than age, gender, grade, and stage in the 3-year ROC of the risk model (Fig. 5E).
Gene set enrichment analyses and investigation of immune status Firstly, with the GSEA software, we explored the KEGG pathways in low-and high-risk groups to investigate the biofunction differences. Of the top 10 significantly enriched pathways, the high-risk group was enriched in only one pathway related to steroid biosynthesis. In comparison, the low-risk group was enriched in four immune-associated KEGG pathways (all FDR < 0.25) (Fig. 6A). Then we analyzed the relevance of immune cells to the low-and high-risk groups using data from multiple platforms. Results showed the majority of immune cells were related to the low-risk group, which generally agrees with the GSEA results (Fig. 6B). The content of B cells, follicular helper T cells (cTfh), was negatively related to the risk score, suggesting they were more abundant in the low-risk group (Fig. 6C), which had larger immune score and ESTIMATE score (Fig. 6D). Furthermore, the low-risk group had more immune cell infiltration such as B cells, CD8+ T cells, and dendritic cells (DCs), and more prosperous immune functions (Fig. 6E). The above results imply that the lowrisk group is under a higher immune infiltration state.

LINC00861 inhibited the proliferation of nasopharyngeal carcinoma cells
We selected LINC00861 for further analysis since it had the highest correlation with the risk signature. The qRT-PCR indicated that LINC00861 was differentially expressed in three different nasopharyngeal carcinoma cell lines ( Fig. 8A). We transfected the LINC00861-PCDNA3.1 plasmid in CNE1 and CNE2 cells to overexpress it and verified the transfection efficiency by qRT-PCR (Fig. 8B). CCK-8, Edu, and SA-β-gal Staining assays showed that LINC00861 overexpression considerably inhibited the growth and increased the senescence of CNE1 and CNE2 cells (Figs. 8C-8E).

Discussion
Surgery, radiotherapy, and standard chemotherapy are commonly used for HNSCC. Recently, medicine including cetuximab, an anti-EGFR, and pembrolizumab, an anti-PD-1 therapy, was accepted for late-stage, recurrent, or metastatic HNSCC [33,34]. Regardless of therapeutic progress, the 5-year OS rate of HNSCC is still low (<50%) due to its ease of invasion, metastasis, development of chemo-resistance, and the fact that most cases are diagnosed at late stages [35]. Moreover, the 3-year OS rate is around 80% for HPV+ HNSCC and 55% for HPV-HNSCC [36]. Recently, a variety of novel biomarkers such as PITX2 methylation and C1GALT1 autoantibody that forecast prognosis, and B7-H3 that predicts immunotherapeutic prognosis and reaction in HNSCC were identified gradually [37-39]. However, the prognostic models and biomarkers need to be updated since more new drugs, including ICP inhibitors, are put into clinical practice. Here, a new model based on ARLs was constructed for the first time in HNSCC, which can both predict prognosis and differentiate the TIME.
Aging is closely connected with the occurrence of many chronic diseases and tumors [5]. Cell senescence is a physiological process whereby cells lose the proliferative ability forever. Cell senescence is triggered by various stimuli and the DNA damage response (DDR), leading to stimulation of the p53 and/or p16INK4A pathways, and its important regulatory factors include the cell cycle regulators p53, p21, p16, and pRb proteins [40]. The increased activity of lysosomal senescence-related β-galactosidase (SA-β-Gal) is an important marker of senescent cells [41]. Cell senescence is a double-edged sword in malignant tumors, as it prevents tumor development and promotes HNSCC progression through the secretion of various factors, called senescence-associated secretory phenotype (SASP). Notch, mTOR, and NF-κB pathways are engaged in SASP management [42]. The majority of studies conclude that cellular senescence acts primarily as a cancer suppressor rather than a cancer promoter [43-50]. Silencing peroxiredoxin1(Prx1)-induced cellular senescence blocks malignant conversion and inhibits the formation of oral squamous carcinoma [51]. In nasopharyngeal carcinoma, senescent cells promote tumor cell invasion by increasing the secretion of matrix metalloproteinase 9 (MMP 9) [52]. Bioinformatics analyses reveal the critical function of AGs in predicting the survival prognosis and revealing the TIME landscape of HNSCC [9,53].
Recently, many studies show that lncRNAs are critical in HNSCC occurrence, and can both promote and suppress cancer. More importantly, some lncRNAs can serve as biomarkers and treatment targets for HNSCC and can predict patient prognosis [2,54,55]. However, as the significant member involved in transcriptional and posttranscriptional regulation, the adjusting role of aging by lncRNAs in HNSCC is unknown. Therefore, we analytically validated the link of lncRNAs with aging and the prognosis of HNSCC, aiming to make some contribution to subsequent studies.
By LASSO and multi-cox regression analysis, we identified nine lncRNAs associated with eight AGs in HNSCC, all of which were positively regulated with the genes (Fig. 3E). Three lncRNAs (LINC01508, AL359502.2, AC243960.1) were related with a high risk of patients with poor prognosis, and the remaining 6 lncRNAs (AC246787.2, LINC00996,`CDKN2A-DT`, Z97653.1, LINC00861, LINC02384) were associated with low risk. These nine lncRNAs were enrolled into a prediction model. Through validation, the model was proved to be reliable and stable (Figs. 4, 5A-5E). By comparing with published prognostic models of HNSCC constructed based on ferroptosis-related lncRNAs, we found that these articles had high AUC predictive values for 3-year survival (0.83, 0.811, 0.726, and 0.687, respectively) [56][57][58][59], the 3-year AUC in our results was 0.779 for the train set, so collectively it seems that our proposed model shows a good accuracy in predicting prognosis. Reportedly, LINC01508 upexpression inhibited cisplatin tolerance of ovarian cancer cells by suppressing the Hippo-YAP pathway [60]. In cervical cancer cells, LINC00861 inhibits tumor progression and EMT by acting as a ceRNA for miR-513b-5p and modulating the PTEN/ AKT/mTOR pathway [61]. The biological functions of the rest lncRNAs used in the modeling are unknown, about which researchers can carry out future experiments.
Aging not only regulates cell proliferation, but also alters the TIME through the secretion of SASP, which contains MMP, growth factors (VEGF and GM-CSF), inflammatory factors (IL-6 and IL-8), and the accumulation of lipofuscin deposits [62]. GSEA implied that the low-risk group was enriched with immune pathways such as the B cell receptor (BCR) and T cell receptor (TCR) pathways, and had higher immune infiltration and immune score, which may be related to their better prognosis (Figs. 6A and 6D). Considering the function of cell senescence in regulating TME, we conducted ssGSEA to explore the immune state in the low-and high-risk groups. The immune cells (aDCs, B cells, CD8+ T cells, DCs, iDCs, mast cells, neutrophils, NK cells, pDCs, T helper cells, Tfh cells, Th1 cells, Th2 cells, TIL, Tregs) and immune roles (APC (or T cell) cosuppression, APC (or T cell) co-stimulation, CCR, checkpoint, cytolytic activity, HLA, proinflammation, type II IFN response) were more active in the low-risk group (Figs. 6B and 6E). ICPs including PD-1, CTLA4, CD200, CD27, CD86, CD80, CD70, and TIGIT were more expressed in the low-risk group (Fig. 7A). These results imply that we can use the risk scores of ARLs to well distinguish the immune infiltration status of HNSCC and which is higher in the low-risk group. SASP factors released by senescent cells recruit innate immune cells (neutrophils, NK cells) and adaptive immune cells (CD8+ T cells) to modulate the removal of senescent tumor cells [63]. Therefore, the low-risk group may have more active cell senescence and thus is can be more cured by immunotherapy.
Radio-, chemo-, immuno-, and targeted therapies can induce cellular senescence and play an important role in TME through SASP [64]. Most HNSCC patients receive radiotherapy, which alters TME and activates the immune response of tumor cells [65,66] and may account for better outcomes in the low-risk group. As for drug sensitivity, the low-risk group was more sensitive to six inhibitors of the PI3K/AKT/mTOR pathway, such as rapamycin, which not only enhanced immunotherapeutic effects [18,22-25,29] but also blocked SASP-induced tumor progression (Fig. 7B). One study found that rapamycin prevented radiationinduced secretion of NFκB-driven pro-inflammatory SASP from prostate cancer cells and fibroblasts, thereby inhibiting tumor progression [67]. As a result, our study can guide clinical treatment and provide evidence for improving future immunotherapy and seeking appropriate target populations. Moreover, LINC00861 was identified as a protective factor against HNSCC progression in vitro (Fig. 8).
However, our study still has some limitations. Firstly, both our training and validation sets were from TCGA, which may cause bias, so the results may be more plausible if external validation was performed. Secondly, we only performed qRT-PCR, CCK-8, Edu, and SA-β-gal staining assays validation of LINC00861 in cell lines, but did not do further functional experiments or tissue validation, so further experiments are still needed. Finally, more clinical follow-up data shall be used to demonstrate the value of our prognostic model.
In conclusion, we built a robust prognostic prediction model for HNSCC using nine ARLs. It has a good predictive value in survival time, immune infiltration, immune checkpoint expression, and sensitivity to multiple drugs. LINC00861 acts as a protective factor against HNSCC progression. Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.