MSK-QLL1 and SCC15 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM): Nutrient Mixture F-12 (Gibco) supplemented with 10% fetal bovine serum (Gibco) and 1% penicillin/streptomycin (Gibco). All cells were cultured in a cell culture incubator at 37 °C with 5% CO2. The cells were irradiated with 6 Gy each time, then cultured for 2–3 days until they reached 70%–80% confluency, followed by re-irradiation. This procedure was repeated 10 times. Finally, a clonogenic assay was conducted under 4 Gy irradiation to compare parental cells with repeatedly irradiated cells to confirm the acquired radioresistance of the repeatedly treated cells (RR-MSK-QLL1 and RR-SCC15).

DMEM/F12 was treated using an NTP jet device. The device was composed of a mass flow controller, gas distributor, DC (direct current) power supply, plasma inverter, and NTP reactor (Fig. 1). The mass flow controller controlled the gas flow with helium (4 standard liters per minute) and oxygen (1 standard cubic centimeter per minute). The two gases were mixtures in a gas distributor and entranced in the NTP reactor. The NTP reactor was powered by a DC power supply and plasma inverter, delivering energy to two electrodes at a few kilovolts and 20 kHz, with the discharge power ranging from 10 W to 24 W. The NTP reactor comprised a quartz tube (outer diameter, 6 mm; inner diameter, 4 mm) with two electrodes, an inner stainless-steel tube electrode, and an outer ground ring electrode. The gas was converted into a plasma jet via electrodes and directed vertically into the DMEM/F12 media inside the flask. The distance between the plasma device and the medium was maintained at approximately 1 cm. The plasma’s emission spectra were reported in our previous study with hydrogen peroxide and nitric oxide (NO₃^–, NO₂^–) [18]. Before each experiment, DMEM/F12 was exposed to the plasma jet for 2 hours to generate fresh NTP-activated media (NTPAM). Then, radioresistant cells were treated with 100% NTPAM for 6 and 12 hours, with the untreated control as the baseline for comparison.

Parent and RR-HNSCC cells were seeded in 60 mm dishes with 3×10⁴ cells. On the second day, the cells were irradiated with 4 Gy. After irradiation, irradiated cells and no-irradiated control were reseeded at 1,000 cells/well in 12-well plates. Each well was supplemented with 2 mL of medium, which was replaced every 2–3 days. After 10–14 days, cells were washed with cold PBS, fixed with 4% paraformaldehyde (Thermo Scientific), and stained with 1% crystal violet (Sigma-Aldrich). Finally, images were captured, and cell counting was performed using the ImageJ software. The vitality of the cells was tested to demonstrate their radiation resistance.

Parent and RR-HNSCC cells were seeded into 96-well plates at a density of 6×10³ cells/well. After the radioresistant cells underwent a 12-hour treatment with 100% NTPAM, they were washed twice with cold PBS, and the medium was replaced with fresh medium. Then, each well was added 10 μL of Cell Counting Kit-8 solution (CCK-8, Dojindo). The culture plates were then incubated for 1–2 hours, and absorbance was measured at 450 nm using a spectrophotometer.

The NTPAM-intervened (12-hour) and untreated control RR-HNSCC cell lines were harvested, and RNA was isolated using the Qiagen RNeasy Kit (Qiagen) following the manufacturer’s protocol. All experiments were conducted under clean conditions, and the equipment was preautoclaved. Sequencing was performed by Macrogen, Inc. Briefly, the quality of extracted RNA was evaluated using an Agilent 2100 Bioanalyzer RNA Nano Chip (Agilent Technologies). RNA libraries were constructed with 1 μg of RNA using the TruSeq Stranded mRNA Sample Preparation Kit v2 (Illumina) according to the manufacturer’s protocol. Library quality was analyzed using an Agilent 2100 Bioanalyzer equipped with an Agilent DNA 1000 kit. All samples were sequenced on an Illumina HiSeq 2500 platform (Illumina), yielding an average of 38 million paired-end 100 base-pair reads. The study included four groups, with three replicates: RR-MSK-QLL1, RR-SCC15, NTP-MSK-QLL1, and NTP-RR-SCC15 (after 12-hour NTPAM intervention). Transcriptomic data were analyzed using the Deseq2 method to identify differentially expressed genes (DEGs). R (version 4.3.2) was used to generate volcano plots and differential heat maps.

