Zuming Liang1,2†
Zhihao Zhou1†Jing Wang1†
Lingjun Shen3†
Yiran Li4
Litong Zhu5Gengrui Hong1Qiwen Li1
Dong Xiao1*
Xiaolin Lin6*
Taoyan Lin7*- 1Laboratory Animal Management Center, Cancer Research Institute, School of Basic Medical Sciences, Southern Medical University, Guangzhou, China
- 2Division of Life Science, State Key Laboratory for Nervous System Disorders, The Hong Kong University of Science and Technology, Hong Kong, Hong Kong SAR, China
- 3Department of Tuberculosis, Yunnan Clinical Medical Center for Infectious Diseases, The Third People’s Hospital of Kunming (The Sixth Affiliated Hospital of Dali University), Kunming, China
- 4Institute of Antibody Engineering, School of Laboratory Medicine and Biotechnology, Southern Medical University, Guangzhou, China
- 5State Key Laboratory of Organ Failure Research, National Clinical Research Center for Kidney Disease, Guangdong Provincial Institute of Nephrology, Guangdong Provincial Key Laboratory of Renal Failure Research, Division of Nephrology, Nanfang Hospital, Southern Medical University, Guangzhou, China
- 6Southern Medical University Hospital of Integrated Traditional Chinese and Western Medicine, Southern Medical University, Guangzhou, China
- 7Department of Pharmacy and Clinical Pharmacy Center, Nanfang Hospital, Southern Medical University, Guangzhou, China
Background: Nasopharyngeal carcinoma (NPC) is an aggressive malignancy with endemic prevalence in southern China. Emerging evidence highlights the critical function of the N6-methyladenosine (m6A) methylation in NPC progression, where sustained cytokine activity contributes to immunosuppression and immune evasion.
Methods: Single-cell RNA sequencing (scRNA-seq) and bulk RNA sequencing datasets were obtained from the GEO database. High-dimensional downstream analyses, including hdWGCNA, cell–cell communication, and pseudotime analysis was performed to characterize cellular interactions and transcriptional programs. Machine learning algorithms and immune infiltration profiling were integrated with MeRIP-seq to identify the key m6A-regulated gene. The post-transcriptional regulatory role was further validated in NPC cells with overexpression or knockdown of METTL14, or in cells with YWHAH silenced.
Results: B cells were identified as the primary senders to the TNF pathway, with epithelial and myeloid cells acting as influencers and receivers. YWHAH emerged as a key TNF-associated, m6A-regulated gene, with elevated expression in naive/GC B cells, interleukin (IL)-1β+ tumor-associated macrophages, and differentiated epithelial cells. METTL14-deficient increased YWHAH transcript abundance and RNA stability, whereas TNF-α stimulation further enhanced YWHAH expression. Conversely, YWHAH knockdown impaired NPC cell migration and upregulated IL−6/IL−8 expression, effects that were partially rescued by TNF−α treatment.
Conclusion: Integration of multi-omics data facilitated the identification of YWHAH as a METTL14-regulated gene, which plays a pivotal role in the NPC immune microenvironment. The elevated expression of YWHAH indicates its role in regulating immune balance. Together these findings suggest a potential regulatory link between TNF-α, YWHAH and METTL14 in the context of NPC.
1 Introduction
Nasopharyngeal carcinoma (NPC) is a head and neck malignancy strongly associated with Epstein–Barr virus (EBV) infection. It demonstrates notable ethnic and geographic disparities. The incidence rate ranges from 15 to 50 cases per 100,000 people in southern China and Southeast Asia, compared to approximately 0.4 per 100,000 in Western countries. NPC is characterized by its insidious onset and prevalence of cervical lymph node metastasis (1, 2). Constitutive, constitutive activation of distinct NF-κB complexes has been observed in nearly all EBV-positive NPC tumors (3).
Epigenetic modifications have been demonstrated to exert a pivotal role in the initiation and progression of most cancers. Among these modifications, N6-methyladenosine (m6A) has been identified as the most prevalent internal modification of mammalian messenger RNA. The dynamic regulation of m6A is orchestrated by a coordinated interplay of methyltransferases (“writers”), demethylases (“erasers”), and m6A-binding proteins (“readers”). The core writer complex, comprised of METTL3 (4), METTL14 (5), and WTAP (6), has been shown to catalyze the methylation of adenosine residues at specific consensus motifs (e.g., RRACH), thereby modulating RNA metabolism at multiple levels, including stability, splicing, nuclear export, and translation (7–10). Furthermore, the “writer” complex exhibit cooperative interactions that promote invasion, metastasis, and resistance to radiotherapy in NPC by modulating DNA damage repair pathways (11).
In our group’s previous study in the role of METTL14 in regulating epithelial–mesenchymal transition (EMT), migration, and invasion in NPC, we unexpectedly found that METTL14 exerted a marked influence on the expression of inflammation- and immunity-associated genes. The RNA sequencing data further validated that METTL14 modulates diverse inflammation-related pathways, including key components of the TNF signaling axis (e.g., TNFRSF12A, TNFRSF9, TNFSF15, TRAF1, TRAF2) and interleukin signaling members (e.g., IL7, CXCL10, NFKB2). These findings suggest a potential role for METTL14 in modulating inflammatory cascades that contribute to NPC progression (12).
NPC is distinguished by the infiltration of immune cells within and surrounding the tumor lesions, indicative of a complex tumor microenvironment (TME). This distinctive immunological landscape underscores the potential of immunotherapeutic approaches in the management of malignancy. Indeed, the blockade of the PD-1/PD-L1 immune checkpoint axis has proven to be clinically beneficial for a subset of NPC patients. However, the overall response rate to anti-PD-1 therapy remains limited, ranging from only 20% to 30% (13–15), highlighting an urgent need for deeper mechanistic insight into the TME to identify novel therapeutic targets and reliable stratification biomarkers. It is imperative to acknowledge the limitations of current genomic and transcriptomic analyses of NPC, which are predominantly constrained to bulk tissue samples with minimal cellular diversity. This limitation results in an inadequate degree of spatial and single-cell resolution, which is necessary to fully capture the heterogeneity and complexity of the TME.
Furthermore, the NF-κB signaling pathway, as a central regulatory axis of inflammatory responses, performs a pivotal role in shaping the immune microenvironment of NPC. TNF, a canonical upstream activator of the NF-κB pathway, has been shown to mediate the activation of the IKK complex through TNFR1 or TNFR2. This ultimately leads to nuclear translocation of the RelA/p50 heterodimer and subsequent transcription of target genes (16). Recent studies have indicated that NPC tumor cells may secrete TNF to stimulate surrounding myeloid cells, triggering NF-κB activation and promoting the expression of immunosuppressive mediators, thereby fostering an immune-privileged TME (17). Consistent with this, our experimental data revealed that reduced METTL14 expression was associated with elevated TNFAIP3 levels in certain NPC cell lines, suggesting activation of an intrinsic negative feedback loop aimed at dampening excessive NF-κB signaling.
The advent of single-cell RNA sequencing (scRNA-seq) has profoundly impacted the field of molecular biology, offering researchers a novel approach to delineating transcriptional profiles at single-cell resolution. This technological advancement has facilitated a more comprehensive understanding of the cellular heterogeneity present within tumors, providing researchers with a more detailed perspective on the complex biological processes underpinning cancer. Furthermore, high-dimensional weighted gene co-expression network analysis (hdWGCNA) provides an unbiased systems biology approach to identify co-regulated gene modules and hub genes (18). Utilizing this framework, we constructed a functional atlas based on scRNA-seq data, thereby enabling in-depth characterization of the phenotypic interplay between B cells and macrophages and their interactions with malignant epithelial populations within the NPC immune microenvironment.
Machine learning algorithms are increasingly employed in the analysis of scRNA-seq data to address high-dimensional feature extraction and identification of key regulatory drivers. In this study, we integrated bulk transcriptomic and single-cell data and applied both random forest (RF) and least absolute shrinkage and selection operator (LASSO) regression to identify robust feature genes potentially involved in driving immune cell heterogeneity. The incorporation of MeRIP-seq data enabled the further prioritization of feature genes that exhibited a significant response to m6A modification.
In the present analysis, we examined the dynamic differentiation trajectories of immune cell subsets within NPC and proposed YWHAH as a potential epitranscriptomic effector that contributes to tumor progression. It was observed that the m6A writer METTL14 modulates YWHAH expression. This finding was further validated through RNA stability assays, TNF-α treatment and functional assays.
2 Materials and methods
2.1 Data collection
The data presented in the study are deposited in the GEO repository, accession number are GSE162025, GSE118719 and GSE68799. GSE162025 includes 176,447 cells from 10 NPC tumor-blood pairs. GSE118719 contains 7 NPC biopsy specimens and 4 normal nasopharyngeal mucosal specimens. GSE68799 contains mRNA profiles of 42 Chinese NPC patients and 4 non-NPC tissues. Moreover, the GSE118719 dataset was designated as the training dataset, and GSE68799 was utilized as the test dataset (Supplementary Table S1).
2.2 Single-cell data quality control and batch effect
The comprehensive analysis of scRNA-seq data was performed using the “Seurat” R package (version 5). To define cellular populations with similar expression profiles and to ensure the exclusion of low-quality data resulting from cell damage or library preparation failures, we conducted quality control according to the following criteria: (1) cells expressing fewer than 200 genes were excluded, and genes detected in fewer than three cells were excluded; (2) cells with fewer than 500 or more than 2,500 expressed genes were discarded; (3) only cells with a mitochondrial gene proportion < 15%, ribosomal gene proportion > 3%, and hemoglobin gene proportion < 0.1 were retained; (4) MALAT1 housekeeping genes and mitochondrial genes were filtered out (Supplementary Figure S1A).
The gene count was normalized by the total counts for each cell, multiplied by 10,000, and then log-transformed. Among these genes, 2,000 high-variance genes were identified by the Seurat “FindVariableFeatures” function. (Supplementary Figure S1B).
Principal component analysis (PCA) was conducted to reduce dimensionality (Supplementary Figure S1C). To mitigate the impact of doublets, artefactual libraries arising from the encapsulation of multiple cells, doublets were identified and removed using the “DoubletFinder” R package (version 2.0.4; Supplementary Figure S1D). “DoubletFinder” applies an artificial nearest neighbor-based algorithm to detect cells with mixed transcriptional signatures. Doublet detection was performed independently for each sample using an expected doublet rate of 0.005, with default parameters applied otherwise. Cells retained after filtering were used for all downstream analyses (19). Finally, the “Harmony” R package (version 1.2.3) was used for batch correction.
2.3 Annotation of cell types
Cell clusters were identified using the “FindClusters” function. To annotate clusters 0-18 (dim = 25, resolution = 0.5; Supplementary Figures S1E-G), we employed the CellMarker 2.0 (http://www.bio-bigdata.center/) database in combination with the Annotation of Cell Types (ACT) tool (20, 21). This annotation was manually refined using previously established marker genes of GSE162025 (22) (Supplementary Figure S1H).
2.4 Identification of malignant tumor cells based on InferCNV analysis
To distinguish malignant tumor cells from normal cells, copy number variation (CNV) inference was performed using the “inferCNV” package on the expression matrix of distinct cellular subpopulations. Expression values were normalized to log2 (TPM + 1), where TPM was calculated as the proportion of unique molecular identifier counts per gene relative to the total UMI count per cell multiplied by 1,000,000.
CNV scores were computed using a sliding window approach with a window size of 101 genes. High-confidence normal cells from normal samples were selected as reference controls, and CNV profiles were normalized by subtracting the average reference signal. To reduce outlier effects, values greater than 1 were capped at 1, values less than −1 were capped at −1, and values between −0.3 and 0.3 were set to zero. Clusters exhibiting widespread chromosomal amplifications or deletions were classified as malignant.
2.5 High-dimensional weighted gene correlation network analysis
The hdWGCNA was employed to construct cell type-specific co-expression networks and to identify gene modules and co-expressed genes within the network (23). Modules were defined using hierarchical clustering combined with dynamic tree cutting, where distinct branches correspond to different gene modules. Hub genes were selected based on gene significance (GS) and membership degree (MM) in the degree (24).
2.6 AddModuleScore pathway scores and cellchat analysis
The “AddModuleScore” function was used to calculate the scores of functional modules for the cell cluster. This function computes the average expression of each input gene set for individual cells and contrasts it with the average expression of randomly selected gene sets with similar expression levels, yielding a relative score that reflects the activity of the corresponding functional module in each cell.
The “CellChat” R package (version 1.6.1) provides an effective analytical framework for studying intercellular interactions and communication networks (25). To investigate the potential cell–cell communications between any two different cell types in NPC, we used CellChat to calculate the interaction probabilities and statistical significance between cell pairs. The resulting interactions were visualized using chord diagrams and bubble plots, allowing for the inference of key biological interactions within the immune microenvironment of NPC.
2.7 Functional enrichment analysis
Differential expression genes analysis in the GSE118719 RNA-seq data was assessed using the “limma” (version 3.64.1), “DESeq2” (version 1.48.1) and “edgeR” (version 4.6.2) R packages to compare NPC biopsy specimens and normal nasopharyngeal mucosal specimens. Functional and mechanistic insights were further explored through Gene Ontology (GO) and KEGG pathway enrichment analyses using the “clusterProfiler” R package (version 4.16.0).
2.8 MeRIP-seq and machine learning algorithms
Methylated RNA immunoprecision sequencing (MeRIP-seq) procedures were performed according to the standard protocols of Wuhan SeqHealth Tech Co., Ltd (https://www.seqhealth.cn/sy).
Machine learning algorithms including Random Forest (RF) and least absolute shrinkage and selection operator (LASSO) regression were applied to evaluate the diagnostic potential of the putative hub genes. The GSE118719 dataset was designated as the training dataset to predict disease status and identify prognostically relevant variables. For RF analysis, we utilized the “randomForest” R package (26) alongside the “caret” R package to rank gene importance, using a threshold of importance score >2 for feature selection (27). Candidate feature genes identified from this step were subsequently input into a LASSO logistic regression model, implemented with the “glmnet” package (28), and optimal feature genes were selected based on the lambda.min parameter. Model performance was independently validated using the GSE68799 dataset as a test dataset. The predictive capacity of the optimal feature genes was assessed by receiver operating characteristic (ROC) curve analysis, and corresponding area under the ROC curve (AUC) values were calculated. Statistical significance was considered at p < 0.05.
2.9 Immune cell infiltration analysis and drug sensitivity analysis
CIBERSORT was employed to estimate immune cell compositions in the GSE118719 RNA-seq dataset (29). FPKM values were selected due to their reported advantages in deconvolution-based analyses (30). The resulting gene expression matrix, comprising 11 samples, was uploaded as the mixture file. CIBERSORT was executed using both relative and absolute quantification modes, the LM22 signature gene file, 100 permutations, and with quantile normalization disabled. While 100 permutations were applied as the minimum recommended threshold, increasing the number to 1000 did not alter the estimated absolute proportions of immune cell types, indicating robust convergence of the model (31).
Drug sensitivity analysis was performed using the oncoPredict package (version 0.2), which integrates GLDS, calcPhenotype, and IDWAS modules (32). The training dataset used for pharmacogenomic prediction (https://osf.io/c6tfx/) comprises drug response profiles and gene expression matrices from two large-scale repositories: the Cancer Therapeutics Response Portal (CTRP) and the Genomics of Drug Sensitivity in Cancer (GDSC). The “GLDS” module was used to identify key molecular features across cell lines, while “calcPhenotype” employs a ridge regression model trained on extensive transcriptomic and drug screening data to predict chemotherapeutic responses in independent datasets. Gene–drug association were further assessed using the “IDWAS” module which quantifies interactions between drug response and genomic alterations, incorporating somatic mutations or copy number variations (CNV) (33).
2.10 Developmental trajectory inference
To delineate the functional changes and potential lineage relationships among distinct immune cell populations, pseudotime analysis was performed on B cells, myeloid cells, and epithelial cells integrating CytoTRACE2 (version 1.1.0) and Monocle3 (version 1.3.7). CytoTRACE2 was applied to infer cellular developmental potential based on scRNA-seq data, while Monocle3 was used to reconstruct cell-state trajectories through dimensionality reduction, graph-based learning and pseudotime ordering.
2.11 Cell culture and reagents
A panel of NPC cell lines—including CNE2, SUNE1, HONE1-EBV, S18, 5-8F and HK1-EBV—were kindly supplied by Prof. S.-W. Tsao, Prof. Qiao Tao, Prof. Yixin Zeng, Prof. Musheng Zeng, and Dr. Dengke Li, as previously documented in the literature (34). NPC cell lines were maintained in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS; PAN Biotech, Cegrogen Biotech) at 37 °C in a humidified atmosphere containing 5% CO2. NP69 cells were cultured in keratinocyte serum-free medium (Invitrogen), and HEK293T cells were propagated in DMEM with 10% FBS (VivaCell Biotech) under identical conditions. The identity of all cell lines was verified through short tandem repeat (STR) profiling, performed by GUANGZHOU IGE BIOTECHNOLOGY Co., Ltd. (Guangzhou, China). TNF-α was purchased from Sangon Biotech (Cat. #C600021-0010).
2.12 Plasmids, lentivirus production and transduction
A 1,455 bp fragment encoding human METTL14 was amplified from the pENTER vector (Vigene Biosciences Co., Ltd., Jinan, China) and subcloned into the XbaI and BamHI restriction sites of the lentiviral backbone pCDH-EF1-MCS-GFP-Puro (System Biosciences, Cat. #CD550A-1), yielding the METTL14-overexpression construct pLV-METTL14 (pCDH-EF1-METTL14-GFP-Puro). For METTL14 silencing, short hairpin RNA (shRNA) sequences were designed and inserted into a modified pLKO.1-puro plasmid (Dahong Biosciences Co., Ltd., Guangzhou, China). All shRNA sequences are provided in Supplementary Table S2.
2.13 siRNA transient transfection
For gene silencing, specific siRNAs targeting YWHAH or negative control siRNAs that purchased from Guangzhou Tenuo Biotechnology Co., Ltd. (Guangzhou, China) were introduced into cells. Transient transfection was carried out using Lipofectamine 2000 (Thermo Fisher Scientific) according to the manufacturer’s recommendation. All siRNA sequences are provided in Supplementary Table S2.
2.14 RNA stability assay
Transcript stability was assessed by treating METTL14-knockdown or METTL14-overexpressing HONE1-EBV and 5-8F cells with actinomycin D, followed by quantification of mRNA levels over time using qRT-PCR. Half-lives were calculated by fitting the data to a one-phase exponential decay model, with normalization to the 0 h time point.
2.15 RNA extraction, library preparation, and sequencing
Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, Cat# 15596026) in accordance with the manufacturer’s protocol. Genomic DNA was removed by DNase I treatment (NEB, Cat# M0303L). RNA purity and concentration were assessed by A260/A280 ratio using a NanoDrop™ OneC spectrophotometer (Thermo Fisher Scientific) and Qubit 3.0 Fluorometer with the Qubit™ RNA Broad Range Assay Kit (Cat# Q10210). RNA integrity was evaluated with the LabChip GX Touch system (Revvity).
2.16 RNA isolation, reverse transcription and quantitative real−time PCR
Total RNA extraction (Takara, Cat# 9109), reverse transcription (Cat# R433, Vazyme Biotech Co., Ltd), and qRT-PCR (Cat# Q711, Vazyme Biotech Co., Ltd) were performed following established protocols (34–36). Primer sequences employed in the qRT-PCR assays are listed in Supplementary Table S2. GAPDH served as the internal reference gene for normalization, and relative expression levels were calculated using the 2−△△Ct method.
2.17 Western blotting assay
Western blot analysis was performed following established protocols (34–36). Membranes were probed with rabbit polyclonal antibodies targeting METTL14 (Sigma-Aldrich, HPA038002; 1:1,000) or GAPDH (Proteintech, 10494-1-AP; 1:10,000), the latter serving as a loading control. Signal detection was carried out using enhanced chemiluminescence (ECL; Millipore, WBKLS0500) according to the manufacturer’s instructions.
2.18 Transwell migration assays
Transwell assays were conducted using 24-well chambers with 8-μm pore polycarbonate membranes (Cat# 353504, Falcon®). For migration assays, 1 × 104 cells in serum-free medium with or without TNF-α (10 ng/mL) were seeded into the upper chamber (Cat# 353097, Falcon®), while the same conditions were applied to the medium in the lower chamber. Cells were allowed to migrate or invade for 14–24 hours (migration). Non-migratory cells on the upper surface were removed, and cells on the lower surface were fixed, stained with crystal violet, and counted under a light microscope across five random fields per insert.
2.19 Statistical analysis
All statistical analyses of high-throughput sequencing data were conducted in R (version 4.2.3). Differential expression between tumor and control groups was evaluated using the Wilcoxon rank-sum test (37). Correlation analyses were performed using Pearson’s correlation coefficient for normally distributed variables. P-values were adjusted for multiple testing using the Benjamini–Hochberg method.
Quantitative experimental data are presented as mean ± standard deviation (SD). Comparisons between two groups were conducted using a two-tailed Student’s t-test. One-way ANOVA followed by Tukey’s post hoc test was used for multi-group comparisons. GraphPad Prism (v9.0) was used for all statistical plots and tests. Statistical significance was defined as: P < 0.05; P < 0.01; P < 0.001; NS, not significant.
3 Results
3.1 Comprehensive transcriptomic mapping of NPC tumors using single-cell and bulk RNA-seq
The overall workflow of this study is illustrated in Figure 1. After performing quality control and filtering low-quality cells from the GSE162025 scRNA-seq data (Figures 1A, B), clusters were annotated into nine distinct cell types (Figure 1C). Genes highly associated with nasopharyngeal carcinoma (NPC) were identified using hdWGCNA (Figure 1D). Differentially expressed genes within malignant epithelial cells were further screened based on inferCNV (Figure 1E). These DEGs were then subjected to machine learning algorithms and MeRIP-seq analysis to identify a single gene most relevant to m6A modification (Figures 1F, G). Pseudotime analysis revealed distinct temporal expression pattern of the key gene across subpopulations of three cell types (Figure 1H), with experimental validation of its functional role by RNA stability assay, TNF-α treatment, and transwell migration assay with si-YWHAH (Figures 1I, J).
Figure 1. The complete workflow of this study. (A) Nasopharyngeal carcinoma research at single cell resolution. (B, C) scRNA-seq profiling. (D) hdWGCNA analysis. (E) inferCNV analysis. (F) Machine learning-based feature genes selection. (G) MeRIP–seq analysis of m6A methylation. (H) Pseudotime analysis. (I) qRT–PCR Validation of YWHAH RNA expression, RNA stability assay and TNF-α treatment. (J) Functional assays with siRNA-mediated silencing of YWHAH.
3.2 Unbiased single-cell analysis reveals transcriptional heterogeneity and immune cell diversity in NPC
A total of 176,447 cells were obtained from GSE162025, comprising matched tumor biopsies and peripheral blood mononuclear cells (PBMC) from 10 NPC patients. Violin plots show the distribution of total RNA molecules per cell (nCount_RNA), the number of detected genes (nFeature_RNA), and the proportion of mitochondrial transcripts (percent_mt) within both tumor and PBMC compartments, confirming data quality and uniformity across compartments (Figure 2A).
Figure 2. Single-cell analysis reveals the cellular composition of NPC. (A) Total transcript counts, gene counts, and mitochondrial gene percentages per sample. (B) Uniform manifold approximation and projection (UMAP) visualization of scRNA-seq clusters. (C) Identification of distinct cell populations. (D) Dot plot showing the expression of established cell-type markers. Dot color and size denote normalized gene expression and percentage of expressed cells, respectively. (E) Enrichment analysis.
Following a rigorous quality control filtration process, 158,400 cells were recovered. Unsupervised clustering using uniform manifold approximation and projection (UMAP) identified nineteen clusters and nine major cell populations: T cells (29.7%), CD8+ T cells (28.1%), NK cells (15.8%), B cells (13%), Treg cells (6.2%), myeloid cells (3.9%), NKT cells (1.6%), epithelial cells (1.2%) and plasma cells (0.5%). All identified populations were broadly represented in tumor tissues, reflecting the immune and structural heterogeneity characteristic of the NPC TME (Figures 2B, C; Supplementary Table S3). They were annotated based on their specific markers and extracted for subsequent analysis (Figure 2D).
GO and KEGG enrichment analysis further highlighted functional distinctions among cell types. Myeloid cells were enriched for chemotactic and migratory programs, epithelial cells showed enrichment in cytoskeletal organization and intermediate filament dynamics, indicative of their proliferative potential and epithelial–mesenchymal transition (EMT) activity. B cells and plasma cells were associated with B cell receptor signaling and immunoglobulin production. NKT and Treg cells were linked to immune regulatory pathways and T cell receptor (TCR) signaling, whereas CD8+ T cells exhibited enrichment in pathways related to T cell differentiation and effector responses (Figure 2E).
3.3 Identification of epithelial cell–associated transcriptional modules in NPC through inferCNV and hdWGCNA
To identify potential malignant cell populations, we applied the “inferCNV” function and defined T cell subsets as the reference group including T cells, Treg cells, NKT cells and CD8+ T cells. Epithelial, myeloid, and B cells were designated as observation groups. CNV states were inferred by comparing the average gene expression across genomic intervals. Epithelial cells exhibited pronounced copy number gains on chromosomes 6, 7, and 12, and losses on chromosomes 1, 11, and 14 (Figure 3A). Moreover, CNV scores were significantly elevated in epithelial cells compared with other cell types (Figure 3B), which strongly supports the malignant nature of the epithelial cells.
Figure 3. InferCNV analysis of cell clusters and construction of hdWGCNA. (A) Inferred copy number variation levels across 22 autosomes. (B) Comparison of overall CNV scores among different cell clusters. (C) Selection of optimal soft-thresholding power in hdWGCNA. (D) Dendrogram illustrating gene modules identified by co-expression network analysis using the selected soft-threshold. Seven distinct modules are shown. (E) kME values of the top 10 hub genes within each module. (F) UMAP visualization colored by module eigengene (ME) expression. (G) Module activity across nine cell clusters.
To further characterize transcriptional differences between epithelial cells and other cells, we performed differential expression analysis using the wilcox.test, identifying 1,764 differentially expressed genes (DEGs), including 1,393 upregulated and 371 downregulated genes (Supplementary Table S4).
Given that EBV is recognized as a central etiological factor in NPC that initially targets B cells and subsequently reshapes the TME, we applied hdWGCNA to identify co-expression modules. During the construction of an unsigned primary epithelial cell network, we determined that a soft-thresholding power (β) of 9 produced a scale-free topology fitting index of 0.80, indicating optimal network topology (Figure 3C). This analysis revealed seven distinct co-expression modules (Figure 3D).
Module eigengenes (MEs) were calculated to evaluate gene connectivity and the top 10 hub genes in each key module were ranked by intra-modular connectivity (Figure 3E). UMAP revealed that the yellow and black modules were highly expressed in both epithelial and myeloid cells, whereas the green module was predominantly activated in B and myeloid cells (Figure 3F). Consistently, bubble plot demonstrated strong positive correlations between these three modules and epithelial cell identity (Figure 3G and Supplementary Table S4).
3.4 B cell-mediated TNF may orchestrate intercellular communication in NPC
Cells were scored using the “AddModuleScore” function, and B cells were stratified into UP_Bcells and DOWN_Bcells groups based on median expression level (Supplementary Figure S1I). Cell-cell communication analysis revealed that UP_Bcells exhibited stronger and more frequent interactions with epithelial, myeloid, T, and Treg cell populations (Figure 4A). Heatmap indicated that UP_Bcells were predicted to be prominent senders of TNF signals, while receiving input signals enriched for ADGRE5, BAFF, and APRIL, compared with DOWN_Bcells (Figure 4B). Consistency, bubble plots showed that UP_Bcells displayed significantly enhanced activity in inflammatory signaling, including TNF–TNFRSF1A/1B, MIF, LGALS9–CD44/CD45, MHC-II, and CD99 (Figure 4C). Notably, increased TNF–TNFRSF1B signaling in UP_Bcells suggested potential engagement of the NF-κB signaling.
Figure 4. CellChat network analysis. (A) Network diagrams depicting cell–cell communication interactions and corresponding weighted networks across all cell types. (B) Heatmaps showing signaling pathways emitted and received by all cell types. (C) Bubble plot illustrating differentially expressed ligand–receptor pairs between UP_B cells and DOWN_B cells. (D) Heatmap depicting the functional roles of distinct cell clusters. (E) Hierarchical diagram illustrating source and target relationships in the TNF signaling pathway.
Within TNF signaling network, B cells emerged as the primary senders and influencers, whereas epithelial cells, myeloid cells, and Tregs predominantly functioned as receivers (Figures 4D, E). Collectively, this analysis suggests that a subset of activated B cells may coordinate intercellular communication in the NPC microenvironment primarily through TNF signaling, which could contribute to tumor immune modulation.
3.5 Machine learning–based analysis identifies YWHAH as a potential METTL14-regulated gene in NPC
To identify transcriptional programs closely associated with epithelial, B cells, and myeloid populations, we intersected hub genes from inferCNV-based epithelial DEGs with hdWGCNA. Specifically, 558 hub genes from the yellow, black, and green modules overlapped with 1,764 epithelial DEGs, yielding 192 potential feature genes (Figure 5A; Supplementary Table S4). They were further analyzed in the GSE118719 RNA-seq dataset, comprising 3,201 upregulated and 3,267 downregulated genes (Figure 5B).
Figure 5. Machine learning-based identification of hub genes. (A) Candidate feature genes identified by the intersection of hdWGCNA modules and DEGs. (B) Venn diagram illustrates the overlap of DEGs. (C) Key genes selected by the random forest model. (D) LASSO regression analysis for feature selection. (E) Coefficient profiles of features derived from LASSO regression. (F) ROC curve of the classifier evaluated on the independent test set GSE68799.
Machine learning algorithms were performed using GSE118719 as the training dataset. A RF classifier was first employed to model the expression profiles of 192 candidate genes, ranking their relative importance based on Gini impurity. The top 50 genes with the highest discriminative value were retained for further analysis (Figure 5C). To refine feature selection and mitigate overfitting, we subsequently applied LASSO regression. Under the optimal regularization parameter (λ), the model achieved the minimum mean squared error, ultimately identifying nine feature genes: CENPF, NUCKS1, CENPP, CCDC34, MORF4L2, LMAN1, RHEB, CDK16, and YWHAH (Figures 5D, E).
Validation in an independent test dataset GSE68799, confirmed the strong predictive capacity of the model, yielding an area under the ROC curve (AUC) of 0.9762 (Figure 5F). This result indicates excellent model performance and demonstrates the robust predictive capability of the feature genes across both datasets.
3.6 Negative feedback regulation between METTL14 and YWHAH modulates the immune landscape of NPC
MeRIP-seq analysis of SUNE1 cells with RNAi-mediated METTL14 knockdown revealed a set of genes and loci with significantly altered m6A methylation (Supplementary Table S5). Among the nine feature genes, YWHAH exhibited one of the most pronounced and consistent changes in m6A enrichment, indicating direct regulation by METTL14 in the context of NPC.
Immune infiltration analysis revealed that YWHAH expression positively correlates with Macrophages M1, CD8 T cells, and NK cells resting, all of which play a key role in innate and adaptive anti-tumor immunity. Conversely, negative correlations were observed with naïve B cells and resting CD4 memory T cells, which are typically associated with less activated or immunologically quiescent states (Figure 6A). These associations suggest a link between YWHAH expression and immune activation states.
Figure 6. Comprehensive tumor microenvironment analyses. (A) Heatmap showing the correlation between YWHAH expression and immune cell infiltration. (B) Correlation between YWHAH expression and 48 immune checkpoint genes. (C) Scatter plots depicting the associations of YWHAH with CD274 and NCR3. (D) m6A methylation levels of YWHAH in SUNE1 cells as visualized using IGV. (E) Association between YWHAH expression and IC50 values of relevant compounds. (F) Expression levels of four specific genes across distinct cell types.
In parallel, immune checkpoint analysis showed that YWHAH expression positively correlates with the immunosuppressive ligands CD70 and CD274 (PD-L1) and negatively correlates with NCR3, an activating receptor for NK cells (Figures 6B, C). These results imply that YWHAH may be involved in immune evasion in NPC.
To investigate the epitranscriptomic regulatory mechanism of YWHAH, we analyzed MeRIP-seq data from SUNE1 cells using the Integrative Genomics Viewer (IGV). Our analysis identified a significant enrichment of m6A peaks within the 3′ untranslated region (3′UTR) of YWHAH (chr22: 31,957,213–31,957,364). Notably, cells with METTL14 knockdown (LV-shMETTL14) exhibited markedly reduced methylation signals compared to controls (LV-shSCR), suggesting that METTL14 may directly mediates m6A modification of YWHAH (Figure 6D). The decrease in peak intensity suggests that METTL14-mediated methylation at this 3′UTR site may be functionally relevant, potentially affecting YWHAH mRNA stability, localization, or translational efficiency.
Furthermore, drug sensitivity analysis showed that YWHAH expression was positively correlated with predicted IC50 values for PCI-34051, JQ1, Daporinad, and Sabutoclax (Figure 6E). These associations suggest that YWHAH may contribute to drug resistance in NPC, highlighting its potential utility as a predictive biomarker for therapeutic response.
Finally, having established the functional role of YWHAH in NPC, we further investigated the GSE162025 data. Previous studies have identified a κB motif within the promoter region of YWHAH, and inhibition of NF-κB has been shown to most effectively suppress the activity of YWHAH (38). Consistent with this, violin plot (Figure 6F) revealed an inverse expression trend between METTL14 and the genes TNF, NFKB2, and YWHAH within the Epi and Myeloid clusters, suggesting a potential negative feedback regulatory relationship.
3.7 B cell subtypes diversity in NPC
Given the prominent role of B cells in TNF-mediated intercellular communication, we next sought to delineate their subpopulation diversity. A total of 20,567 B cells were identified and clustered into six subsets, including four for naïve B cells (B_C1_TCL1A, B_C3_IFITM3, B_C4_ISG15, and B_C5_HSPA1A), one for germinal center (GC) B cells (B_C6_LRMP), and one for memory B cells (B_C2_FCRL3) (Figure 7A; Supplementary Figures S2A, B). Violin plots illustrated the expression profiles of these four genes (Figure 7B), while the UMAP representation of CytoTRACE2 scores suggested that all B cell subtypes exhibited signatures of cellular maturation or terminal differentiation (Figure 7C).
Figure 7. Expression and developmental trajectories of Bcells subpopulations. (A) UMAP visualization of B cells classified into six distinct subtypes. (B) Violin plots depicting expression levels of four specific genes. (C) UMAP projection showing differentiation potency scores of B cell subtypes. (D) Heatmap of YWHAH expression across B cell subtypes. (E) Composition and differentiation trajectories of six B cells subtypes. (F) Pseudotime analysis of B cell developmental progression.
We next examined the developmental trajectory of YWHAH across B cell subpopulations. YWHAH expression was prominently elevated in the highly proliferative naïve B cell subsets B_C4_ISG15 and B_C5_HSPA1A, as well as in GC B cells (B_C6_LRMP), with expression localized predominantly at the early pseudotime locus. In contrast, YWHAH expression was diminished in the later-stage memory B cells (Figures 7D–F). Collectively, YWHAH expression exhibited limited variation among B cell subtypes, suggesting that B cells may primarily function as TNF-α signal senders rather than as key sites of YWHAH-mediated regulation (Figure 4).
3.8 Trajectory analysis suggests a transitional YWHAH+ activation state in myeloid cells
To characterize the transcriptional changes of YWHAH across distinct developmental stages of myeloid cells, we performed pseudotime analysis on ten subpopulations (Figure 8A; Supplementary Figures S2C, D). TNF, NFKB2, and YWHAH exhibited the highest expression levels in the Mac_C1_IL1B subset, whereas METTL14 showed an opposite trend (Figure 8B). CytoTRACE analysis further indicated that the differentiation potential of monocytes was markedly higher than that of dendritic cells, macrophages, and mast cells (Figures 8C, D), consistent with the canonical developmental hierarchy of the myeloid lineage.
Figure 8. Expression and developmental trajectories of myeloid cells subpopulations. (A) UMAP visualization of myeloid cells classified into ten subtypes. (B) Violin plots depicting expression levels of four specific genes. (C) UMAP projection showing differentiation potency scores of myeloid cell subtypes. (D) Boxplots illustrate developmental potential across myeloid cells subpopulations. (E) Heatmap of YWHAH expression across myeloid cells subtypes. (F) Composition and differentiation trajectories of ten myeloid cells subtypes. (G) Pseudotime analysis of myeloid cells developmental progression. (H) Heatmap presenting pathway enrichment analyses of myeloid cells subtypes.
Furthermore, DC_C3_LAMP3 and Mac_C1_IL1B occupied opposing termini of the pseudotime trajectory. Notably, in-depth analysis of macrophages revealed that within the Mac_C1_IL1B subset, YWHAH expression progressively declined toward the terminal trajectory (Figures 8E, F, G). Enrichment analyses of DC_C3_LAMP3 and Mac_C1_IL1B revealed a distinct pattern of signaling enrichment. Canonical oncogenic pathways, including Hedgehog, Wnt/β-catenin and mTORC1 signaling—were prominently enriched. Simultaneous enrichment of inflammatory and immune response pathways, such as IL6-JAK-STAT3, TNFα-NFκB signaling, interferon responses, and complement activation, reflected a highly pro-inflammatory state. Metabolic reprogramming was also evident, characterized by upregulation of glycolysis, cholesterol homeostasis, heme metabolism, and hypoxia-related programs, suggesting adaptation to the nutrient- and oxygen-deprived tumor milieu Figure 8H.
YWHAH, a member of the 14-3–3 family, showed dynamic expression in the interleukin (IL)-1β+ tumor-associated macrophages (TAMs), implicating a role in modulating TNF receptor signaling and macrophage polarization, similar to other 14-3–3 isoforms. This trajectory, likely driven by chronic exposure to TGF-β, metabolic stress and hypoxia, highlights the plasticity of macrophages in NPC.
3.9 Epithelial trajectory suggests YWHAH associated tumor progression
To further explore the expression dynamics of YWHAH within epithelial cells, the epithelial compartment was reclassified into distinct subpopulations. Given that EPCAM is a marker gene associated with stem-like, developmentally immature epithelial states, cells were segregated into Epi_C1 and Epi_C2_EPCAM subsets (Figure 9A; Supplementary Figures S2E-F), a classification that was supported by CytoTRACE analysis (Figure 9B). Violin plots revealed that TNF, NFKB2, and YWHAH were markedly upregulated in the Epi_C1 subset (Figure 9C).
Figure 9. Expression and developmental trajectories of epithelial cells subpopulations. (A) UMAP visualization of epithelial cells classified into two subtypes. (B) UMAP projection showing differentiation potency scores of epithelial cell subtypes. (C) Violin plots depicting expression levels of four specific genes. (D) Heatmap of YWHAH expression across epithelial cells subsets. (E) Composition and differentiation trajectories of two epithelial cells subtypes. (F) Pseudotime analysis of epithelial cells developmental progression.
Subsequently, a pseudotime analysis was conducted, and, in accordance with biological relevance, designated Epi_C2_EPCAM was designated as the root of differentiation. This analysis revealed one major branching point leading to three distinct lineages, with increasing cellular maturation observed as cells progressed away from the origin (Figures 9D, E, F).
Finally, pathway heatmaps revealed that Epi_C1 cells were enriched in multiple cancer-associated pathways, including those regulating cell cycle progression (G2M checkpoint, E2F targets, mitotic spindle) and genome stability (DNA repair). It is noteworthy that the upregulation of the unfolded protein response and protein secretion pathways further indicates adaptation to heightened biosynthetic stress, a hallmark of malignant epithelial cells (Supplementary Figure S2G). Collectively, these findings lend support to a state of functional reprogramming, with YWHAH acting as a potential regulatory node that facilitates tumor progression, metabolic remodeling, and immune escape within the NPC microenvironment.
3.10 Loss of METTL14 promotes TNF-α-associated YWHAH expression by enhancing RNA stability
To investigate whether METTL14 regulates YWHAH expression via m6A modification, stable METTL14-overexpressing and -knockdown cell models were first established across multiple NPC cell lines. Western blot analysis confirmed a significant increase in METTL14 protein levels in the overexpression model and effective suppression in the knockdown model (Figure 10A; Supplementary Figure S3A). Subsequently, we observed consistent upregulation of YWHAH protein expression in three METTL14-knockdown NPC cell lines (CNE2, HONE-EBV, and SUNE1). Notably, in SUNE1 cells, the expression of NFKB2, a key component of the NF-κB pathway, was also synchronously elevated, suggesting that YWHAH may be involved in linking TNF-α signaling to NF-κB pathway activation (Figures 10B, C).
Figure 10. Validation of the effect of METTL14 on YWHAH expression and function. (A) METTL14 expression was analyzed by Western blot in NPC cells which were transduced with lentivirus expressing METTL14 or shMETTL14. (B) qRT-PCR measurement of METTL14, NFKB2 and YWHAH expression in CNE2 shMETTL14-expressing cells. (C) qRT-PCR measurement of METTL14, NFKB2 and YWHAH expression in HONE1-EBV and SUNE1 shMETTL14-expressing cells. (D) shSCR and shMETTL14 HONE1-EBV were treated with Act.D (5 μg/mL) for the indicated times. (E) CON and METTL14-WT/-R298P 5-8F were treated with Act.D (5 μg/mL) for the indicated times. (F) METTL14 overexpression in SUNE1 and 5-8F counteracted the stimulatory effect of TNF-α on YWHAH expression. (G) YWHAH regulates the expression of TNF-α–related genes. (H) YWHAH silencing reduces NPC cell migration and is partially rescued by TNF-α.
To further elucidate the post-transcriptional regulatory mechanism of METTL14 on YWHAH, we performed actinomycin D-mediated mRNA stability assays. The results showed that the half-life of YWHAH mRNA was significantly prolonged in METTL14-knockdown HONE-EBV cells, indicating enhanced stability. Conversely, in 5-8F cells overexpressing wild-type METTL14, YWHAH mRNA stability was markedly reduced. In contrast, the catalytically inactive mutant METTL14-R298P lost this regulatory capacity, even resulting in YWHAH stability exceeding that of the control group (Figures 10D, E). These findings collectively suggest that METTL14 negatively regulates YWHAH mRNA stability in an m6A methyltransferase activity-dependent manner. Next, we treated METTL14-overexpressing SUNE1 and 5-8F cells with recombinant human TNF-α. TNF-α stimulation induced a time-dependent upregulation of YWHAH expression; however, sustained METTL14 overexpression significantly attenuated this inductive effect, consistent with the observed decrease in mRNA stability. Importantly, the METTL14-R298P mutant effectively restored TNF-α-induced YWHAH expression (Figure 10F), further confirming the central role of METTL14 catalytic activity in this regulatory function.
To validate the role of YWHAH in NPC cell migration and its association with TNF signaling, we performed YWHAH knockdown experiments. The results showed that silencing YWHAH significantly reduced its mRNA levels and inhibited cell migratory capacity, while TNF-α treatment partially rescued the downregulation of YWHAH expression and the impaired migration caused by YWHAH depletion (Figures 10G, H; Supplementary Figure S3B). Concurrently, the expression levels of the downstream inflammatory cytokines IL-6 and IL-8 were correspondingly altered. Given that IL-6 and IL-8 are well-established transcriptional targets of the NF-κB signaling pathway, their expression changes consistently reflect the modulation of NF-κB activity. Taking together, these results demonstrate that YWHAH contributes to the regulation of NPC cell migration, with METTL14/m6A acting as a key downstream effector of TNF-α and YWHAH regulation.
4 Discussion
Rapid advances in m6A RNA methylation studies have demonstrated its impact on nearly every aspect of RNA metabolism, including RNA expression, splicing, nuclear export, translation, decay, and RNA–protein interactions (39–41). More than one research have highlighted the critical role of m6A modifications in diverse physiological processes such as DNA damage response (42), pluripotency, and cellular reprogramming (43). However, the expression and regulatory mechanisms of the three types of m6A regulators across different tissues remain poorly understood.
Emerging evidence indicates that m6A modifications can suppress tumorigenesis in various cancers. For instance, decreased m6A levels on ADAM19 in glioblastoma (GBM) enhance its expression, promoting glioma stem cell growth and self-renewal, thereby facilitating tumor development (44). Approximately 70% of endometrial tumors exhibit reduced m6A leading to AKT pathway activation, which increases cell proliferation and tumorigenicity (45). In hepatocellular carcinoma (HCC), diminished m6A on pri-miR126 disrupts its maturation and tumor-suppressive function, accelerating disease progression (46). Moreover, FTO-mediated downregulation of m6A on ASB2 and RARA reduces these tumor suppressors’ expression, contributing to leukemogenesis (47).
Chronic inflammation is widely recognized as a key contributor to cancer initiation and progression. As a prototypical pro-inflammatory cytokine, TNF orchestrates the remodeling of the tumor microenvironment in a manner that facilitates oncogenic initiation, progression, and metastasis. Notably, all three classes of m6A regulators showcase altered expression in response to TNF-α stimulation, ultimately impacting disease development. For example, TNF-α reduces FTO expression in mesenchymal stromal cells (MSCs), shortening Nanog mRNA half-life and impairing MSC differentiation into sweat gland cells (48). During endothelial inflammation and atherosclerotic plaque formation, METTL14 enhances FOXO1 transcription through m6A modifications at VCAM-1 and ICAM-1 promoters, promoting translation through YTHDF3 and contributing to pathogenesis. Conversely, METTL14 is upregulated in diabetic nephropathy, where its overexpression increases serum TNF-α secretion and induces glomerular endothelial apoptosis by downregulating α-klotho in an m6A-dependent manner (49).
In this study, integrative single-cell analyses highlight TNF signaling as a central feature and suggest that B cells may serve as key coordinators of TNF-mediated intercellular communication, contributing to immune remodeling in NPC. By combining epitranscriptomic profiling with immune analyses, we further identify YWHAH as a critical node within this regulatory network. MeRIP-seq data indicate that YWHAH is subject to METTL14-dependent m6A modification, linking inflammatory signaling to post-transcriptional regulation. Consistent with this, YWHAH expression is associated with distinct immune states, characterized by increased M1 macrophage infiltration and reduced naïve B cells and resting CD4+ memory T cells. Moreover, the complex correlations between YWHAH and immune checkpoint molecules suggest a dual regulatory role, in which YWHAH may simultaneously promote T cell activation through co-stimulatory pathways while engaging compensatory immune inhibitory mechanisms, including PD-L1 upregulation. Similar findings in colorectal cancer show YWHAH enhances ac4C modification stability, driving CD8+ T cell exhaustion and immune escape, closely linked with immune checkpoint pathways (50). Drug sensitivity analysis indicated a significant positive correlation (p < 0.05) between the epigenetic regulator JQ1 and YWHAH expression. JQ1 inhibits BET family proteins, particularly BRD4, promoting nuclear export of METTL14 and functionally mimicking its knockdown effect, suggesting indirect regulation of YWHAH expression via m6A-related pathways (51).
Pseudotime analysis suggested that high YWHAH expression in naïve and GC B cells during early developmental stages, indicating a role in early B cell development and germinal center selection. Given that germinal center reactions rely on NF-κB/PI3K/AKT signaling, elevated YWHAH expression may contribute to the regulation of somatic hypermutation and affinity maturation (52). As a canonical 14-3–3 family member, YWHAH may support GC B cell proliferation and antigen presentation.
Within the myeloid lineage, macrophage-derived inflammatory mediators are essential for tissue protection and repair, yet cytokines such as IL-1β can cooperate with oncogenic programs to drive the transformation of mutant clones (53, 54). LAMP3+ dendritic cells promote CD8+ T-cell exhaustion via the NECTIN2–TIGIT (55) axis and recruit regulatory T cells, thereby establishing immune suppression and implicating YWHAH in shaping dendritic cell function. IL-1β+ tumor-associated macrophages (TAMs) display high expression of inflammatory and repair-associated genes but limited capacity to stimulate cytotoxic immunity, identifying them as key drivers of cancer-promoting inflammation. In human PDAC, IL-1β+ TAMs cluster near tumor subpopulations enriched for IL-1–responsive programs and aggressive traits, including oncogenic signaling, hypoxia and epithelial-to-mesenchymal transition (56). A similar spatial pattern is observed in renal cancer, where IL-1β+ TAMs accumulate at invasive margins coinciding with EMT activity (57). These findings indicate that IL-1β+ TAMs propagate inflammation, immune evasion and malignant progression. In the Mac_C1_IL1B subset, IL-1β+ macrophages assume dual roles—sustaining inflammation yet vulnerable to tumor hijacking through PGE2/TNF and IL-1β-driven feedback loops that enforce tolerance and exhaustion (58). Notably, early YWHAH upregulation may integrate NF-κB and metabolic inputs to initiate immune activation, whereas its decline in later stages reflects suppression, suggesting YWHAH may act as a nodal regulator balancing pro-inflammatory activity and immune dysfunction.
In epithelial cells, YWHAH expression markedly increases from EBV+ undifferentiated EPCAM+ progenitors to mature Epi_C1 cells at pseudotime trajectory termini, indicating a role in epithelial maturation and malignant transformation. Epi_C1 cells are enriched for G2/M checkpoint regulation, DNA repair, unfolded protein response, and elevated protein synthesis pathways, reflecting heightened proliferation and translational stress. Under these conditions, YWHAH likely facilitates malignant phenotype stabilization. EBV infection via LMP1 activates NF-κB (59) and mTORC1 (60) signaling, promoting epithelial phenotype shifts and tumorigenesis while enhancing immune evasion.
In the context of “cold tumors” such as NPC, which are characterized by low immune infiltration and poor responsiveness to immune checkpoint inhibitors, YWHAH may represents a potential regulator of tumor immune remodeling. In this study we demonstrate that METTL14 negatively regulates YWHAH expression through m6A-dependent modification of its 3’UTR, thereby promoting mRNA decay. This epitranscriptomic control links TNF-associated inflammatory signaling with post-transcriptional regulation and may contribute to both immunological and metabolic reprogramming within the NPC microenvironment.
A limitation of this study is the lack of in-depth investigation into m6A-binding proteins, including YTHDF1, YTHDF 2, YTHDF 3and YTHDC1/2 (61). m6A modifications exert their regulatory functions through interactions with distinct reader proteins, which determine RNA fate by modulating translation, degradation and nuclear export (62, 63). Previous studies have shown that YTHDF1 promotes translation by binding m6A modification site within the 3′ UTR of target mRNAs and recruiting eukaryotic initiation factor eIF3, thereby enabling cap-independent translation (64). In neuroblastoma (NB), METTL14-mediated methylated YWHAH transcripts, particularly within the 5′ UTR, were also specifically recognized by the YTHDF1, enhancing translational efficiency and promoting NB cell activity (65). In addition, YTHDF1 has been implicated in head and neck squamous cell carcinomas (HNSCCs) by regulating iron metabolism through m6A-dependent control of TFRC mRNA translation via interactions with both 3′ and 5′ UTR (66). Therefore, the interplay between METTL14 and YTH domain–containing reader proteins at the 3′ UTR of YWHAH warrants further investigation.
Collectively, the differential expression of YWHAH between immune and tumor cell compartments suggests a potential strategy for immune-selective intervention. In light of our observations that METTL14 loss activates the TNF−α pathway via YWHAH, therapeutic targeting of this regulatory axis may hold promising potential. Future studies employing epigenetic modulators such as JQ1 will further clarify the mechanistic and therapeutic implications.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.
Author contributions
ZL: Writing – original draft, Writing – review & editing. ZZ: Writing – review & editing. JW: Writing – review & editing. LS: Writing – review & editing. YL: Writing – review & editing. LZ: Writing – review & editing. GH: Writing – review & editing. QL: Writing – review & editing. DX: Writing – review & editing. XL: Writing – review & editing. TL: Writing – review & editing.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This work was supported by the National Natural Science Foundation of China (Grant No. 82203251, to TY Lin; Grant No. 82573835, 82173299 and 81872209, to D Xiao; Grant No. 81870602, to XL Lin), the Guangzhou Basic and Applied Basic Research Foundation (Grant No. 202201010909, to TY Lin), the Science and Technology Program of Guangzhou(Grant No. 2025A04J5133, to TY Lin), the Natural Science Foundation of Guangdong Province of China (Grant No. 2024A1515013024 and 2022A1515012477 to D Xiao; Grant No. 2021A1515111127, to TY Lin; Grant No. 2023A1515011844, to XL Lin), Youth S&T Talent Support Programme of Guangdong Provincial Association for Science and Technology(Grant No. SKXRC2025139, to TY Lin) and the Basic Research Foundation of Yunnan Province Local Universities (Grant No. 202401BA070001-063, to LJ Shen).
Acknowledgments
We are grateful to Prof. Andrea Ventura (Memorial Sloan Kettering Cancer Center) and Prof. Jeng-Shin Lee (Harvard Gene Therapy Initiative, Harvard Medical School) for their support in providing plasmids. We would like to acknowledge the GEO database for providing data (GSE162025, GSE118719 and GSE68799). Additionally, we extend our thanks to all members of Dr. Dong Xiao’s lab for their insightful discussions and helpful suggestions.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2026.1717039/full#supplementary-material
References
1. Chen Y-P, Chan ATC, Le Q-T, Blanchard P, Sun Y, and Ma J. Nasopharyngeal carcinoma. Lancet. (2019) 394:64–80. doi: 10.1016/S0140-6736(19)30956-0
2. Chua MLK, Wee JTS, Hui EP, and Chan ATC. Nasopharyngeal carcinoma. Lancet. (2016) 387:1012–24. doi: 10.1016/S0140-6736(15)00055-0
3. Chung GT-Y, Lou WP-K, Chow C, To K-F, Choy K-W, Leung AW-C, et al. Constitutive activation of distinct NF-kappaB signals in EBV-associated nasopharyngeal carcinoma. J Pathol. (2013) 231:311–22. doi: 10.1002/path.4239
4. Sendinc E and Shi Y. RNA m6A methylation across the transcriptome. Mol Cell. (2023) 83:428–41. doi: 10.1016/j.molcel.2023.01.006
5. Guan Q, Lin H, Miao L, Guo H, Chen Y, Zhuo Z, et al. Functions, mechanisms, and therapeutic implications of METTL14 in human cancer. J Hematol Oncol. (2022) 15:13. doi: 10.1186/s13045-022-01231-5
6. Zeng C, Huang W, Li Y, and Weng H. Roles of METTL3 in cancer: mechanisms and therapeutic targeting. J Hematol Oncol. (2020) 13:117. doi: 10.1186/s13045-020-00951-w
7. Wen T, Li T, Xu Y, Zhang Y, Pan H, and Wang Y. The role of m6A epigenetic modifications in tumor coding and non-coding RNA processing. Cell Communication Signaling. (2023) 21:355. doi: 10.1186/s12964-023-01385-w
8. Shi T, Zhang H, and Chen Y. The m6A revolution: transforming tumor immunity and enhancing immunotherapy outcomes. Cell Biosci. (2025) 15:27. doi: 10.1186/s13578-025-01368-z
9. Sun T, Wu R, and Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother. (2019) 112:108613. doi: 10.1016/j.biopha.2019.108613
10. Pan Y, Ma P, Liu Y, Li W, and Shu Y. Multiple functions of m6A RNA methylation in cancer. J Hematol Oncol. (2018) 11:48. doi: 10.1186/s13045-018-0590-8
11. Wang YM, Peng ZY, Zhang LY, Wang YX, Fan RH, Zhang H, et al. N6-methyladenosine RNA modification landscape in the occurrence and recurrence of nasopharyngeal carcinoma. World J Oncol. (2022) 13:205–15. doi: 10.14740/wjon1491
12. Zhou Z, Wang J, Shen L, Han L, Li Q, Wu A, et al. METTL14 regulates the expression of genes related to interferon, interleukin and MHC class I in nasopharyngeal carcinoma cells. Cancer Med. (2025) 14:e71371. doi: 10.1002/cam4.71371
13. Hsu C, Lee S-H, Ejadi S, Even C, Cohen RB, Le Tourneau C, et al. Safety and antitumor activity of pembrolizumab in patients with programmed death-ligand 1-positive nasopharyngeal carcinoma: results of the KEYNOTE-028 study. J Clin Oncol. (2017) 35:4050–6. doi: 10.1200/JCO.2017.73.3675
14. Ma BBY, Lim W-T, Goh B-C, Hui EP, Lo K-W, Pettinger A, et al. Antitumor activity of nivolumab in recurrent and metastatic nasopharyngeal carcinoma: an international, multicenter study of the mayo clinic phase 2 consortium (NCI-9742). J Clin Oncol. (2018) 36:1412–8. doi: 10.1200/JCO.2017.77.0388
15. Fang W, Yang Y, Ma Y, Hong S, Lin L, He X, et al. Camrelizumab (SHR-1210) alone or in combination with gemcitabine plus cisplatin for nasopharyngeal carcinoma: results from two single-arm, phase 1 trials. Lancet Oncol. (2018) 19:1338–50. doi: 10.1016/S1470-2045(18)30495-9
16. Chen G and Goeddel DV. TNF-R1 signaling: A beautiful pathway. Science. (2002) 296:1634–5. doi: 10.1126/science.1071924
17. Hayden MS and Ghosh S. Regulation of NF-kappaB by TNF family cytokines. Semin Immunol. (2014) 26:253–66. doi: 10.1016/j.smim.2014.05.004
18. Langfelder P and Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. (2008) 9:559. doi: 10.1186/1471-2105-9-559
19. Su M, Pan T, Chen Q-Z, Zhou W-W, Gong Y, Xu G, et al. Data analysis guidelines for single-cell RNA-seq in biomedical studies and clinical applications. Mil Med Res. (2022) 9:68. doi: 10.1186/s40779-022-00434-8
20. Quan F, Liang X, Cheng M, Yang H, Liu K, He S, et al. Annotation of cell types (ACT): a convenient web server for cell type annotation. Genome Med. (2023) 15:91. doi: 10.1186/s13073-023-01249-5
21. Hu C, Li T, Xu Y, Zhang X, Li F, Bai J, et al. CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data. Nucleic Acids Res. (2023) 51:D870–D6. doi: 10.1093/nar/gkac947
22. Liu Y, He S, Wang X-L, Peng W, Chen Q-Y, Chi D-M, et al. Tumor heterogeneity and intercellular networks of nasopharyngeal carcinoma at single cell resolution. Nat Commun. (2021) 12:741. doi: 10.1038/s41467-021-21043-4
23. Han J, Zhou Y, Zhang C, Feng J, Wang J, Guo K, et al. Intratumoral immune heterogeneity of prostate cancer characterized by typing and hub genes. J Cell Mol Med. (2022) 27:101–12. doi: 10.1111/jcmm.17641
24. Li R, Zhao M, Miao C, Shi X, and Lu J. Identification and validation of key biomarkers associated with macrophages in nonalcoholic fatty liver disease based on hdWGCNA and machine learning. Aging (Albany NY). (2023) 15:15451–72. doi: 10.18632/aging.205374
25. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan C-H, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9
26. Wang H, Yang F, and Luo Z. An experimental study of the intrinsic stability of random forest variable importance measures. BMC Bioinf. (2016) 17:60. doi: 10.1186/s12859-016-0900-5
27. Tian Y, Yang J, Lan M, and Zou T. Construction and analysis of a joint diagnosis model of random forest and artificial neural network for heart failure. Aging. (2020) 12:26221–35. doi: 10.18632/aging.202405
28. Friedman J, Hastie T, and Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Software. (2010) 33:1–22. doi: 10.18637/jss.v033.i01
29. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337
30. Jin H, Wan Y-W, and Liu Z. Comprehensive evaluation of RNA-seq quantification methods for linearity. BMC Bioinf. (2017) 18:117. doi: 10.1186/s12859-017-1526-y
31. Craven KE, Gökmen-Polar Y, and Badve SS. CIBERSORT analysis of TCGA and METABRIC identifies subgroups with better outcomes in triple negative breast cancer. Sci Rep. (2021) 11:4691. doi: 10.1038/s41598-021-83913-7
32. Maeser D, Gruener RF, and Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Briefings Bioinf. (2021) 22:bbab260. doi: 10.1093/bib/bbab260
33. Xu W, Zhang W, Zhao D, Wang Q, Zhang M, Li Q, et al. Unveiling the role of regulatory T cells in the tumor microenvironment of pancreatic cancer through single-cell transcriptomics and in vitro experiments. Front Immunol. (2023) 14:1242909. doi: 10.3389/fimmu.2023.1242909
34. Lin T-Y, Jia J-S, Luo W-R, Lin X-L, Xiao S-J, Yang J, et al. ThermomiR-377-3p-induced suppression of Cirbp expression is required for effective elimination of cancer cells and cancer stem-like cells by hyperthermia. J Exp Clin Cancer Res. (2024) 43:62. doi: 10.1186/s13046-024-02983-3
35. Qin Y-J, Lin T-Y, Lin X-L, Liu Y, Zhao W-T, Li X-Y, et al. Loss of PDK4 expression promotes proliferation, tumorigenicity, motility and invasion of hepatocellular carcinoma cells. J Cancer. (2020) 11:4397–405. doi: 10.7150/jca.43459
36. Li J, Lin T-Y, Chen L, Liu Y, Dian M-J, Hao W-C, et al. miR-19 regulates the expression of interferon-induced genes and MHC class I genes in human cancer cells. Int J Med Sci. (2020) 17:953–64. doi: 10.7150/ijms.44377
37. Jin W, Liu Z, Zhang Y, Che Z, and Gao M. The effect of individual musculoskeletal conditions on depression: updated insights from an irish longitudinal study on aging. Front Med. (2021) 8:697649. doi: 10.3389/fmed.2021.697649
38. Qiu Y, Dai Y, Zhang C, Yang Y, Jin M, Shan W, et al. Arsenic trioxide reverses the chemoresistance in hepatocellular carcinoma: a targeted intervention of 14–3-3η/NF-κB feedback loop. J Exp Clin Cancer Res. (2018) 37:321. doi: 10.1186/s13046-018-1005-y
39. He J, Zhou M, Yin J, Wan J, Chu J, Jia J, et al. METTL3 restrains papillary thyroid cancer progression via m6A/c-Rel/IL-8-mediated neutrophil infiltration. Mol Ther. (2021) 29:1821–37. doi: 10.1016/j.ymthe.2021.01.019
40. Wang Z, Shang J, Qiu Y, Cheng H, Tao M, Xie E, et al. Suppression of the METTL3-m6A-integrin β1 axis by extracellular acidification impairs T cell infiltration and antitumor activity. Cell Rep. (2024) 43:113796. doi: 10.1016/j.celrep.2024.113796
41. Xi Q, Yang G, He X, Zhuang H, Li L, Lin B, et al. M6A-mediated upregulation of lncRNA TUG1 in liver cancer cells regulates the antitumor response of CD8+ T cells and phagocytosis of macrophages. Adv Sci. (2024) 11:e2400695. doi: 10.1002/advs.202400695
42. Mu X, Zhao Q, Chen W, Zhao Y, Yan Q, Peng R, et al. IL-37 confers anti-tumor activity by regulation of m6A methylation. Front Oncol. (2021) 10:526866. doi: 10.3389/fonc.2020.526866
43. Johnson DE, O’Keefe RA, and Grandis JR. Targeting the IL-6/JAK/STAT3 signaling axis in cancer. Nat Rev Clin Oncol. (2018) 15:234–48. doi: 10.1038/nrclinonc.2018.8
44. Chang Q, Bournazou E, Sansone P, Berishaj M, Gao SP, Daly L, et al. The IL-6/JAK/stat3 feed-forward loop drives tumorigenesis and metastasis. Neoplasia. (2013) 15:848–IN45. doi: 10.1593/neo.13706
45. Hung C-H, Wu S-Y, Yao C-ID, Yeh H-H, Lin C-C, Chu C-Y, et al. Defective N-glycosylation of IL6 induces metastasis and tyrosine kinase inhibitor resistance in lung cancer. Nat Commun. (2024) 15:7885. doi: 10.1038/s41467-024-51831-7
46. Orange ST, Leslie J, Ross M, Mann DA, and Wackerhage H. The exercise IL-6 enigma in cancer. Trends Endocrinol Metab. (2023) 34:749–63. doi: 10.1016/j.tem.2023.08.001
47. Chan L-C, Li C-W, Xia W, Hsu J-M, Lee H-H, Cha J-H, et al. IL-6/JAK1 pathway drives PD-L1 Y112 phosphorylation to promote cancer immune evasion. J Clin Invest. (2019) 129:3324–38. doi: 10.1172/jci126022
48. Kerur N, Veettil Mohanan V, Sharma-Walia N, Bottero V, Sadagopan S, Otageri P, et al. IFI16 acts as a nuclear pathogen sensor to induce the inflammasome in response to kaposi sarcoma-associated herpesvirus infection. Cell Host Microbe. (2011) 9:363–75. doi: 10.1016/j.chom.2011.04.008
49. Cai H, Yan L, Liu N, Xu M, and Cai H. IFI16 promotes cervical cancer progression by upregulating PD-L1 in immunomicroenvironment through STING-TBK1-NF-kB pathway. Biomed Pharmacother. (2020) 123:109790. doi: 10.1016/j.biopha.2019.109790
50. Li Q, Yuan Z, Wang Y, Zhai P, Wang J, Zhang C, et al. Unveiling YWHAH: A potential therapeutic target for overcoming CD8+ T cell exhaustion in colorectal cancer. Int Immunopharmacol. (2024) 135:112317. doi: 10.1016/j.intimp.2024.112317
51. Lu X, Peng L, Ding J, Li Y, Li Q, Mengchen R, et al. A deregulated m 6 A writer complex axis driven by BRD4 confers an epitranscriptomic vulnerability in combined DNA repair–targeted therapy. Proc Natl Acad Sci. (2023) 120:e2304534120. doi: 10.1073/pnas.2304534120
52. Guldenpfennig C, Teixeiro E, and Daniels M. NF-kB’s contribution to B cell fate decisions. Front Immunol. (2023) 14:1214095. doi: 10.3389/fimmu.2023.1214095
53. Yuan S, Almagro J, and Fuchs E. Beyond genetics: driving cancer with the tumor microenvironment behind the wheel. Nat Rev Cancer. (2024) 24:274–86. doi: 10.1038/s41568-023-00660-9
54. Weeden CE, Hill W, Lim EL, Gronroos E, and Swanton C. Impact of risk factors on early cancer evolution. Cell. (2023) 186:1541–63. doi: 10.1016/j.cell.2023.03.013
55. Wang Z, Ji X, Zhang Y, Yang F, Su H, Zhang H, et al. Interactions between LAMP3+ dendritic cells and T-cell subpopulations promote immune evasion in papillary thyroid carcinoma. J ImmunoTher Cancer. (2024) 12:e008983. doi: 10.1136/jitc-2024-008983
56. Caronni N, La Terza F, Vittoria FM, Barbiera G, Mezzanzanica L, Cuzzola V, et al. IL-1beta(+) macrophages fuel pathogenic inflammation in pancreatic cancer. Nature. (2023) 623:415–22. doi: 10.1038/s41586-023-06685-2
57. Li R, Ferdinand JR, Loudon KW, Bowyer GS, Laidlaw S, Muyas F, et al. Mapping single-cell transcriptomes in the intra-tumoral and associated territories of kidney cancer. Cancer Cell. (2022) 40:1583–99 e10. doi: 10.1016/j.ccell.2022.11.001
58. Caronni N, La Terza F, Frosio L, and Ostuni R. IL-1β+ macrophages and the control of pathogenic inflammation in cancer. Trends Immunol. (2025) 46:403–15. doi: 10.1016/j.it.2025.03.001
59. Dawson CW, Tramountanis G, Eliopoulos AG, and Young LS. Epstein-barr virus latent membrane protein 1 (LMP1) activates the phosphatidylinositol 3-kinase/akt pathway to promote cell survival and induce actin filament remodeling. J Biol Chem. (2003) 278:3694–704. doi: 10.1074/jbc.M209840200
60. Su ZY, Siak PY, Leong C-O, and Cheah S-C. The role of Epstein–Barr virus in nasopharyngeal carcinoma. Front Microbiol. (2023) 14:1116143. doi: 10.3389/fmicb.2023.1116143
61. Zhu Y, Li J, Yang H, Yang X, Zhang Y, Yu X, et al. The potential role of m6A reader YTHDF1 as diagnostic biomarker and the signaling pathways in tumorigenesis and metastasis in pan-cancer. Cell Death Discov. (2023) 9:34. doi: 10.1038/s41420-023-01321-4
62. Gao Y, Pei G, Li D, Li R, Shao Y, Zhang QC, et al. Multivalent m(6)A motifs promote phase separation of YTHDF proteins. Cell Res. (2019) 29:767–9. doi: 10.1038/s41422-019-0210-3
63. Zaccara S and Jaffrey SR. A Unified Model for the Function of YTHDF Proteins in Regulating m(6)A-Modified mRNA. Cell. (2020) 181:1582–95 e18. doi: 10.1016/j.cell.2020.05.012
64. Liu T, Wei Q, Jin J, Luo Q, Liu Y, Yang Y, et al. The m6A reader YTHDF1 promotes ovarian cancer progression via augmenting EIF3C translation. Nucleic Acids Res. (2020) 48:3816–31. doi: 10.1093/nar/gkaa048
65. Wang J, Yin H, Li G, Wu D, Xu Y, Chen Y, et al. METTL14 promotes neuroblastoma formation by inhibiting YWHAH via an m6A-YTHDF1-dependent mechanism. Cell Death Discov. (2024) 10:186. doi: 10.1038/s41420-024-01959-8
Keywords: Nasopharyngeal carcinoma, METTL14, single-cell RNA sequencing, machine learning, TNF signaling, YWHAH, immune regulation
Citation: Liang Z, Zhou Z, Wang J, Shen L, Li Y, Zhu L, Hong G, Li Q, Xiao D, Lin X and Lin T (2026) METTL14 modulates the nasopharyngeal carcinoma microenvironment via m6A-modified YWHAH identified through single-cell and machine learning analyses. Front. Immunol. 17:1717039. doi: 10.3389/fimmu.2026.1717039
Received: 01 October 2025; Accepted: 19 January 2026; Revised: 12 January 2026;
Published: 02 February 2026.
Edited by:
Chao Yang, Zhejiang Ocean University, ChinaReviewed by:
Yingkun Xu, Shandong University, ChinaZhuangzhuang Geng, Children’s Hospital of Philadelphia, United States
Copyright © 2026 Liang, Zhou, Wang, Shen, Li, Zhu, Hong, Li, Xiao, Lin and Lin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Dong Xiao, eGlhb2RvbmdAc211LmVkdS5jbg==; Xiaolin Lin, bHhsMDY0NEBzbXUuZWR1LmNu; Taoyan Lin, bGludF95YW5AMTI2LmNvbQ==
†These authors have contributed equally to this work