RNA-sequencing (RNA-seq) transcriptomic data and the accompanying clinical information for HNSCC and normal control samples were obtained from the Cancer Genome Atlas (TCGA) (https://gdc.cancer.gov/) (23). RNA-seq data were normalized using the Expectation-Maximization method. A total of 502 HNSCC and 44 healthy subjects were subjected to randomized analysis. TCGA transcriptomic data were used to obtain accurate and effective gene expression matrices using Perl language. An expression matrix of mRNAs and lncRNAs was additionally obtained using Perl language. The ‘limma’ package in ‘R’ was utilized to determine differential expression of ferroptosis genes and mRNAs using the following screening conditions: P<0.05, false discovery rate <0.05 and log2 fold change >1. Volcano plots and heatmaps of ferroptosis-related mRNAs were generated. ‘Corfile >0.4, P<0.01’ was set as the screening standard and the relationship between differentially expressed mRNAs and ferroptosis genes determined using R 4.1.3 software (The R Foundation) (24).

To determine the action sites of key mRNAs, a co-expression network of ferroptosis genes and mRNAs was established, with the aim of screening mRNAs significantly related to ferroptosis-related genes. Screening conditions were ‘P<0.01, coef >0.4’ and the network was established using R software.

The differentially expressed mRNA matrix was merged with clinical data downloaded from TCGA and a prognostic model constructed based on differentially expressed mRNAs. Patients with HNSCC were divided into training and test groups. The training group, accounting for 50% of the patients was used to establish a prognostic model, which was validated in the test group. Univariate Cox proportional hazards regression analysis was performed and survival-related mRNAs screened based on P-values <0.05. A least absolute shrinkage and selection operator (LASSO) regression model was generated using the screened mRNAs and 10-fold cross-validation performed with 1,000 cycles to identify the point with the smallest cross-validation error. Significantly expressed mRNAs were identified via LASSO regression and used to generate a Cox model.

A total of 15 pairs of tumor tissue (T) and adjacent non-tumor tissue (N) samples were obtained from the Department of Stomatology of the Affiliated Hospital of Xuzhou Medical University from patients with HNSCC without preoperative treatment (Xuzhou, China). Among these 15 patients, there were 9 female patients and 6 male patients, all aged between 40–70 years-old. The tissues were cut and subsequently stored in liquid nitrogen. Samples were independently identified by two pathologists and detailed clinicopathological features were collected. All tissue samples were obtained from July 2020 to June 2022. Informed consent was obtained from each patient prior to any specimen-related studies and the study was approved by the Ethics Committee of the Affiliated Hospital of Xuzhou Medical University (approval no. XYFY2020-KL216-01; Xuzhou, China). Among them, 7 pairs of cancer tissues and 3 cases of adjacent tissues were used for transcriptomics; 15 pairs of cancer tissues and adjacent tissues were applied for reverse transcription-quantitative PCR and western blot experiments.

Total RNA was isolated using TRIzol^® reagent from cancer tissues of seven patients with HNSCC and three adjacent tissue samples. Detection of the concentration and purity of extracted RNA was performed via spectrophotometry on a Nanodrop 2000 instrument (Agilent Technologies, Inc.); only high-quality RNA sample (OD 260/280=1.8~2.2, OD 260/230 ≥2.0, RIN ≥6.5, 28S:18S ≥1.0, >1 µg) was used to construct sequencing library. Agarose gel electrophoresis was used to determine the integrity of RNA, and RNA integrity number values were established with the Agilent 2100 Bioanalyzer system (Agilent Technologies, Inc.). RNA-seq transcriptome library was prepared following TruSeqTM RNA sample preparation Kit from Illumina, Inc. using 1 µg of total RNA. Shortly, messenger RNA was isolated according to polyA selection method by oligo(dT) beads and then fragmented by fragmentation buffer firstly. Secondly double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen; Thermo Fisher scientific, Inc.) with random hexamer primers (Illumina, Inc.). Then the synthesized cDNA was subjected to end-repair, phosphorylation and ‘A’ base addition according to Illumina's library construction protocol. Libraries were size selected for cDNA target fragments of 300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quantified by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq xten/NovaSeq 6000 sequencer (2×150 bp read length). according to the manufacturer's protocol. Transcriptomes from all samples were merged to reconstruct a comprehensive transcriptome using perl script. Normalization was performed using R software (v4.1.3) and differentially expressed genes were identified by their P-value and expression fold change.

