Tissue samples were obtained from the Sun Yat-sen University Cancer Center (SYSUCC; Guangzhou, China), including 9 paraffin-embedded specimens from patients who experienced tumor recurrence within 5 years after radiotherapy, as well as 9 samples from patients without recurrence (Human ethics approval code: No. B2022-313-01 Date:2022-05-18). The demographic characteristics and clinical data of the participants were retrieved from medical records (Supplementary Table S1). The study protocol was approved by the Institutional Ethical Review Board of Sun Yat-sen University Cancer Center (No. B2022-313-01), and all patients provided written informed consent for the use of their tissue samples upon admission. All procedures, including the use of human tissue specimens and analysis of clinical data, were conducted in strict accordance with the guidelines of the Declaration of Helsinki.

The untargeted metabolomics analysis was performed using ultra-high-performance liquid chromatography coupled with high-resolution mass spectrometry (UHPLC-OE-MS). The following high-purity reagents were used: methanol (LC-MS grade, CNW Technologies, CAS: 67-56-1), acetonitrile (LC-MS grade, CNW Technologies, CAS: 75-05-8), hydrochloric acid (analytical grade, Titan, CAS: 7647-01-0), sodium chloride (analytical grade, Sangon Biotech, CAS: 7647-14-5), ultrapure water (Watsons), acetic acid (LC-MS grade, Sigma-Aldrich, CAS: 64-19-7), and 2-propanol (LC-MS grade, CNW Technologies, CAS: 67-63-0).

Sample preparation and metabolite extraction were conducted using a JXFSTPRP-24 homogenizer (Shanghai Jingxin Technology Co., Ltd.) and a PS-60AL ultrasonic processor (Shenzhen Leidebang Electronics Co., Ltd.). Centrifugation was performed using a Heraeus Fresco17 refrigerated centrifuge (Thermo Fisher Scientific). The extracts were concentrated using an LGJ-10C freeze dryer (Sihuan Furuike Instrument Technology Development Co., Ltd.).

Chromatographic separation was carried out on a Vanquish UHPLC system (Thermo Fisher Scientific), and high-resolution mass spectrometric detection was performed using an Orbitrap Exploris 120 mass spectrometer (Thermo Fisher Scientific). Sample weights were accurately measured using a BSA124S-CW analytical balance (Sartorius). All procedures were conducted following standard protocols to ensure reproducibility and analytical accuracy.

The QZD solution was centrifuged at 12,000 rpm (RCF = 13,800 (× g), R = 8.6 cm) for 15 min at 4°C. 300 μL of the supernatant liquids was taken, mixed with 1,000 μL of extraction solution (MeOH:ACN:H2O, 2:2:1 (v/v/v)) which contain deuterated internal standards. The mixed solutions were vortexed for 30 s, sonicated for 5 min in 4°C water bath and incubated for 1 h at −20°C to precipitate proteins. Then the samples were centrifuged at 12,000 rpm (RCF = 13,800 (× g), R = 8.6 cm) for 15 min at 4°C. The supernatants were transferred to fresh glass vials for analysis. The quality control sample was prepared by mixing an equal aliquot of the supernatant of samples.

LC-MS/MS analyses were performed using an UHPLC system (Vanquish, Thermo Fisher Scientific) with a Phenomenex Kinetex C18 (2.1 mm × 100 mm, 2.6 μm) coupled to Orbitrap Exploris 120 mass spectrometer (Orbitrap MS, Thermo). The mobile phase A:0.01% acetic acid in water; mobile phase B:IPA:ACN (1:1,v/v). The auto-sampler temperature was 4°C, and the injection volume was 2 μL. The Orbitrap Exploris 120 mass spectrometer was used for its ability to acquire MS/MS spectra on information-dependent acquisition mode in the control of the acquisition software (Xcalibur, Thermo). In this mode, the acquisition software continuously evaluates the full scan MS spectrum. The ESI source conditions were set as following: sheath gas flow rate 50 Arb, Aux gas flow rate 15 Arb, capillary temperature 320°C, full MS resolution 60,000, MS/MS resolution 15,000, collision energy: SNCE 20/30/40, spray voltage 3.8 kV (positive) or −3.4 kV (negative).

