The original GEO datasets GSE37991 and GSE65858 were downloaded from NCBI GEO databases. The data of GSE37991 and GSE65858 were based on GPL6883 (Illumina HumanRef-8 v3.0 expression beadchip) and GPL10558 (Illumina HumanHT-12 V4.0 expression beadchip), respectively. For the TCGA HNSCC cohort, the raw data and corresponding clinical information were downloaded from the TCGA data portal. The tumor located in lips, tongue, oral cavity, oropharynx, larynx, and hypopharynx were selected. The format of the TCGA HNSCC downloaded data was HTseq-FPKM, and the formats of GSE37991 and GSE65858 were normalized microarray data.

The RNA‐seq data of the TGCA HNSCC cohort included 521 HNSCC samples and 44 normal control samples. Forty‐three out of 44 normal samples were matched to the HNSCC tumor samples. Only one normal sample from the salivary gland was not matched to the tumor samples. Two hundred and seventy HNSCC tumor samples were included in the microarray data of GSE65858, which did not include normal control samples. The microarray data of the GSE37991 cohort included 40 HNSCC samples and 40 paired normal samples. Regarding the inclusion/exclusion criteria for the enrolled patients, both tumor and normal control samples in TCGA HNSCC and GSE37991 cohorts were included for the differential expression analysis. For the construction of risk signature, only the tumor samples in TCGA HNSCC, GSE37991 and GSE65858 cohorts were considered, and all normal control samples were excluded. The patients without clinical characteristics such as follow-up time, follow-up status, age, gender, clinical stage, T stage or N stage were excluded. Besides, the patients with the TX stage and NX stage were also excluded due to the disturbance to grouping. The clinical characteristics including age, gender, stage of the HNSCC tumor samples in the TCGA HNSCC discovery cohort and GSE65858 validation cohort were summarized in Tables S1 and S2 , respectively.

Briefly, the probes that have no expression in most of the samples were excluded. The FPKM data of the TCGA HNSCC cohort was converted to TPM format for further analysis. For GSE37991, data normalization was achieved with GeneSpring GX software (Agilent Technologies). For GSE65858, the data was first processed within the R/Bioconductor. Then the expression values were subjected to log₂-transformation and the normalization was performed using Robust Spline Normalization (RSN). Batch effects of expression BeadChips were corrected using ComBat.

The DEMEs in GSE37991 and TCGA HNSCC cohort were identified by the edgeR package. Absolute log₂FC>1 and p<0.05 were selected as the demarcation criteria based on Benjamini & Hochberg (BH) procedure. The common differentially expressed genes (DEGs) between GSE37991 and TCGA HNSCC cohort were identified by the intersect function in R.

Gene Ontology (GO) enrichment analysis was performed using the DAVID (the Database for Annotation, Visualization and Integration Discovery).

TCGA HNSCC cohort and GSE65858 were used as discovery and validation cohorts, respectively. DEMEs were demonstrated to be correlative with the OS of HNSCC by univariate Cox proportional hazards regression analysis, which adopted the statistics from the TCGA HNSCC cohort. Subsequently, the most optimum DEMEs were chosen by the approach of the least absolute shrinkage and selection operator (LASSO) regression. The independent DEME-based prognostic signature was determined by the multivariate Cox proportional hazards regression analysis. Afterward a risk score model was constructed with this independent DEME-based prognostic signature. The summation of each DEME’s score was computed as a risk score for each HNSCC patient: risk score = Σ i = 1 n   β i × E i (14). On basis of the median value of the risk scores, the TCGA HNSCC cohort was divided into a low-risk group and a high-risk group. We compared the OS between low- and high-risk groups and assessed the differential survival distinguished by clinicopathological parameters between low- and high-risk groups. Similarly, the GSE65858 validation cohort was classified into low- and high-risk groups according to the above risk score model built up by the TCGA HNSCC cohort. Besides, the OS and the survival distinction were likewise differentiated by clinicopathological parameters between low- and high-risk groups.