CAL-27 (RRID: CVCL_0030) and HOK (RRID: CVCL_YE19) cells purchased from the Procell Life Science & Technology Co., Ltd. were cultured in DMEM/F12 (1:1) medium (Gibco; Thermo Fisher Scientific, Inc.) containing 10% heat-inactivated FBS (Gibco; Thermo Fisher Scientific, Inc.) and 1% penicillin and streptomycin. Cells were cultured in a constant temperature incubator under 5% CO₂ at 37°C. Erastin (cat. no. HY-15763) was obtained from MedChemExpress.

(RT-qPCR) was used to detect Senp1 transcription levels. Total RNA was extracted from cells and tissues using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc). Extracted RNA was reverse-transcribed into cDNA with HiScript^® II QRT Super Mix for qPCR (Vazyme Biotech Co., Ltd.) and RT-qPCR was performed using SYBR Green qPCR Master Mix (MedChemExpress) on the Roche Light Cycler 480 II instrument (Roche Diagnostics) using β-actin as an internal standard. Relative gene expression was normalized to β-actin mRNA expression and analyzed using the comparative Cq method (ΔΔCq) (25). The same thermocycling profile conditions were used for all primers sets: 95°C for 10 min, 40 cycles at 95°C for 15 sec, and 60°C for 1 min. The primers used were as follows: β-actin forward, 5′-CTCCATCCTGGCCTCGCTGT-3′ and reverse, 5′-GCTGTCACCTTCACCGTTCC-3′; SENP1 forward, 5′-CCAGCATTTTAACTAACCAGGAAC-3′ and reverse, 5′-GCTAAGTTATCTGGCTGATGTGG-3′; and ACSL4 forward, 5′-CATCCCTGGAGCAGATACTCT-3′ and reverse, 5′-TCACTTAGGATTTCCCTGGTCC-3′.

Proteins were extracted from tissues and CAL-27 cells using RIPA lysis buffer (cat. no. 89900; Thermo Fisher Scientific, Inc.) with Halt™ Protease Inhibitor Cocktail (cat. no. 78429; Thermo Fisher Scientific, Inc.). Subsequently, the BCA method was used to quantify the tested protein. Premade adhesive gels (4–20%) were prepared and 30 µg protein samples were loaded per lane, marked on both sides to indicate protein molecular weight and to distinguish left and right; then separated via sodium dodecyl sulfate-polyacrylamide gel electrophoresis and transferred to polyvinylidene fluoride membranes. Next, membranes were blocked with 5% non-fat milk at room temperature for 2 h, followed by overnight incubation at 4°C with anti-SENP1 (1:1,000; cat. no. 25349-1-AP), anti-GPX4 (1:1,000; cat. no. 67763-1-Ig), anti- solute carrier family 7 member 11 (SLC7A11) (1:1,000; cat. no. 26864-1-AP) and anti-β-actin (1:5,000; cat. no. 20536-1-AP; all from Proteintech Group, Inc.) antibodies. After washing, membranes were incubated with goat anti-mouse (1:5,000; cat. no. ab216776) or goat anti-rabbit secondary antibody (1:5,000; cat. no. ab216773; both from Abcam) for 2 h at room temperature. An enhanced chemiluminescence reagent (Wuhan Servicebio Technology Co., Ltd.) was used to visualize the protein bands. Protein bands were scanned with an Odyssey infrared scanner (LI-COR Biosciences) and band density was analyzed with ImageJ v1.8.0 software (National Institutes of Health).