The raw data were converted to the mzXML format using ProteoWizard and processed with an in-house program. which was developed using R and based on XCMS, for peak detection, extraction, alignment, and integration. The R package and the BiotreeDB (V3.0) were applied in metabolite identification.

The disease-related genes for NPC were retrieved from the GeneCards database using “nasopharyngeal carcinoma” as the keyword. The intersection between the drug target genes and NPC-related genes was identified using Venny 2.1, resulting in a set of potential therapeutic targets for NPC.

The relationships among the herbal medicine, active components, target genes, and NPC were integrated and visualized using Cytoscape 3.9.1. The NetworkAnalyzer tool was employed to calculate the topological properties of the network, and the final “Drug-Component-Target-Disease” network was generated after optimization.

The potential therapeutic targets for NPC were imported into the STRING database (version 11.5) with the following parameters: species set to “Homo sapiens,” minimum interaction score >0.4, and removal of disconnected nodes. The resulting PPI network was exported in TSV format and imported into Cytoscape 3.9.1 for further analysis. The topological properties of the network were calculated using the NetworkAnalyzer tool. Six topological parameters, including Betweenness Centrality (BC), Closeness Centrality (CC), Degree Centrality (DC), Eigenvector Centrality (EC), Network Centrality (NC), and Local Average Connectivity (LAC), were computed using the cytoNCA plugin. Nodes with higher parameter values were considered more critical in the network, and the core target genes were identified accordingly.

Genomic DNA was extracted from tissue specimens utilizing the Allprep DNA/RNA Kit (Qiagen). Genome-wide methylation profiling was conducted using the Infinium HumanMethylation850 BeadChips (850K array, Illumina). Data preprocessing, normalization, and β-value computation were executed with the R package minfi (v1.26.2). Gene annotations were segmented into six regions: TSS1500, TSS200, 5′UTR, 1st Exon, gene body, and 3′UTR. CpG island annotations were divided into five regions: N shelf, N shore, CpG Island, S shore, and S shelf.

Quality control measures included: (i) removal of probes with detection P-values ≥0.01 in over 5% of samples; (ii) exclusion of probes on the X or Y chromosomes; (iii) elimination of probes overlapping SNPs; and (iv) discarding probes mapping to multiple genomic loci. Post-quality control, 866,092 probes were retained for subsequent analysis.

CpG probe annotations were sourced from the ENCODE Project database. Each probe was annotated with corresponding gene, genomic region, CpG island-associated regions, and functional regions.

Differentially methylated CpG sites and regions were identified using the IMA package in R. The analysis was conducted based on both gene annotation categories (TSS1500, TSS200, 5′UTR, 1st Exon, gene body, and 3′UTR) and CpG island annotation categories (N shelf, N shore, CpG Island, S shore, and S shelf). For each region, the signal values of all probes were averaged, and differential methylation regions (DMRs) were identified using the pooled t-test method. DMRs were selected based on the following criteria: P < 0.05 and |Beta.Difference| > 0.14.

The Beta value (β), which ranges between 0 and 1, represents the methylation level at a given CpG site. Values closer to 1 indicate higher methylation, while values closer to 0 indicate lower methylation. The Beta value was calculated using the following formula: beta i = max ⁡ y i , methy , 0 max y i , methy , 0 + max y i , unmethy , 0 + 100

The Beta value density curves of the samples were generated to visualize the overall methylation distribution of all CpG sites within each sample, enabling the identification of potential differences in global methylation patterns across high-risk locally advanced NPC sample groups.

The potential therapeutic targets for NPC were subjected to Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the DAVID database. For GO analysis, the top 10 biological processes with the lowest P-values were selected and visualized using bar plots. For KEGG analysis, the top 20 significantly enriched pathways were identified and visualized using bubble plots.