Gene Set Variation Analysis (GSVA) assessed ROS and cell death pathway changes using the “GSVA” R package. The ROS pathway data were downloaded from MSigDB (https://www.gsea-msigdb.org/gsea/index.jsp). The gene set for cell death was based on relevant literature [23,24]. Radiotherapy resistance gene sets were referenced in studies [25-29]. Gene sets were downloaded from MSigDB (https://www.gsea-msigdb.org/gsea/index.jsp). Subsequently, Gene Set Enrichment Analysis (GSEA) was performed on the transcriptomic data of NTPAM treatment using the “fgsea” R package.

The process for extracting total cellular RNA is detailed in section “Transcriptome analysis.” The complementary DNA (cDNA) were synthesized with 5 μg of total RNA and TOPscript RT DryMIX (Enzynomics) according to the manufacturer’s instructions. Amplification was performed using the SYBR Green qPCR Master Mix (Thermo Fisher Scientific). The real-time quantitative polymerase chain reaction reactions were performed for 40 cycles of 95 °C for 15 seconds, 60 °C for 1 minute, and 72 °C for 1 minute. Melting curve analysis was performed. The qPCR primer sequences were as follows: human ABCC3-F: 5´-TAC GTG GAC CCA AAC AAT GTG-3´, ABCC3-R: 5´-GGA ATT GCT GGA TCC GTT TCA-3´; human DUSP16-F: 5´-GTT CTC TCG TTG TTT CCC TGG-3´, DUSP16-R: 5´-CAG CTC CTT GTT GAG GAC ATC-3´; human PDGFB-F: 5´-CCA TTC CCG AGG AGC TTT ATG-3´, PDGFB-R: 5´-GGT CAT GTT CAG GTC CAA CTC-3´; human RAF1-F: 5´-TCG ATG TCA GAC TTG TGG CTA-3´, RAF1-R: 5´-TCA AAG AAG GTA GTG CTG GGA-3´; and human THBS1-F: 5´-CAT CAG TGA GAC CGA TTT CCG-3´, THBS1-R: 5´-ATA ACC TAC AGC GAG TCC AGG-3´. We used the 2^−ΔΔCt model for relative quantification of real-time qPCR fold change. The cycle threshold values obtained by real-time qPCR were imported into Microsoft Excel. The formulas for 2^−ΔΔCt method are as follows.

ΔCt_ctrl=(Ct_ctrl−Ct_GAPDH), ΔCt_NTPAM=(Ct_NTPAM−Ct_GAPDH)

ΔΔCt_ctrl=ΔCt_ctrl−ΔCt_ctrl average

ΔΔCt_NTPAM=ΔCt_NTPAM−ΔCt_ctrl average

The relative expression values are calculated as 2^{−ΔΔCtctrl}, 2^{−ΔΔCtNTPAM}.

The Cancer Genome Atlas (TCGA)-HNSCC dataset was downloaded from the official TCGA website (https://portal.gdc.cancer.gov/). It includes transcriptomic data from 520 tumor tissues and 46 adjacent non-cancerous tissues. The TCGA clinical dataset was downloaded using the R package “TCGAbiolinks.” Based on the information from “additional radiation therapy,” “radiation therapy,” and “follow-up treatment success,” the outcomes of the patients who were treated with radiation can be determined. According to the outcomes of “complete remission/response” (CR), “partial remission/response” (PR), “stable disease” (SD), “persistent disease” (PED), and “progressive disease” (PRD), the patients were divided into two groups: “responders” (CR/PR) and “non-responders” (SD/PED/PRD) (PED patients were originally CR but developed new neoplasms, so they were grouped with SD and PRD) [30]. Survival analysis for specific selected genes was conducted using the “survival” and “survminer” R packages. The Kaplan-Meier method and log-rank test were used for survival analysis. The GSE213862 dataset includes 46 normal tissues, 43 negative-margin tissues, and 46 tumor tissues. This database was downloaded from the National Center for Biotechnology Information (NCBI) website (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi). According to the gene response to NTPAM intervention, we classified the genes into immediate-early genes (DUSP16, PDGFB, and THBS1) and delayed primary response genes (ABCC3 and RAF1) [31]. To verify the pathways associated with immediate-early genes, we categorized them into high- and low-expression groups based on the median expression levels of the genes [32]. Specifically, we performed a differential analysis between the DUSP16^low+PDGFB^high+THBS1^high and DUSP16^high+PDGFB^low+THBS1^low groups to identify the pathways associated with specific immediate-early genes. Differential analysis of TCGA-HNSCC and Gene Expression Omnibus data was performed using the R package “Deseq2.” The R package “clusterProfiler,” “ReactomePA,” “fgsea” were utilized for GSEA analysis, applying gene sets from Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/), Gene Ontology (GO; https://geneontology.org/), and Reactome (https://reactome.org/).

A flowchart of the study is shown in Fig. 2. Significant gene pathways were selected using GSEA. We then identified the genes associated with these pathways. Genes with P<0.05, |log2(fold change)|>1.5 genes underwent survival analysis in TCGA-HNSCC. Genes showing statistically significant OS and disease-specific survival (DSS) were identified. Finally, considering the clinical information and NTPAM treatment outcomes, genes consistent with the selection criteria for both risk and protective factors were chosen.

Data were analyzed using SPSS version 26 (IBM Corp.) and R software version 4.3.2 (R Foundation for Statistical Computing). Normality was assessed for the pre-analyzed data. If the data exhibited a normal distribution, comparisons between the two groups were conducted using the t-test, whereas comparisons among multiple groups were performed using analysis of variance. Pairwise comparisons of normally distributed data were conducted using Dunnett’s t-test to compare multiple experimental groups with the control group. Non-normally distributed data were analyzed using non-parametric tests such as the Wilcoxon test. The results are expressed as mean±standard deviation. Statistical significance was set at α=0.05.

We successfully established two RR-HNSCC cell lines by subjecting the parental HNSCC cells to chronic radiation stimulation. The clonogenic assay results indicated that the RR-MSK-QLL1 and RR-SCC15 cell lines exhibited higher cell viability than their respective parental cell lines following exposure to a 4 Gy radiation dose (Fig. 3A and B). To assess the impact of NTPAM on RR-HNSCC cells, we evaluated the viability of both untreated and NTPAM-treated RR-MSK-QLL1 and RR-SCC15 cells using the CCK8 assay. The treatment with NTPAM successfully inhibited the proliferation of these cell lines (Fig. 3C and D). Several studies have shown that NTPAM induces ROS-dependent cell death in various cancers. Consequently, we explored the changes in ROS and cell death signaling pathways in RR-HNSCC cell lines following NTPAM treatment through bioinformatics analysis. GSVA was employed to evaluate changes in ROS and cell death pathways based on RNA sequencing results. Using gene sets from the ROS pathway, NTPAM therapy was found to activate the ROS pathway in both RR-HNSCC cell lines (Fig. 3E and F). Post-NTPAM treatment, the RR-MSK-QLL1 cell line showed an increase in ROS-related pathways according to both the GO and Hallmark databases. In contrast, the RR-SCC15 cell line only demonstrated upregulation of ROS pathways in the GO database. The results of cell death analysis revealed that NTPAM treatment significantly activated the ferroptosis (t=10.15, P<0.001) and autophagy pathway (t=4.02, P=0.16) in RR-MSK-QLL1 cells (Fig. 3G). In contrast, NTPAM treatment cells had higher apoptosis (t=4.23, P=0.01), ferroptosis (t=5.51, P=0.005), and alkaliptosis (t=6.54, P=0.003) scores in RR-SCC15. (Fig. 3H). In RR-HNSCC cells, NTPAM treatment activated the ROS pathway and promoted cell death.

Combining 58 radiotherapy resistance pathways (TCGA) [25] with seven additional gene sets [26-28], along with a set of 21 radiotherapy-resistant genes [29], resulted in a total of 66 pathways associated with radiotherapy resistance. To investigate the regulation of pathways associated with radiotherapy resistance following interaction with NTPAM, GSEA was performed on the RNA sequencing data of RR-HNSCC cell lines that received or did not receive NTPAM treatment, focusing on 66 pathways associated with radiotherapy resistance. After NTPAM treatment of the RR-MSK-QLL1 cell line (Fig. 4A), 15 pathways related to radiotherapy resistance exhibited significant changes, whereas in the RR-SCC15 cell line, 8 pathways displayed significant results (Fig. 4B). In both cell lines, the KEGG ribosome, KEGG MAPK signaling, Reactome peptide ligand-binding receptors, and KEGG cytokine-cytokine receptor interaction pathways were upregulated after NTPAM treatment. NTPAM treatment reduced the expression of the Normalized Enrichment Score of Reactome heparan sulfate/heparin (HS-GAG) metabolism signaling pathway. The GSEA analysis revealed that 18 pathways related to radiotherapy resistance were correspondingly influenced by NTPAM treatment.

Subsequently, to identify hub genes affected by NTPAM treatment in RR-HNSCC cells, we analyzed genes that were significantly altered in pathways associated with resistance to radiation therapy. This analysis encompassed a total of 1,924 genes. DEG analysis showed that, of these genes, 290 were significantly upregulated and 184 were downregulated in the RR-MSK-QLL1 cell line (Fig. 4C and E). In the RR-SCC15 cell line, NTPAM treatment resulted in the significant upregulation of 181 genes and downregulation of 209 genes from the same pool of 1,924 genes (Fig. 4D and F). Across both RR-HNSCC cell lines, there were 117 genes commonly upregulated and 73 commonly downregulated following NTPAM therapy.

The TCGA-HNSCC database was utilized to perform survival regression analysis using 117 upregulated and 73 downregulated DEGs from both cell lines. This analysis helped identify key genes responsive to NTPAM therapy in RR-HNSCC cell lines. We obtained survival analysis results for eight genes that were consis tent with the effects of NTPAM treatment. The genes ABCC3 (hazard ratio for overall survival [HR_OS]=1.35, hazard ratio for DSS [HR_DSS]=1.67), PDGFB (HR_OS=1.88, HR_DSS=1.85), and THBS1 (HR_OS=1.67, HR_DSS=2.13) were identified as risk factors for poor OS and DSS in patients with head and neck cancer (Table 1). After NTPAM treatment, the expression of these genes was downregulated in RR-MSK-QLL1 and RR-SCC15 cells (Table 1). Similarly, the protective genes DUSP16 (HR_OS=0.69, HR_DSS=0.6), DUSP8 (HR_OS=0.76, HR_DSS=0.65), LTA (HR_OS=0.74, HR_DSS=0.67), PIK3R3 (HR_OS=0.64, HR_DSS=0.59), and RAF1 (HR_OS=0.74, HR_DSS=0.68) were upregulated in the NTPAM-treated RR-HNSCC cell lines (Table 1).

Subsequently, we evaluated the expression of eight genes in tumor and adjacent normal tissues using the TCGA-HNSCC database. The analysis revealed that the expression levels of ABCC3, PDGFB, and RAF1 were significantly different from those in normal tissues, aligning with the findings from the survival regression analysis (Fig. 5A, E, and G). Although DUSP16 (P=0.374) and THBS1 (P=0.098) did not reach statistical significance (Fig. 5B and H), the trends observed in the survival analysis were consistent with our data (Table 1). In contrast, despite the survival regression results, the expression of DUSP8, LTA, and PIK3R3 was higher in cancer tissues than in adjacent normal tissues (Fig. 5C, D, and F). Consequently, our analysis primarily focused on the expression of ABCC3, PDGFB, RAF1, DUSP16, and THBS1.

Next, we explored the relationship between these five genes and radiotherapy prognosis. In the TCGA-HNSCC cohort, 309 patients underwent radiotherapy (CR, 163; PR, 5; SD, 4; PED, 6; PRD, 57; missing data, 74). Comparing the expression of these five genes between responders (CR and PR) and non-responders (SD, PED, and PRD), high expression levels of PDGFB and THBS1 were associated with a poor prognosis in radiotherapy for head and neck cancer (Fig. 6C and E) [30]. However, no differences in ABCC3, DUSP16, and RAF1 were observed among radiotherapy patients (Fig. 6A, B and D). These findings suggested that the PDGFB and THBS1 genes, which NTPAM regulated, were also associated with clinical radioresistance and prognosis in HNSCC. In the radioresistance pathway, ABCC3 and PDGFB were enriched in the WP NRF2 pathway. DUSP16, PDGFB, RAF1, and THBS1 were enriched in the KEGG MAPK signaling pathway, KEGG melanoma, KEGG focal adhesion, and REACTOME hemostasis pathways (Fig. 6F). We subsequently derived a risk-scoring formula using LASSO Cox regression. The risk ratio between the high-risk-score group and the low-risk-score group was 2.26 (P=0.001) (Fig. 6F). Additionally, the expression of ABCC3, PDGFB, and THBS1 was positively correlated with the risk score, whereas that of DUSP16 and RAF1 was negatively correlated (Fig. 6G).

To explore whether the therapeutic effects of NTPAM on RR-HNSCC cells depended on the time of manifestation, we treated RR-HNSCC cells with NTPAM for 6 and 12 hours and measured the mRNA expression of the five genes using PCR. Within the first 6 hours, NTPAM regulated the expression of DUSP16, PDGFB, and THBS1, and this regulation continued until 12 hours. However, the expression of ABCC3 and RAF1 only changed 12 hours after NTPAM treatment (Fig. 7). Consistent with previous results, the tumor candidate protective genes DUSP16 and RAF1 showed increased expression levels following the NTPAM intervention, whereas NTPAM suppressed the potential oncogenes ABCC3, PDGFB, and THBS1.

According to the PCR results, NTPAM intervention regulated the expression of DUSP16, PDGFB, and THBS1 in the primary response, whereas it regulated the expression of ABCC3 and RAF1 in the secondary response [31]. NTPAM may regulate subsequent cellular responses via DUSP16, PDGFB, and THBS1. To explore the pathways related to immediate-early genes, we divided the samples into two groups (DUSP16^high+PDGFB^low+ THBS1^low and DUSP16^low+PDGFB^high+THBS1^high) based on their gene expression levels. These groups were then compared to predict the early regulation of signaling pathways by NTPAM. In the DUSP16^low+PDGFB^high+THBS1^high group and the DUSP16^high+ PDGFB^low+THBS1^low group of tumors, the GSEA-GO results indicated an association between upregulation of the extracellular matrix (ECM) signaling pathway and downregulation of the T-cell receptor signaling pathway (Fig. 8A). In the GSEA-KEGG analysis, upregulation of the cellular processes-focal adhesion pathway was observed, accompanied by activation of ECM, phosphoinositide 3-kinase (PI3K)/protein kinase B (AKT), and transforming growth factor-beta (TGF-β) signaling pathways, whereas some metabolism pathways, including arachidonic acid metabolism, linoleic acid metabolism, and arginine biosynthesis were downregulated (Fig. 8B). Consistent with the KEGG results, GSEA-REACTOME analysis revealed the upregulation of the ECM signaling pathway and activation of the mesenchymal-epithelial transition factor (MET) signaling pathway. Additionally, the downregulation of fatty acid metabolism was observed in these high-risk tumors (Fig. 8C). Overall, in the primary response of RR-OSCC to NTPAM, the genes DUSP16, PDGFB, and THBS1 were activated. The GSEA results indicated that the ECM, fatty acid metabolism, and PI3K-AKT, TGF-β, and MET signaling pathways may be responsible for the changes in the NTPAM immediate-early genes.