A nomogram model including the risk score and other clinicopathological indices was constructed. The calibration curves were used to assess the accuracy for predicting 1-year OS and 3-year OS of HNSCC.

Both SCC1 and SCC23 were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco, USA) supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin. Cells were grown in a 37°C, 5% CO₂ cell incubator in a humidified atmosphere. The cells were transfected with siTXNRD1 and siCTRL using Lipofectamine^® RNA iMAX Transfection Reagent (Invitrogen, Carlsbad, CA, USA).

The protein samples were separated on 4-20% SDS-PAGE gels and then transferred to polyvinylidene difluoride (PVDF) membranes. Following by blocking in 5% skimmed milk for 1 hr at room temperature, the blots were then probed with primary antibody against TXNRD1 (1:1000, Proteintech, Chicago, IL, USA) at 4°C overnight. The corresponding HRP-conjugated secondary antibody was used to incubate the membranes for 1 hr at room temperature. The signal was detected with ECL kits.

The cells were seeded into the 96‐wells of the plate at a density of 3,000 cells per well in 200 μl cell culture media. MTT solution (20 μl,5 mg/ml) was added to each well at the indicated time points and incubated for 4 hrs at 37°C. Dimethyl sulfoxide (DMSO) was added to dissolve the precipitate after removing the supernatant. A microculture plate reader (Tecan, Mannedorf, Switzerland) was used to measure the absorbance at 570 nm.

According to the manufacturer’s instructions, the Click‐iT™ 5‐Ethynyl‐2′‐deoxyuridine (EdU) Cell Proliferation Kit for Imaging (Invitrogen) was used to perform the EdU assay. Briefly, EdU was added to the cells and incubated at 37°C for 2 hrs. Subsequently, 3.7% formaldehyde was used to fix the cells at room temperature for 20 min. After washing three changes of PBS, 0.5% Triton X‐100 was added to increase the permeability of the cellular membrane. 1× Click‐iT reaction cocktail was used to stain the cells in the dark at room temperature for 30 min. Then, the cell nucleus was stained by Hoechst 33342 dye. A confocal laser scanning microscope (Olympus, Center Valley, PA) was used to capture at least four random images per well.

The volcano plot and heatmaps were drawn by the “ggplot2” package of R software. The univariate and multivariate Cox regression analyses incorporated the clinical characteristics including age, gender, clinical stage, and risk score. The independent prognostic factors for HNSCC were identified by the univariate and multivariate Cox regression analyses. The Kaplan‐Meier method and log‐rank test were performed to calculate the OS distinguish between different groups. A p-value less than 0.05 is considered statistically significant.

The volcano plot was used to visualize the distribution of metabolic enzymes between cancer and normal tissues from the GSE37991 and TCGA HNSCC cohort. The significantly downregulated or upregulated metabolic enzymes were represented as green or red dots, respectively ( Figure 1A ). In total, 478 (402 upregulated and 76 downregulated) and 223 (102 upregulated and 121 downregulated) significantly changed metabolic enzymes were identified in GSE37991 and TCGA HNSCC cohort, respectively. The detailed information of the significantly changed metabolic enzymes was summarized in Tables S3 and S4 . The common DEMEs (23 upregulated and 9 downregulated) between GSE37991 and TCGA HNSCC cohort were shown in Figure 1B .

Gene ontology (GO) analysis of the DEMEs showed that small molecule catabolic process, alpha-amino acid metabolic process, glycosyl compound metabolic process, cellular amino acid metabolic process, organophosphate catabolic process, protein tetramerization, nucleoside metabolic process, nucleobase-containing small molecule catabolic process, aspartate family amino acid metabolic process, purine-containing compound catabolic process were the enriched ( Figure 2A ).