A total of 15 primary NPC tumor samples were collected for single-cell RNA sequencing (scRNA-seq). Fresh biopsy samples designated for scRNA-seq were thoroughly washed with phosphate-buffered saline (PBS; Thermo Fisher Scientific) and minced into fragments of less than 1 mm³. The tissue fragments were subsequently transferred into 5 mL of Dulbecco’s Modified Eagle Medium (DMEM) supplemented with collagenase IV (1 μg/mL; Thermo Fisher Scientific) and incubated at 37°C for over 40 min, with manual agitation every 10 min to facilitate enzymatic digestion. Following digestion, the cell suspension was filtered using 5-mL round-bottom polystyrene test tubes equipped with cell strainer snap caps (BD-Falcon). The filtrate was then centrifuged at 800 rpm for 5 min, after which the resulting cell pellet was collected and resuspended in a cell preservation solution (CELLBANKER 2). Cell viability and concentration were subsequently determined via Trypan blue staining using a hemocytometer. The entire single-cell processing workflow was completed within 90 min to ensure optimal cell integrity.

The prepared single-cell suspensions were loaded onto the 10X Chromium system for single-cell RNA sequencing, utilizing either the Chromium Single Cell 3′ Library & Gel Bead Kit v2 or v3 (10X Genomics) or the 5′ Library & Gel Bead Kit (10X Genomics). Each well was designed to capture between 8,000 and 14,000 cells. Library preparation was conducted following the manufacturer’s protocols, ensuring high-quality gene expression profiling. Sequencing was performed on an Illumina HiSeq X Ten platform using a 150-bp paired-end read configuration.

For quality control, unique molecular identifiers (UMIs) and mitochondrial genes were quantified using Seurat (version 4.3.0). Cells with more than 100 UMIs and less than 15% of UMIs derived from mitochondrial genes were selected for analysis. In this study, the top 20 principal components and the top 2000 highly variable genes were selected. Cells were filtered based on nFeature_RNA, retaining those with more than 100 and fewer than 5000 features. The “ScaleData” function was employed to regress out the effects of UMI counts and the percentage of mitochondrial-derived UMIs. Subsequently, major cell clusters were identified using the “FindClusters” function in Seurat. Unbiased cell type identification was achieved through Uniform Manifold Approximation and Projection (UMAP). To characterize specific marker genes for individual cell clusters, the “FindMarkers” function in Seurat was utilized to compare the gene expression profiles of cells within a given cluster against those of all other cells in the dataset. This function employs the Wilcoxon rank-sum test to identify differentially expressed genes (DEGs) between the two groups. The resulting p-values were then adjusted for multiple testing using the Bonferroni correction, accounting for the total number of genes analyzed. Marker genes were defined as those with an adjusted p-value below 0.05 and exhibiting at least a two-fold higher average expression level in the cluster of interest compared to all other clusters. To validate the annotation of each cell cluster, we relied on the expression of canonical marker genes. Specifically, T/NK cells were identified by the expression of CD3D, CD3E, and NKG7; memory B cells were characterized by CD79A and MS4A1. Myeloid cells were identified through the expression of CD68, CD163, CD14, LYZ, CD74, CLEC9A, and CD1C. Stromal cells were distinguished by high expression of COL1A1, ACTA2, and TAGLN. On the other hand, tumor cells were identified by the expression of CCR7, LEF1, and TCP7. This comprehensive analysis, leveraging statistical tests and the expression of established marker genes, enabled accurate annotation and characterization of each cell cluster. The expression levels and distribution of genes were visualized using FeaturePlot and ggplot2.

The cell lines used in this study were authenticated and provided by Dr. Zeng from the Cancer Center of Sun Yat-sen University. The human nasopharyngeal carcinoma (NPC) cell line, HONE1, was cultured in Roswell Park Memorial Institute (RPMI)-1640 medium (Invitrogen, Carlsbad, CA, USA), supplemented with 10% fetal bovine serum (FBS; Yingjie Biotechnology).

Cell Proliferation Assay: Cell proliferation was assessed using the Cell Counting Kit-8 (CCK-8) (Dojindo Laboratories, Kumamoto, Japan). A total of 1 × 10³ cells per well were seeded into 96-well plates, and cell viability was measured every 24 h for a period of 5 days, following the manufacturer’s instructions.

Transwell Migration and Invasion Assays: For the migration and invasion assays, 5 × 10⁴ cells (for migration) or 1 × 10⁴ cells (for invasion) were seeded in the upper chamber of Transwell inserts (BD Biosciences, San Jose, CA, USA) in serum-free medium. For the invasion assay, the chamber was pre-coated with Matrigel (BD Biosciences) to mimic the extracellular matrix. The lower chamber was filled with medium containing 10% FBS as a chemoattractant. After 24 h of incubation, cells that had migrated or invaded through the membrane were fixed, stained with hematoxylin, and visualized under an inverted microscope.