The CCK-8 assay (Beyotime Biotechnology Inc.) was employed to detect cell activity. In brief, CAL-27 cells were seeded into a 96-well plate at a density of 1×10⁴ cells/well and cultured for 3 days. Next, cells were transfected with Senp1 short hairpin RNA (shRNA) or overexpression plasmid and survival rates were assessed via CCK-8 colorimetry at 72 h. In brief, CCK-8 was added to each well and was incubated for 1 h at 37°C. CCK-8 solution accounted for ~10% of the volume of the culture medium. OD values were measured at 450 nm using a microplate reader (BioTek Instruments, Inc.).

Transfected cells were cultured to ~90% confluence in 12-well plates and a wound area introduced using a sterile 200-µl pipette tip. After washing twice with PBS, cells were cultured in DMEM/F12 without FBS for 48 h. Cell images were obtained using fluorescence microscopy (IX73; Olympus Corporation; scale bars, 100 µm). Wound sizes were measured with ImageJ and calculated as a percentage (from 0 h to 48 h).

Cells were collected and resuspended in medium (~1×10⁵ cells/ml). An aliquot (200 µl) of the cell suspension was placed in the upper chamber of a Transwell (8-µm pore size; MilliporeSigma), which was inserted into a 24-well plate. After incubation in RPMI-1640 containing 10% (v/v) FBS for 24 h at 37°C, cells on the upper membrane of the insert were removed gently with a cotton swab while those migrating to the bottom surface were fixed with 4% cold methanol for 30 min and stained with 0.1% crystal violet for 1 h. The number of cells on the bottom surface was counted and quantified using fluorescence microscopy (IX73; Olympus Corporation; scale bars, 100 µm).

The TBA method kit (cat. no. A003; Nanjing Jiancheng Bioengineering Institute) was used to measure the concentration of MDA in cells. Cells were mixed well in RIPA lysis buffer. The reagent was added to the centrifuge tube according to the manufacturer's instructions and thoroughly mixed, followed by boiling in a 95°C water bath for 40 min and cooling at room temperature. After centrifugation at 1,600 × g for 10 min at 4°C, absorbance of the supernatant was measured at 532 nm to calculate malondialdehyde (MDA).

An iron assay kit (cat. no. ab83366; Abcam) was employed to evaluate the relative intracellular iron levels. According to the requirements of the reagent supplier, the prepared iron probe, reducing agent, standard diluent (2–10 nmol/well) and test sample (cell lysate) were mixed. Next, 5 µl iron reducing agent was added to each standard well and incubated at 37°C for 30 min, followed by addition of 100 µl iron probe to each well at 37°C and incubation in the dark for 60 min. The sample was immediately evaluated with a colorimetric enzyme marker (593 nm). The total iron content of test samples was determined from the standard curve.

Slices were immunohistochemically stained as experimental procedure. Coverslips seeded with tumor cells or tissue cryo-sections were fixed with 4% paraformaldehyde for 30 min, and then permeabilized with 0.5% Triton X-100 for 10 min. After 2 h of incubation with goat serum (cat. no. 16210064; Thermo Fisher Scientific, Inc.) at 37°C, Consecutively, 4-µm-thick sections were analyzed using primary antibodies against SENP1 (1:1000; cat. no. 25349-1-AP; Proteintech Group, Inc.) at 4°C. Then, sections were incubated with Alexa Fluor (cat. no. R37116; Thermo Fisher Scientific, Inc.) secondary antibody for 1 h at 37°C. Images were obtained using fluorescence microscopy (IX73; Olympus Corporation; scale bars 200 µm).