The survival-related DEMEs in the TCGA HNSCC cohort were identified by the univariate Cox regression, and ASNS, CYP27B1, TXNRD1, GATM, PLOD2, FUT6, and HPRT1 were harvested. Based on the HRs, GATM and FUT6 were protective genes, while ASNS, CYP27B1, TXNRD1, PLOD2, and HPRT1 were risky genes ( Figure 2B ). The LASSO regression analysis identified six optimal DEMEs including ASNS, CYP27B1, TXNRD1, PLOD2, FUT6, and HPRT1. The risk score for each patient was computed as follows: risk score = (0.308) *ASNS + (0.228) * CYP27B1 + (0.284) * TXNRD1 + (0.174) * PLOD2 + (-0.092) * FUT6 + (0.362) * HPRT1. Subsequently, the HNSCC patients were divided into high- and low-risk groups based on the median value of risk scores ( Figure 3A ). The survival time and survival status of each HNSCC patient in TCGA HNSCC were presented in a scatter plot ( Figure 3B ). The expression levels of ASNS, CYP27B1, TXNRD1, PLOD2, FUT6, and HPRT1 in each HNSCC patient were shown with a heatmap ( Figure 3C ). Besides, survival analysis demonstrated that the OS of the high-risk group was significantly lower compared to the low-risk group ( Figure 3D ). We then stratified the HNSCC patients with different clinical indices including age, gender, clinical stage, T stage, and N stage. As shown in Figures 4A–E , the low-risk group got remarkably better OS than the high-risk group for the HNSCC patients with age>60 (p=0.003), or with female gender (p=0.023) or with male gender (p=0.022), or at the clinical stage III-IV (p<0.001), or at the stage T3-4 (p=0.003), or with node metastasis (p=0.005).

Similarly, the HNSCC patients in GSE65858 were divided into high- and low-risk groups using the same median risk score in the TCGA HNSCC cohort ( Figure 5A ). The survival status, survival time, and expression level of prognosis-related DEMEs in each HNSCC patient were shown in Figures 5B, C . More importantly, the HNSCC patients in the high-risk group suffered a significantly poorer OS than those in the low-risk group ( Figure 5D ). As displayed in Figures 6A–E , the low-risk group got remarkably better OS than the high-risk group for the HNSCC patients with age ≤ 60 (p=0.026), or with male gender (p=0.009), or at the clinical stage III-IV (p=0.027), or the stage T3-4 (p=0.030).

In the TCGA HNSCC cohort, the univariate Cox regression analysis showed that age (p =0.014, HR=1.020), gender (p =0.045, HR=0.695) and risk score (p<0.001, HR=2.116) were significantly associated with survival ( Figure 7A ). The multivariate Cox regression analysis revealed that only risk score (p<0.001, HR=2.047) was the independent prognostic factor for HNSCC ( Figure 7B ). Similarly, in GSE65858, the univariate Cox regression analysis showed that age (p =0.013, HR=1.027), clinical stage (p=0.001, HR=1.615) and risk score (p=0.024, HR=1.715) were significantly correlated with survival ( Figure 7C ). The multivariate Cox regression analysis showed that age (p =0.015, HR=1.027), clinical stage (p=0.001, HR=1.637) and risk score (p=0.023, HR=1.710) were independent prognostic factors for HNSCC ( Figure 7D ).

The risk score, age, gender and clinical stage were incorporated into the nomogram model to forecast the clinical outcome of HNSCC ( Figure 8 ). A total nomogram-based score was obtained for each HNSCC patient derived from the clinicopathological parameters and their corresponding points. The 1-year OS or 3-year OS of HNSCC patients was forecasted with the nomogram model. The calibration curves showed that the nomogram model we built up exhibited excellent conformance for predicting the 1-year OS or 3-year OS of HNSCC ( Figures 9A, B ). The C indices of the nomogram model are 0.658 or 0.616 for predicting the 1-year OS or 3-year OS of HNSCC, respectively.

The expression level of TXNRD1 was significantly reduced in HNSCC cells following siTXNRD1 treatment ( Figure 10A ). Compared to the siCTRL treated cells, the MTT assay showed that the optical density (OD) values were markedly lower in TXNRD1-knockdown cells at indicated time points ( Figure 10B ). Similarly, the EdU assay showed the percentage of EdU‐positive cells was dramatically lower in HNSCC cells treated with siTXNRD1 ( Figure 10C ).