Colony Formation Assay: To assess clonogenic survival, single-cell suspensions (200–10,000 cells per well) were seeded into 6-well plates and treated with QZD (10 mg/mL). After 10–14 days, once colonies had formed, the wells were washed with PBS, fixed with methanol, and stained with crystal violet. Colonies consisting of more than 50 cells were counted to evaluate the treatment effect.

Total protein was extracted from HONE1 cells using RIPA buffer (with protease and phosphatase inhibitors). After lysis, proteins were separated by SDS-PAGE and transferred to PVDF membranes. The membranes were blocked with 5% non-fat milk in TBST, then incubated overnight with primary antibodies (e.g., anti-AKT1, 1:1000, Cell Signaling Technology) at 4°C. After washing, membranes were incubated with HRP-conjugated secondary antibodies (1:2000, Thermo Fisher Scientific) for 2 h at room temperature. Proteins were detected using chemiluminescence (ECL, GE Healthcare) and imaged with the FlourChen M system. Band intensity was quantified using ImageJ software (version 1.53c).

Total RNA was extracted from cells using TRIzol reagent (Invitrogen) following the manufacturer’s protocol. Complementary DNA (cDNA) was synthesized using the HiScript III RT SuperMix (Vazyme, Nanjing, China) according to the supplier’s instructions. Quantitative real-time PCR (qPCR) was performed using the CFX96 Touch Real-Time PCR System (Bio-Rad, Hercules, CA, USA) with SYBR Green qPCR SuperMix-UDG (Invitrogen). Each sample was analyzed in triplicate, and the threshold cycle (Cq) values were recorded. Gene expression levels were normalized to glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as an internal control. The primer sequences used for qPCR analysis are listed in Table 1.