A GSH colorimetric assay kit (cat. no. CS0260; MilliporeSigma) was used to determine relative GSH levels. This experiment used 5% of 5-thiosalicylic acid solution (SSA), GSH standard solution, working mixture, NADPH solution and test sample (cell extract). SSA (5%) was used as the reagent blank, a total of 10 µl GSH standard solution and 10 µl of test sample were added to a separate well. An aliquot of 150 µl working mixture was added to each well and was incubated at room temperature for 5 min. Finally, a total of 50 µl NADPH solution was added to each well and mixed. After absorbance of each well was determined at 412 nm, the GSH concentration was calculated according to the standard curve.

To construct stable SENP1-knockdown CAL-27 cell line, shRNA specifically targeting SENP1 was designed by the Thermo Fisher Scientific, Inc. online tool (https://rnaidesigner.thermofisher.com/rnaiexpress/construct.do?pid=−431479708000205999) and cloned into the lentiviral pLVX-sh1 vector. SENP1 overexpressing CAL27 cells were cloned into the lentivirus pLVX–IRES-puro vector. Lentiviruses were synthesized by Genechem Shanghai Co., Ltd. Cells were seeded in cell culture dishes and used to transfect lentiviral with different MOI values (10, 20, 30) at 30% density. On the second day of transfection, cells were cultured in DMEM containing 10% FBS. After 3 days, the transfection effect was observed under fluorescence microscopy (IX73; Olympus Corporation). Then, transfected cells were passaged and cultured with puromycin (5 ng/µl; Shanghai Yeasen Biotechnology Co., Ltd.) used to selected cells for 1 week. Subsequently, stable cell lines were cultured in DMEM containing 10% FBS and puromycin (2 ng/µl). The sequences were as follows: sh-SENP1, 5′-CTCGATGTCTTAGTTCCAGTA-3′; and sh-NC, 5′-GGGACTCGATTCCCACGAATC-3′.

The cells were lysed with 250 µl NP-40 lysis buffer (cat. no. P0013F; Beyotime Institute of Biotechnology). The cell lysates were centrifuged at 4,800 × g for 5 min to remove the cell fragments. The protein stock solution, with anti-SENP1 or ACSL4 antibody bound to 25 µl magnetic beads, was incubated overnight at 4°C. The small pallet was washed three times using the pyrolysis method and sodium dodecyl sulfate samples were boiled for 5 min. SENP1, ACSL4 and SUMO1 levels were analyzed via western blot analysis using anti-SENP1 (1:1,000; cat. no. 25349-1-AP), anti-ACSL4 (1:1,000; cat. no. 22401-1-AP) and anti-SUMO1 (1:1,000; cat. no. 10329-1-AP; all from Proteintech Group, Inc.) antibodies, respectively.

All data were analyzed with Graph Pad Prism 7.0 software (GraphPad Software; Dotmatics). Continuous data are presented as the mean ± standard deviation and differences between experimental groups were analyzed using one-way ANOVA or unpaired two-tailed Student's t-test. In addition, to estimate the risk differences between the main outcome groups, a post hoc analysis was conducted using Bonferroni's test following ANOVA. Survival analysis was performed using the Kaplan-Meier method(https://ualcan.path.uab.edu/). The screening of key genes comes from the VENN diagram method (https://bioinfogp.cnb.csic.es/tools/venny/index.html). Differences were considered statistically significant at P<0.05. Experiments were repeated three times and representative experiments are shown.

Data from a total of 44 normal and 502 tumor samples were obtained from TCGA database. Based on the co-expression patterns of 40 ferroptosis genes (Table SI) and differentially expressed mRNAs (|log2FC|>1 and P<0.01), a total of 457 ferroptosis-related mRNAs were obtained (correlation coefficients >0.4 and P<0.001). All mRNAs and ferroptosis genes were positively correlated. Volcano plots and heat maps were used to examine expression patterns in normal and tumor groups (Fig. 1A and C) and the relationship between mRNAs and ferroptosis genes was ultimately visualized via network graphs (Fig. 1B).

Using univariate Cox proportional hazards regression analysis, 32 significant ferroptosis-related mRNAs among the 457 were identified and a heatmap was generated (Fig. 2A and B). The significant mRNAs obtained via univariate Cox analysis were combined with clinical data to construct LASSO regression curves (Fig. 2C). Cross-validation was performed and the point of the smallest cross-validation error determined. Ultimately, 13 ferroptosis-related mRNAs were screened (Fig. 2D). Following optimization of the Cox model, seven mRNAs (SENP1, PBX4, SOCS1, DNAJB11, ART5, TCEAL5 and GRIA4) were selected. Simultaneous transcriptome analysis of HNSCC and corresponding paracancerous tissues led to the identification of 427 differential genes (199 upregulated and 228 downregulated) (Fig. S1A). When these 427 genes were analyzed with the VENN map with the aforementioned seven genes, a key gene, SENP1, was obtained (Fig. S1B). Thus, SENP1 was selected to start further verification.

According to the risk score, a Kaplan-Meier curve of overall survival (OS) probability of individuals with ferroptosis-related genes was constructed for the training, test and all groups. The survival rate was significantly lower in high risk than low risk patients for all groups (Fig. 3D). Similar results were obtained for age, sex, stage, type and tumor-node-metastasis (TNM) stage (Fig. S2). Since HNSCC is prone to lymph node metastasis, ‘N’ was specifically taken out as an independent clinical prognostic factor. Effective independent prognostic factors were screened via univariate Cox analysis (screening condition, P<0.05). The risk score of the high-risk group was obviously greater whereas survival time was markedly lower than that of the low-risk group in the training, test and all groups. The differences in expression levels of the seven ferroptosis-related mRNAs in prognostic models of the high risk as well as low risk groups were analyzed using heat maps (Fig. 3A-C).

As aforementioned, effective independent prognostic factors were screened via univariate Cox analysis (screening condition, P<0.05). The hazard ratio (HR) for the risk score was 1.732 and 95% confidence interval (CI) was 1.476–2.032. ‘Stage’ and ‘N’ were identified as two additional prognostic parameters (Fig. 4A). In multivariate Cox analysis, HR of the risk score was 1.623 with 95% CI of 1.373–1.918. Age, stage and N could be used as additional prognostic parameters (Fig. 4B). A nomogram was generated according to the clinical risk factors identified, specifically, age, sex, stage, grade, TNM stage and comprehensive risk scores were obtained to assess 1-, 3- and 5-year survival rates of patients with HNSCC (Fig. 4C). Calibration curves for 1-, 3- and 5-year survival rates were additionally generated to evaluate the reliability of the nomogram (Fig. 4D). Receiver operation curve (ROC) analysis was utilized to validate the preciseness of the model for predicting patient survival. The predictive values of area under the curve (AUC) were 0.672, 0.716 and 0.666 for 5-, 3- and 1-year survival, respectively (Fig. 4E). All AUC values were >0.65, validating the accuracy of the prognostic model in survival prediction. Moreover, AUCs were larger than those obtained from other clinical traits, suggesting that accuracy of the proposed model was greater than that of other clinical traits (Fig. 4F). The data clearly suggested that the seven mRNA ferroptosis-related models constructed by the authors were closely associated with the pathogenesis of HNSCC.

Kaplan-Meier method indicated that high SENP1 expression was generally associated with low survival rates (Fig. S3). SENP1 mRNA that displayed the highest association with ferroptosis was selected for further verification (Fig. S4). Subsequently, the expression of SENP1 in HNSCC tissue was examined via IHC, western blotting and RT-qPCR and was revealed to be significantly higher than that in paracancerous tissue (Fig. 5A-D). Simultaneously, the expression of SENP1 in CAL-27 and HOK was examined via RT-qPCR and western blot analyses. Both sets of experiments revealed higher expression of SENP1 in the CAL-27 cell line (Fig. 5E-G). The results clearly indicated that SENP1 was overexpressed in HNSCC.

Stable knockout or overexpression cell lines using shRNA and SENP1 overexpression lentivirus, respectively, were first established with a view to investigate the biological role of SENP1 in HNSCC (Fig. 6A). Overexpression and knockdown efficiency of SENP1 was validated at both protein and transcriptional levels (Fig. 6C-E). Data from the CCK-8 assay revealed increased CAL-27 cell proliferation under conditions of overexpression of SENP1 (Fig. 6B). Conversely, knockdown of SENP1 significantly inhibited the viability of CAL-27 cells, suggesting that SENP1 was an important regulator of HNSCC cell proliferation. Transwell assays revealed that compared with the control group, the number of migratory HNSCC cells transfected with sh-SENP1 decreased while upon overexpression of SENP1, the number of migratory HNSCC cells was significantly increased (Fig. 6F and G). The results of the wound healing assay consistently demonstrated that SENP1 promoted cell migration. Compared with the negative control group, the migration ability of oe-SENP1 CAL-27 cells was significantly increased while sh-SENP1 induced a decrease in cell migration ability (Fig. 6H and I). The collective findings indicated the pivotal role of SENP1 in the progression of HNSCC.

Erastin is a common inducer of ferroptosis. Following activation of ferroptosis of HNSCC cells with erastin, changes in SENP1 expression and its inhibitory mechanisms of action on growth of HNSCC cells were examined. Initial RT-qPCR and western blotting experiments revealed significantly reduced SENP1 expression in cells treated with erastin compared with those treated with DMSO (Fig. 7A-C). To further determine whether inhibition of cell growth was caused by ferroptosis, under conditions of erastin-induced ferroptosis (10 µM, 24 h), the levels of iron death-related indicators (cellular iron, GSH and MDA) in HNSCC cells were detected. Notably, viability of HNSCC cells and levels of GSH were reduced (Fig. 7D and F), whereas levels of iron and MDA, the end product of lipid peroxidation, were significantly increased in cells (Fig. 7E and G). Simultaneously, overexpression and knockdown efficiency of SENP1 at the protein and transcriptional levels were validated under conditions of ferroptosis induction (Fig. 8A-C). Following knockdown of SENP1, compared with the NC group, cell proliferation ability reduced, the end products of cellular iron and lipid peroxidation, MDA, were significantly increased while that of endogenous GSH was significantly decreased (Fig. 8D-G). In addition, the corresponding level of ferroptosis under overexpression of SENP1 was reduced (Fig. 8F). The findings of the aforementioned experiments clearly indicated that SENP1 served as an important regulatory factor of ferroptosis in HSNCC cells.

ACSL4 is a key enzyme in the lipid peroxidation procedure that is highly expressed to promote ferroptosis. While previous research has revealed that SENP1 regulates the process of ferroptosis, the association between SENP1 and ACSL4 has not been fully elucidated. It was proposed by the authors that SENP1 reduces the stability of ACSL4 protein via deSUMOylation in HNSCC. Co-IP experiments in CAL-27 cells revealed a significant direct association between SENP1 and ACSL4 (Fig. 9A). Further TCGA and western blot analysis disclosed a negative association between the two molecules (Fig. 9B and C). To further validate this hypothesis, IP experiments on SUMO1 conjugated to ACSL4 in the control and gene knockout groups were conducted. Notably, the SUMO1 level of ACSL4 was decreased in the gene knockout group (Fig. 9D). Furthermore, the impact of SENP1 knockdown on the iron death pathway was evaluated. Protein immunoblotting was employed to determine the levels of key proteins, GPX4 and SLC7A11, of the iron death pathway in samples from each group following SENP1 knockdown, which revealed that the GPX4 level was decreased while that of SLC7A11 remained unchanged (Fig. 9E-G). These data supported a direct correlation between SENP1 and ACSL4. Knockdown of SENP1 may promote SUMO1 conjugation and stability of ACSL4, further consuming GPX4 levels and ultimately, facilitating ferroptosis.