The 2D structures of the core active components were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) and converted into 3D structures using ChemBioOffice Ultra 13.0.2-Chem3D. Energy minimization was performed to optimize the structures. The 3D structures of the core target proteins were retrieved from the Protein Data Bank (PDB) (https://www.rcsb.org/) and preprocessed using PyMOL 2.5.2 to remove water molecules and irrelevant residues. The docking grid box was defined using AutoDockTools-1.5.7 to encompass the entire protein structure. Semi-flexible molecular docking was performed using AutoDock Vina (version 1.1.2), and the binding affinity (expressed as absolute binding energy) was calculated. Higher absolute binding energy values indicated stronger receptor-ligand interactions. The docking results were visualized and analyzed using PyMOL and Discovery Studio (version 4.5).

All statistical analyses, including Wilcoxon rank-sum test, and ANOVA, were carried out in R (version 4.2.0). Statistical significance was defined as p < 0.05. Data visualization, including bar charts, point charts, line graphs, and other figures, was generated using the ggplot2 package (version 3.4.2), while heatmaps were constructed with the Heatmap package (version 2.14.0).

The QZD sample was obtained from the Traditional Chinese Medicine Pharmacy of Sun Yat-sen University Cancer Prevention and Treatment Center. The QZD was prepared by decocting 1,600 mL of water until it reached a rolling boil, then simmering over low heat until the volume was reduced to 200 mL, with a total decoction time exceeding 2 h. A second extraction was performed by adding 800 mL of water to the remaining dregs, boiling again, and simmering over low heat until the volume was reduced to 100 mL. The two extracts were combined for use.

The primary components of QZD were analyzed using UHPLC-OE-MS (Supplementary Table S2). Representative ion chromatograms are shown in Figures 1a,b.

To explore the pharmacological mechanisms of QZD, we integrated target information from both its active compounds and NPC-related pathogenic genes. The overlapping targets were mapped and visualized in Cytoscape 3.10.0 to construct a “compound-target-disease” network, as shown in Figure 1c. This network consists of 356 nodes and 1408 edges, illustrating the complex multi-component, multi-target interactions between the decoction and the disease. Notably, a single active compound may interact with multiple targets, while different compounds may share the same target. Among them, the five compounds with the highest degree values (degree ranking, indicating their extensive target interactions) were MD5, SDH5, XS5, TZS3, SSB2 and MD6 (Figure 1c; Table 2; Supplementary Table S3).

To identify core therapeutic targets of QZD for NPC treatment, we integrated data from multiple databases, and identified 285 overlapping anti-NPC targets. A protein-protein interaction (PPI) network was then constructed by importing these 285 targets into the STRING database. The resulting network contained 284 nodes and 6901 edges, with the top 14 nodes ranked by degree centrality listed in Table 3. These nodes were densely interconnected, highlighting their potential role as key regulatory hubs in NPC treatment (Figures 1d,e; Table 3).

To investigate the underlying mechanisms of high-risk locally advanced NPC, we first performed differential methylation analysis on methylation sequencing data (Figure 2a). The results revealed that DMGs were significantly enriched in biological processes such as regulation of neuron projection development, regulation of nervous system development, and axonogenesis, which are known to be involved in neural development and cellular signaling pathways. These findings suggest that alterations in cellular structure and signal transduction in high-risk locally advanced tumors may contribute to increased invasiveness and disease progression.

Additionally, DMGs were highly enriched in cellular components related to synaptic membrane, apical part of the cell, cell leading edge, and neuronal cell body, indicating that modifications in these structural components may facilitate tumor cell migration and invasion. This is particularly relevant as enhanced cellular motility is a hallmark of tumor aggressiveness and progression.

At the molecular function level, DMGs were significantly associated with metal ion transmembrane transporter activity, actin binding, triphosphatase regulator activity, and GTPase regulator activity, highlighting the role of cytoskeletal dynamics and intracellular signaling in tumor cell behavior. Dysregulation of these activities may enhance tumor adaptation to the microenvironment and enable evasion of therapy-induced stress (Figure 2b).

Pathway enrichment analysis further confirmed the involvement of multiple signaling pathways, including Axon guidance, Rap1 signaling pathway, Hormone signaling, Regulation of actin cytoskeleton, Calcium signaling pathway, Cytoskeleton in muscle cells, cAMP signaling pathway, and Neuroactive ligand-receptor interaction (Figure 2c). These pathways are known to regulate cell adhesion, motility, and proliferation, all of which are crucial for tumor progression. Notably, the repeated enrichment of actin cytoskeleton regulation underscores its pivotal role in high-risk locally advanced NPC, potentially by enhancing the migration and invasiveness of residual tumor cells following treatment. We then identified key regulatory proteins by intersecting NPC- related genes, QZD treatment targets, and high-risk locally advanced NPC-associated genes, resulting in 159 key proteins. Protein-protein interaction (PPI) analysis was performed to extract core regulatory proteins, revealing that AKT1, MTOR, HIF1A, IL6, SRC, and ESR1 serve as central regulatory targets (Figures 2d,e). These core proteins are well-known mediators of tumor survival, proliferation, and metastasis, suggesting that their dysregulation may play a critical role in high-risk locally advanced NPC progression and resistance to treatment.

We then identified key regulatory proteins by intersecting NPC-related genes, QZD treatment targets, and high-risk locally advanced NPC-associated genes, resulting in 159 key proteins. Protein-protein interaction (PPI) analysis was performed to extract core regulatory proteins, revealing that AKT1, MTOR, HIF1A, IL6, SRC, and ESR1 serve as central regulatory targets. These core proteins are well-known mediators of tumor survival, proliferation, and metastasis, suggesting that their dysregulation may play a critical role in high-risk locally advanced NPC progression and resistance to treatment.

Pathway enrichment analysis of the 159 key proteins identified from the intersection of NPC-related genes, QZD treatment targets, and high-risk locally advanced NPC-associated genes revealed that the primary intervention pathways of QZD in high-risk locally advanced NPC patients include Apoptosis, PI3K-Akt signaling pathway, Th17 cell differentiation, C-type lectin receptor signaling pathway, Thyroid hormone signaling pathway, EGFR tyrosine kinase inhibitor resistance, Proteoglycans in cancer, Endocrine resistance, Prostate cancer, and Lipid and atherosclerosis (Figure 3a).

These findings suggest that QZD may exert therapeutic effects in high-risk locally advanced NPC by modulating multiple oncogenic pathways, particularly through apoptosis regulation and PI3K-Akt signaling inhibition, which are crucial for tumor survival and drug resistance. Additionally, the impact of QZD on immune-related pathways, such as Th17 cell differentiation and C-type lectin receptor signaling, indicates a potential role in modulating the tumor microenvironment to counteract disease progression and immune evasion.

To investigate the tumor immune microenvironment (TIME) in high-risk locally advanced NPC following GP chemotherapy, we analyzed single-cell RNA sequencing (scRNA-seq) data from 30 NPC patient tumor samples obtained from our hospital. Using the well-established 10X Genomics Chromium-based scRNA-seq platform, we successfully captured the transcriptomes of 179,272 single cells. These cells were classified into six major populations: Tumor cells, T & NK cells, B cells, Myeloid cells, Stromal cells, and Other (Figure 3b).

Expression analysis of key regulatory genes (AKT1, MTOR, HIF1A, SRC, and ESR1) revealed their predominant localization within the T & NK cell population, while IL6 was mainly enriched in the Myeloid compartment (Figure 3c). UMAP distribution showed that AKT1 and ESR1 were primarily expressed in tumor cells, SRC and IL6 were concentrated in the Myeloid subpopulation, while MTOR and HIF1A exhibited a broader distribution across multiple cell types (Figure 3d).

These findings suggest that QZDmay influence high-risk locally advanced NPC by modulating key signaling molecules within specific immune and tumor cell populations. The localization of AKT1 and ESR1 in tumor cells, alongside SRC and IL6 in myeloid cells, highlights potential mechanisms through which the decoction impacts tumor-immune interactions, potentially altering immune surveillance and tumor progression in high-risk locally advanced NPC.

Building on our previous pathway enrichment analysis, we identified that QZDprimarily targets apoptosis-related signaling pathways, suggesting its potential anti-tumor effects. To further validate its therapeutic role in NPC, particularly through apoptosis activation and modulation of AKT1, MTOR, HIF1A, SRC, and ESR1, we evaluated cell viability upon QZD treatment. Our results demonstrated a dose-dependent decrease in cell viability with increasing concentrations of QZD (Figure 4a). Additionally, colony formation assays revealed a significant reduction in colony counts, supporting the inhibitory effect of the decoction on NPC cell proliferation (Figure 4b).

To investigate the mechanisms of apoptosis, we examined the expression of Bax, a pro-apoptotic marker, over time. Bax levels initially decreased at 6h and 24h, likely due to drug effect, but rebounded at 48h, suggesting a sustained apoptotic response (Figure 4c). Concurrently, p-PI3K and p-AKT, two critical components of the PI3K-Akt signaling pathway, showed a significant reduction at 6h (Figure 4d), indicating early apoptosis induction.

Moreover, the expression levels of AKT1, MTOR, HIF1A, SRC, and ESR1 significantly declined over time, suggesting their involvement as key regulatory targets of QZD. Interestingly, IL6, known for its dual role as a pro-inflammatory cytokine and an anti-inflammatory myokine, exhibited a significant increase, suggesting a potential immune-modulatory effect of the decoction in the tumor microenvironment (Figures 4e–h).

The binding affinities of the previously obtained compounds with AKT1, MTOR, HIF1A, SRC, and ESR1 ranged from −14.7 to −5.6 (Figure 5). The docking results indicate that these small molecule compounds exhibit strong binding affinity with AKT1, MTOR, HIF1A, SRC, and ESR1.

Our findings demonstrate that QZD exerts its anti-tumor effects in high-risk locally advanced NPC primarily by inducing apoptosis and inhibiting the PI3K-Akt signaling pathway. QZD significantly reduces cell viability and colony formation in a dose-dependent manner while promoting apoptosis, as indicated by the Bax dynamics over time. The downregulation of p-PI3K and p-AKT confirms suppression of survival signaling, while the inhibition of key oncogenic regulators AKT1, MTOR, HIF1A, SRC, and ESR1 suggests broader effects on tumor metabolism, hypoxia adaptation, and growth factor receptor pathways. Additionally, the upregulation of IL6 may indicate an immune-mediated tumor suppression effect, warranting further investigation (Figure 6). These findings highlight QZD’s potential as a therapeutic intervention for high-risk locally advanced NPC.