Your privacy, your choice

We use essential cookies to make sure the site can function. We also use optional cookies for advertising, personalisation of content, usage analysis, and social media.

By accepting optional cookies, you consent to the processing of your personal data - including transfers to third parties. Some third parties are outside of the European Economic Area, with varying standards of data protection.

See our privacy policy for more information on the use of your personal data.

for further information and to change your choices.

Skip to main content

Comprehensive analysis of the relationship between RNA modification writers and immune microenvironment in head and neck squamous cell carcinoma

Abstract

Objectives

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer worldwide. Four types of RNA modification writers (m6A, m1A, A-I editing, and APA) are widely involved in tumorigenesis and the TME. We aimed to comprehensively explore the role of the four RNA modification writers in the progression and immune microenvironment of HNSCC.

Materials and methods

We first obtained transcription profile data and transcriptional variation of the four types of RNA modification writers from The Cancer Genome Atlas (TCGA) database. HNSCC patients in TCGA dataset were divided into different clusters based on the four types of RNA modification writers. Univariate Cox and Least absolute shrinkage and selection operator (LASSO) analyses were performed to conduct a Writer-score scoring system, which was successfully verified in the GSE65858 dataset and our clinical sample dataset. Finally, we evaluated the relationship between different RNA modification clusters (Writer-score) and immunological characteristics of HNSCC.

Results

Two different RNA modification clusters (A and B) were obtained. These RNA modification clusters (Writer-score) were strongly associated with immunological characteristics (immunomodulators, cancer immunity cycles, infiltrating immune cells (TIICs), inhibitory immune checkpoints, and T cell inflamed score (TIS)) of HNSCC.

Conclusions

This study identified two different RNA modification clusters and explored the potential relationship between RNA modification clusters (Writer-score) and immunological characteristics, offering a new theoretical basis for precision immunotherapy in patients with HNSCC.

Peer Review reports

Background

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common malignancy worldwide, with distant metastases and/or locoregional recurrence in 50% of cases [1,2,3]. Existing evidence indicates that epigenetic alterations are core components of tumorigenesis and progression, including HNSCC [4, 5]. RNA modification, an essential part of epigenetic alterations, plays an important role in the malignant progression of HNSCC. RNA modification is widely regarded as a relatively static trimmer of RNA structure and function [6]. In this study, we highlight several critical adenine-related RNA modifications, namely m1A, m6A, A-to-I and APA, which are mainly produced by the activity of ‘writer’ enzymes.

m1A refers to the methylation that occurs on the first nitrogen atom of adenosine in RNA, which had been found to exist in mRNAs, rRNAs, tRNAs, long non-coding RNAs (lncRNAs) and mitochondrial transcripts [7, 8]. The m1A ‘writer’, ‘reader’, and ‘eraser’ play various roles the posttranscriptional regulation of mRNAs and ncRNAs [9,10,11,12]. M6A, which occurs on the sixth nitrogen atom of adenosine in RNA, is the most prevalent ubiquitous RNA modification. The enzymes that recognize this modification site, ‘reader’, and the ones regulate the level of m1A, the ‘writer’ and ‘eraser’, regulate almost all processes of mRNA metabolism, including RNA transcription, translation, and degradation [13,14,15]. A-to-I is mediated by the adenosine deamination enzyme family, which occurs in double-stranded RNA (dsRNA) [16]. Recent studies have identified that A-to-I editing plays a vital role in mammalian survival, and its dysregulation may contribute to the development and progression of cancer [17]. APA, namely, alternative polyadenylation, is a phenomenon that generates distinct 3' termini on mRNAs and other transcripts generated by RNA polymerase II [18]. The typical core APA regulation complex consists of CPSFs, CSTFs, CFIm, and mammalian cleavage factor complex II (CFIIm) [19]. To fully understand the meaning of posttranscriptional epigenetic modifications, cross-talk between these four types of RNA modifications is urgently needed.

In recent years, with the rise of immunotherapy, immune checkpoint inhibitor (ICI) therapy has become a promising therapeutic method for HNSCC. ICIs have been approved as effective therapeutic strategies for recurrent and metastatic (R/M) HNSCC [20]. However, immune escape, involving perturbation of immune checkpoints, dysregulation of cytokines, and recruitment of inhibitory cell populations, limits the efficacy of ICI therapy [21]. Therefore, improving the efficacy of ICI therapy for HNSCC is challenging. Previous studies have shown that RNA modification writers are novel regulators of the tumor microenvironment [22] [23]. However, due to methodological limitations, these studies were limited to one or two RNA modification writers, and the anti-tumor immune effects of RNA modification writers are highly integrated. Therefore, a comprehensive analysis of the relationship between these four types of RNA modification writers and immune regulation is of great significance for accurate immunotherapy of HNSCC.

Materials and methods

Patients and datasets

The Workflow of this study is shown in Fig. 1.RNA sequencing data, somatic mutation data, Copy Number Variation (CNV) data, and clinical data of The Cancer Genome Atlas (TCGA)-HNSC were downloaded from the UCSC Xena Public Hub (http://xena.ucsc.edu/). As a training dataset in this study, the TCGA-HNSC dataset included 44 normal samples and 502 HNSCC samples after filtering patients without prognostic information. The mRNA-seq data and clinical information of the GSE65858 dataset, which included 270 HNSCC samples, were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/). Sixty-eight HNSCC samples were obtained from the tissue specimen Bank of Shengjing Hospital between 2015 and 2021. This study was approved by the Ethics Committee of Shengjing Hospital of China Medical University, and informed consent was obtained from all patients. All methods were executed in accordance with relevant guidelines and regulations. The GSE65858 and clinical sample datasets were used as the testing datasets. The baseline information for all these datasets is presented in Table S1A-C.

Fig. 1
figure 1

The Workflow of this study. TCGA. The Cancer Genome Atlas; HNSCC, Head and neck squamous cell carcinoma; DE, differentially expressed

Cell culture and transfection

HNSCC cells (SCC25) were obtained from the American Type Culture Collection (Manassas, VA, USA) and cultured in Dulbecco's modified Eagle (DMEM) medium/Ham’s F12 (1:1). Genechem Co., Ltd. (Shanghai, China) assisted in the synthesis of sequences for the knockdown of CYREN (CYREN-interfering sequence,5-TCTGGGAATCCTGATTGAGA-3′) and RDH1 0 (RDH10-interfering sequence,5-TACGATGCTGGAGATTAAT-3). Sixteen hours after transfection, the knockdown lentivirus (10 infections) was screened and a stable transfection cell line was confirmed. Cells were cultured at 37 °C in a 5% CO2 incubator.

Unsupervised clustering of RNA modification writers

RNA modification writers included four m1A enzymes, seven m6A enzymes, three A-to-I enzymes, and 12 APA enzymes. Unsupervised clustering analysis was performed to detect novel RNA modification writer clusters in HNSCC using the “ConsensuClusterPlus” R package based on the four types of RNA modification writers in the TCGA-HNSC dataset [24]. GSVA enrichment analysis was conducted using “GSVA” R package to explore the differential biological function between distinct RNA modification clusters [25].

Construction and validation of a Writer-score scoring system

We identified the RNA modification writer related DEGs using the “limma” R package based on the screening criteria adjusted p-value < 0.05 and |logFC|> 0.5. Univariate Cox regression analysis was conducted to identify DEGs with prognostic value. Next, the Least absolute shrinkage and selection operator (LASSO) regression analysis were performed to establish the Writer-score scoring system to quantify all individuals with HNSCC using the following formula: Writer-score = ∑(Coefi * Expri); here, i means the genes, Coefi means the coefficient of each gene and Expri means the expression level of each gene. Patients were classified into high- and low-Writer-score groups using the “survminer” R package. The robustness of the Writer-score scoring system was verified using the GSE65858 dataset and clinical sample dataset.

Quantitative Real-Time RT-PCR

Total RNA was extracted from HNSCC samples using the TRIzol reagent (Invitrogen, USA). The purity and concentration of the RNA samples were then measured. Subsequently, sample RNA was first generated into cDNA by reverse transcription, and real-time quantitative polymerase chain reaction (RT-qPCR) was performed. The 2-ΔΔCt method was used to calculate the relative expression levels of the genes, with glyceraldehyde-3-phosphate dehydrogenase (GAPDH) used as an internal reference. Primer sequences for the genes are presented in Table S2.

Western blotting

Cell lysates was applied to extract total proteins of HNSCC cells based on the procedure provided by the manufacturer. A bicinchoninic acid (BCA) method assay kit was used for protein quantification. The cell protein lysates were separated on a sodium dodecyl sulfate/polyacrylamide gel, transferred onto a polyvinylidene difluoride membrane, incubated with primary antibodies overnight at 4 °C, and then probed with a horseradish peroxidase-conjugated secondary antibody. The primary antibodies used in this study are listed in Supplementary Table S8. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal reference.

Construction of a nomogram model

A nomogram model was established a to predict the clinical outcome of HNSCC patients at 1-, 3- and 5-year using the “rms” R package. The calibration curve, decision curve analysis (DCA) curve, and receiver operating characteristic curve (ROC) were used to evaluate the performance of the nomogram [26].

Evaluated the correlation of RNA modifications clusters (Writer-score) with clinicopathological features and immunological characteristics

Immunological characteristics in our research included 122 immunomodulators (receptors, MHCs, immune stimulators, and chemokines) (Table S3A), tumor-infiltrating immune cells (TIICs), 36 effector genes of several TIICs (natural killer cells (NK), macrophages, dendritic cells (DCs), Th1 cells, and CD8 + T cells) (Table S3B), cancer immunity cycles, T cell inflamed score (TIS) (Table S3C), and 22 inhibitory immune checkpoints (Table S3D) [27]. The activities of these steps were calculated using single-sample gene set enrichment analysis (ssGSEA) [28]. "ESTIMATE" package in R software was used to estimate the ratio of immune matrix components in Immune microenvironment (TME), which was presented in the form of three scores: ImmuneScore, StromalScore and ESTIMATEScore were positively correlated with immunity, matrix and the sum of the two, respectively, which means that the higher the corresponding score, the greater the proportion of corresponding components in TME.

Results

Landscape of 26 RNA modification writers in TCGA-HNSC dataset

In the TCGA- HNSCC dataset, the somatic mutation frequency of 26 RNA modification writers was not very frequent in HNSCC. Of the 506 HNSCC samples, 60 (11.86%) contained mutations in RNA modification writers (Fig. 2A). Figure 2B shows the chromosomal sites of the 26 RNA-modification writers. We then assessed the CNV alterations of the 26 RNA modification writers and found that CPSF1, CLP1, and PCF11 had a higher frequency of CNV gain, whereas ZC3H13 and RBM15 had a higher proportion of CNV loss (Fig. 2C). Figure 2D indicated that most of the 26 RNA modification writers were more highly expressed in HNSCC samples than in the normal samples. Figure S1 have shown that only the expression of CSTF3 was positively correlated with the CNV gain of CSTF3 among the 26 RNA modification writers, indicating that the vast majority of CNV changes among the 26 RNA modification writers were unrelated to their expression. Pearson’s correlation analysis revealed that the majority of the 26 RNA modification writers were positively correlated with each other (Fig. 2E). The K-M curve revealed that CSTF1, CSTF3, PABPN1, and TRMT61A were poor prognostic factors for HNSCC, whereas ADARB2, ZC3H13, VIRMA, PCF11, and RBM15B were protective prognostic factors (Figure S2).

Fig. 2
figure 2

Landscape of 26 RNA modification writers in TCGA- HNSC dataset. (A) The mutation profiles of 26 RNA modification writers in 506 samples. (B) The chromosomal sites of the 26 RNA modification writers. (C) The copy number variation frequency of 26 RNA modification writers. (D) Differential expression histogram of the 26 RNA modification writers between HNSCC and normal samples. (E) The interactions between 26 RNA modification writers and their prognostic value for HNSCC. The circle size indicates the p value of the log-rank test, and the lines linking the 26 RNA modification writers indicate their interactions. *P < 0.05; **P < 0.01; ***P < 0.001

Unsupervised clustering based on 26 RNA modification writers

Consensus clustering analysis was performed on all HNSCC patient samples in the TCGA HNSCC dataset. K = 2 was found to be the best with the lowest proportion of ambiguous clustering (PAC) measures (Fig. 3A-C). All patients with HNSCC were classified into two clusters (clusters A and B). The patients in cluster A had poorer clinical outcomes than those in clusterB (Fig. 3B). Figure 3C indicated that the majority of the 26 RNA modification writers were more highly expressed in cluster B than in cluster A. PCA analysis showed that the 26 RNA modification writers could effectively separate the cluster A and cluster B group (Fig. 3D). GSVA enrichment analysis predicted that cluster B was markedly enriched in pathways in cancer, T cell receptor signaling pathway, and chemokine signaling pathway. Cluster A was mainly enriched in endocytosis and cytokine receptor interaction (Fig. 3E; Table S4).

Fig. 3
figure 3

Unsupervised clustering based on 26 RNA modification writers in TCGA- HNSC dataset. (A) Consensus matrices of the TCGA- HNSC dataset for k = 2. (B) CDF of consensus clustering at k = 2–10. Numbers next to the colors represent cluster numbers. (C) Relative changes in the area under the CDF curve at k = 2–10. (D) Survival analysis of the overall survival in clusterA and clusterB groups. (E) Expression heat map of the 26 RNA modification writers between clusterA and clusterB. (F) PCA plot of the clusterA and clusterB based on 26 RNA modification writers. (G) GSVA enrichment analysis in clusterA and clusterB groups. The heatmap was used to visualize these biological processes, and red represented activated pathways and green represented inhibited pathways. PCA, Principal Component Analysis. *P < 0.05; **P < 0.01; ***P < 0.001

Construction and validation of a Writer-score scoring system

Among the two clusters, we subsequently identified 9470 RNA modification-related DEGs (Table S5), and enrichment analysis was performed based on the DEGs. Next, 843 DEGs with significant prognostic value were selected to develop a Writer-score scoring system to quantify RNA modification patterns (Table S6). First, we placed the 843 DEGs into LASSO regression analysis, with 24 genes left as lambda = 0.082 and obtained the coefficients of genes in the Writer-score scoring system to calculate the Writer-score for each patient (Fig. 4A; Table S7). The patients in the TCGA-HNSCC dataset were divided into high- or low- groups according to an optimal cutoff value of the Writer-score, and patients in the high group had a poorer prognosis than those in the low- group (Fig. 4B). ROC curves were plotted to show the predictive performance of the Writer-score scoring system, the patients in high-group have poor outcome (p < 0.001). The AUC values at 1, 3- and 5-year were 0.62, 0.77, and 0.79, respectively (Fig. 4C). We then successfully verified our analysis results in the GSE65858 dataset and clinical sample dataset (Fig. 4D-G). Finally, univariate and multivariate Cox regression analyses of the training and testing datasets revealed that the Writer-score scoring system has an independent prognostic value (HR > 1, Writer-score represents risk factors for the prognosis of HNSCC) (Fig. 4H-J).

Fig. 4
figure 4

Construction and validation of a Writer-score scoring system. A LASSO regression was applied to establish the Writer-score scoring system. B Survival analysis of the high- and low- Writer-score groups in TCGA- HNSC dataset. C ROC curve to evaluate the predictive performance of Writer-score scoring system in TCGA- HNSC dataset. D Survival analysis of the high- and low- Writer-score groups in GSE65858 dataset. E ROC curve to evaluate the predictive performance of Writer-score scoring system in GSE65858 dataset. F Survival analysis of the high- and low- Writer-score groups in clinical specimen dataset. G ROC curve to evaluate the predictive performance of Writer-score scoring system in clinical specimen dataset. H Univariate and multivariate Cox regression analyses to evaluate the independent prognostic value of the Writer-score scoring system in TCGA- HNSC dataset. I Univariate and multivariate Cox regression analyses to evaluate the independent prognostic value of the Writer-score scoring system in GSE65858 dataset. J Univariate and multivariate Cox regression analyses to evaluate the independent prognostic value of the Writer-score scoring system in clinical specimen dataset. ROC, receiver operating characteristic; LASSO, least absolute shrinkage and selection operator

Construction of a nomogram model

A nomogram model combining the Writer-score and clinical information, including age, sex, grade, and stage, was constructed (Fig. 5A). The calibration curve at 1, 3- and 5-year clarified that the nomogram-predicted OS was close to the observed OS (Fig. 5B). DCA and ROC curves indicated that the nomogram model had the highest predictive performance compared with the individual features (Fig. 5C, D).

Fig. 5
figure 5

Construction of a nomogram model. A Construction of a nomogram model in TCGA- HNSC dataset. B Calibration curves of the nomogram model of 1-, 3-, 5-year. (C) DCA curves for predicting the net benefit of the nomogram model. D ROC curves for predicting the predictive performance of the nomogram model. DCA, clinical decision curve

Evaluated the correlation of RNA modifications clusters (Writer-score) with clinicopathological features and immunological characteristics

Samples in the low Writer-score group had a higher mutation frequency (91%, Fig. 6A), those in the high Writer-score group (89%, Fig. 6B). We found that the Writer-score had no statistical difference in clinicopathological features (Fig. 6C-H, Table 1). Correlation analysis indicated that Writer-score was positively correlated with mRNAsi (Fig. 6I). The patients in cluster A had a higher Writer-score (Fig. 6J, Table 1), and the Sankey diagram shows the relationship between cluster, Writer-score, fustat (Fig. 6K).

Fig. 6
figure 6

Evaluated the correlation of Writer-score with clinicopathological features. A, B The mutation profiles in high- and low- Writer-score group. C-J Writer-score of the HNSCC patients are classified by age, sex, stage, grade, alcohol_history, Hpv16_status, mRNAsi, cluster. K Sankey diagram shows the relationships among cluster, Writer-score, survival status

Table 1 Relationship between Writer-score and tumor characteristics in ovarian cancer patients

Subsequently, we investigated the relationship between the RNA modification clusters (Writer-score) and immunological characteristics. We found that the majority of immunomodulators were overexpressed in cluster B (low Writer-score) (FiguresS3A, S3B). The comprehensive performance of immunomodulators can directly determine the activities of the cancer immunity cycle. Correspondingly, we found that the majority of the steps in the cancer immunity cycle were activated in cluster B (low Writer-score) (Figs. 7A and 8A). Consequently, the activities of these cancer immunity cycles lead to the recruitment of TIICs (Figs. 7B, 8B, 9B) and expression levels of their effector genes (Fig. 7D and 8D). Moreover, we also observed that significant genes in the Writer-score scoring system had a strong correlation with the TIICs in the HNSCC TME (Fig. 9A). Among them, TRBV7-3 was positively associated with the majority of TIICs, whereas NEAT1was negatively associated with the majority of TIICs. In addition, the patients in the low-risk group, cluster B (low Writer-score), also had higher stromal and immune scores and ESTIMATEScore (Fig. 7C, 8C). We analyzed the relationships between RNA modification clusters (Writer-score) and 22 inhibitory immune checkpoints and found that most of the inhibitory immune checkpoints were overexpressed in cluster B (low Writer-score) (Fig. 7E-G and 8E-G). We also observed that patients in cluster B (low Writer-score) had a higher TIS (Fig. 7H, 8H). Therefore, these findings suggest that cluster B (low Writer-score) may be an inflamed phenotype and may be more sensitive to ICB treatment. In addition, we found that NEAT1 was negatively correlated with the expression of M2-polarization biomarkers (CD163 and CD206), suggesting that NEAT1 is involved in the M2-polarization of HNSCC (Fig. 10A and B). Finally, we investigated the correlation between the Writer-score and CNV gain or loss of the 26 RNA modification writers and found that Writer-score was only positively correlated with the CNV gain of CSTF1(Fig. 10C). Western blotting revealed that the downregulation of CYREN significantly inhibited the expression of DNA damage response (DDR) related proteins (NBS1, Mre11, RAD51, XRCC4, Ku70, and Ku80) (Fig. 11A), whereas the downregulation of RDH10 inhibited the expression of TGF-β/SMAD signaling pathway-related proteins (TGF-β, phosphorylated SMAD2, and phosphorylated SMAD3) (Fig. 11B).

Fig. 7
figure 7

Evaluated the correlation of RNA modifications clusters with immunological characteristics. A Heat map of the cancer immunity cycles between clusterA and clusterB. B Differential histogram of the TIICs between clusterA and clusterB. C Differential histogram of the stromal score, immune score and ESTIMATEScore between clusterA and clusterB. D Expression heat map of the effector genes between clusterA and clusterB. E Differential histogram of 22 immune checkpoints between clusterA and clusterB. F Differential histogram of CD274 between clusterA and clusterB. G Differential histogram of CTLA4 between clusterA and clusterB. H Differential histogram of TIS between clusterA and clusterB

Fig. 8
figure 8

Evaluated the correlation of Writer-score with immunological characteristics. A Heat map of the cancer immunity cycles between high- and low- Writer-score group. B Differential histogram of the TIICs between high- and low- Writer-score group. C Differential histogram of the stromal score, immune score and ESTIMATEScore between high- and low- Writer-score group. D Expression heat map of the effector genes between high- and low- Writer-score group. E Differential histogram of 22 immune checkpoints between high- and low- Writer-score group. F Differential histogram of CD274 high- and low- Writer-score group. G Differential histogram of CTLA4 between high- and low- Writer-score group. H Differential histogram of TIS between high- and low- Writer-score group

Fig. 9
figure 9

Evaluation of the TME between high- and low- Writer-score group. A Correlations between the 23 types of immune cells and the 12 genes in the Writer-score scoring system. B Correlations between Writer-score and 23 types of immune cells. TME, tumor microenvironment

Fig. 10
figure 10

Pearson correlation analysis. A NEAT1 was negatively correlated with the expression of M2-polarization biomarker CD163. B NEAT1 was negatively correlated with the expression of M2-polarization biomarker CD206. C Writer-score was positively correlated with the CNV gain of CSTF1

Fig. 11
figure 11

Western blotting analyzed the functions of CYREN and RDH10 in HNSCC. A Effect of CYREN downregulation on the expression of NBS1, Mre11, RAD51, XRCC4, Ku70, and Ku80. B Effect of RDH10 downregulation on the expression of TGF-β, SMAD2, p-SMAD2, SMAD3, and p-SMAD3

Discussion

Emerging evidence indicates that RNA modification writers play pivotal roles in innate immunity, inflammation, and tumor progression. However, previous publications have focused on single RNA modification writers, and the underlying interrelationships and functions of multiple RNA modification writers are not fully understood. Thus, our study revealed four types of RNA modification writers at the transcriptional level in HNSCC cells.

The writer-scoring system was constructed based on 12 genes (CYREN, NAV1, STING1, WDR41, RDH10, PRDX6, NEAT1, SLC20A1, VPS33B, CELSR3, LINC02487, and TRBV7-3). A review of previous research has shown that most of these genes are associated with tumorigenesis and immunity. The modulator of retrovirus infection (CYREN or MRI) is a 30-kDa protein with a C-terminal XLF-like motif (XLM) and a conserved N-terminal Ku-binding motif (KBM), which plays an important role in the DNA damage response (DDR) [29, 30]. NAV1 (neuron navigator 1), a member of the neuron navigator family, has been reported to be associated with directional neuron migration, diabetes mellitus, calcific aortic valve stenosis, and tumorigenesis [31,32,33,34]. Stimulator of interferon genes (STING) is an important adaptor in innate immunity that can recognize pathogens and self-DNA through a cyclic GMP-AMP synthase (cGAS)-stimulator [35]. Numerous studies have shown that the STING pathway plays an essential role in anti-tumor immunity and tumor immune surveillance by inducing the secretion of inflammatory cytokines [36]. WD-repeat protein 41 (WDR41) is a protein with six WD40‐repeat domains, that is involved in the regulation of dopamine signalling and autophagy-lysosomal pathway [37, 38]. Wang et al. [39]reported that WDR41 is expressed at low levels in breast cancer and promotes tumor characteristics of breast cancer cells, including cell viability, cell cycle, and migration. Retinol dehydrogenase 10 (RDH10), a member of the short-chain dehydrogenase/reductase family, is essential for the visual cycle of retinoid A. Previous studies have shown that RDH10 plays a vital role in the occurrence and development of tumors. Peroxiredoxin 6 (PRDX6), a member of the peroxidase family, is a bifunctional enzyme with peroxidase and phospholipase A2 (PLA2) activities [40]. Zhang et al. [41]found that PRDX6 was downregulated in colon adenocarcinoma (COAD) and was involved in immune infiltration in COAD. Nuclear enriched abundant transcript 1 (NEAT1) is a nuclear paraspeckle located and transcribed from chromatin 11 [42]. It has been reported to correlate with various biological processes, including infection, tumorigenesis, immunity, and neuropathy [43]. In our study, we found that NEAT1 was negatively correlated with the expression of M2-polarization biomarkers (CD163 and CD206) and the infiltration of most immune cell types; therefore, we believe that NEAT1 is involved in M2-polarization and immune escape of HNSCC. Solute carrier family 20 member 1 (SLC20A1) is a sodium/inorganic-dependent phosphate symporter that plays a critical role in cell proliferation and tumor growth [44]. Vacuolar protein sorting-associated 33 B (VPS33B), a member of the Sec1 domain family, plays an important role in platelet activation, in vivo thrombosis, homeostasis, and megakaryocyte biogenesis [45]. Long non-coding RNA LINC02487, located on chromosome 6 and 2557 bp long, can be used as a potential biomarker for the diagnosis and prognosis of oral cancer [46]. However, there have been no reports of CELSR3 and TRBV7-3.

Blocking immune checkpoints can help rescue T-cell responses and promote the development of immunotherapy. Recent studies have revealed that m6A regulation could explain some of the underlying mechanisms of immune regulation. Zhang X et al. discovered that m6A regulator-mediated RNA methylation modification patterns plays a crucial role in immune microenvironment regulation of periodontitis [47]. The deletion of METTL3 in T cells was found to participate in the homeostatic differentiation of T cells [48]. In our research, we also found that HNSCC patients in the cluster B (low Writer-score) group had higher immunological characteristics (immunomodulators, cancer immunity cycles, TIICs, inhibitory immune checkpoints, and TIS). In addition, we found that genes in the Writer-score scoring system were closely correlated with TIICs. These results suggest that cluster B (low Writer-score) may be an inflamed phenotype and be more sensitive to ICB treatment, and the genes in the Writer-score scoring system may be potential targets for improving the response rate to immunotherapy in HNSCC patients.

Our research has some shortcomings: (1) This is a retrospective study and an independent prospective cohort is required to validate the risk model constructed. (2) This study relies heavily on data sets and computational predictions, with little experimental validation. In future, more experimental validation should be done to enhance the robustness and impact of the study.

Conclusions

In conclusion, we established and validated a reliable Writer-score scoring system based on four types of RNA modification writers, which exhibited superior performance across multiple cohorts to predict the patient prognosis and demonstrated close associations with immunotherapy response. This study provides a foundation for prognostic prediction of HNSCC patients through integrating multiomics data and frontier computational algorithms.

Data availability

All data generated or analyzed during this study are included in this published article and its supplementary information files.

References

  1. Alam MM, Fermin JM, Knackstedt M, Noonan MJ, Powell T, Goodreau L, Daniel EK, Rong X, Moore-Medlin T, Khandelwal AR, et al. Everolimus downregulates STAT3/HIF-1α/VEGF pathway to inhibit angiogenesis and lymphangiogenesis in TP53 mutant head and neck squamous cell carcinoma (HNSCC). Oncotarget. 2023;14:85–95.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Machiels JP, René Leemans C, Golusinski W, Grau C, Licitra L, Gregoire V. Squamous cell carcinoma of the oral cavity, larynx, oropharynx and hypopharynx: EHNS-ESMO-ESTRO Clinical Practice Guidelines for diagnosis, treatment and follow-up. Ann Oncol. 2020;31(11):1462–75.

    Article  PubMed  Google Scholar 

  3. Zhang X, Feng H, Li Z, Li D, Liu S, Huang H, Li M. Application of weighted gene co-expression network analysis to identify key modules and hub genes in oral squamous cell carcinoma tumorigenesis. Onco Targets Ther. 2018;11:6001–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Widschwendter M, Jones A, Evans I, Reisel D, Dillner J, Sundström K, Steyerberg EW, Vergouwe Y, Wegwarth O, Rebitschek FG, et al. Epigenome-based cancer risk prediction: rationale, opportunities and challenges. Nat Rev Clin Oncol. 2018;15(5):292–309.

    Article  PubMed  Google Scholar 

  5. Tian X, Shi C, Liu S, Zhao C, Wang X, Cao Y. Methylation related genes are associated with prognosis of patients with head and neck squamous cell carcinoma via altering tumor immune microenvironment. J Dent Sci. 2023;18(1):57–64.

    Article  PubMed  Google Scholar 

  6. Jonkhout N, Tran J, Smith MA, Schonrock N, Mattick JS, Novoa EM. The RNA modification landscape in human disease. RNA. 2017;23(12):1754–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Zhang C, Jia G. Reversible RNA Modification N(1)-methyladenosine (m(1)A) in mRNA and tRNA. Genomics Proteomics Bioinformatics. 2018;16(3):155–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Li J, Zhang H, Wang H. N(1)-methyladenosine modification in cancer biology: Current status and future perspectives. Comput Struct Biotechnol J. 2022;20:6578–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Safra M, Sas-Chen A, Nir R, Winkler R, Nachshon A, Bar-Yaacov D, Erlacher M, Rossmanith W, Stern-Ginossar N, Schwartz S. The m1A landscape on cytosolic and mitochondrial mRNA at single-base resolution. Nature. 2017;551(7679):251–5.

    Article  CAS  PubMed  Google Scholar 

  10. Zhou H, Rauch S, Dai Q, Cui X, Zhang Z, Nachtergaele S, Sepich C, He C, Dickinson BC. Evolution of a reverse transcriptase to map N(1)-methyladenosine in human messenger RNA. Nat Methods. 2019;16(12):1281–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Chen Z, Qi M, Shen B, Luo G, Wu Y, Li J, Lu Z, Zheng Z, Dai Q, Wang H. Transfer RNA demethylase ALKBH3 promotes cancer progression via induction of tRNA-derived small RNAs. Nucleic Acids Res. 2019;47(5):2533–45.

    Article  CAS  PubMed  Google Scholar 

  12. Chujo T, Suzuki T. Trmt61B is a methyltransferase responsible for 1-methyladenosine at position 58 of human mitochondrial tRNAs. RNA. 2012;18(12):2269–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA Modifications in Gene Expression Regulation. Cell. 2017;169(7):1187–200.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Nachtergaele S, He C. Chemical Modifications in the Life of an mRNA Transcript. Annu Rev Genet. 2018;52:349–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Zhao W, Qi X, Liu L, Ma S, Liu J, Wu J. Epigenetic Regulation of m(6)A Modifications in Human Cancer. Mol Ther Nucleic Acids. 2020;19:405–12.

    Article  CAS  PubMed  Google Scholar 

  16. Goncharov AO, Shender VO, Kuznetsova KG, Kliuchnikova AA, Moshkovskii SA. Interplay between A-to-I Editing and Splicing of RNA: A Potential Point of Application for Cancer Therapy. Int J Mol Sci. 2022;23(9):5240.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Ramírez-Moya J, Miliotis C, Baker AR, Gregory RI, Slack FJ, Santisteban P. An ADAR1-dependent RNA editing event in the cyclin-dependent kinase CDK13 promotes thyroid cancer hallmarks. Mol Cancer. 2021;20(1):115.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Tian B, Manley JL. Alternative polyadenylation of mRNA precursors. Nat Rev Mol Cell Biol. 2017;18(1):18–30.

    Article  CAS  PubMed  Google Scholar 

  19. Shi Y, Di Giammartino DC, Taylor D, Sarkeshik A, Rice WJ, Yates JR 3rd, Frank J, Manley JL. Molecular architecture of the human pre-mRNA 3’ processing complex. Mol Cell. 2009;33(3):365–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Gao L, Zhang A, Yang F, Du W. Immunotherapeutic Strategies for Head and Neck Squamous Cell Carcinoma (HNSCC): Current Perspectives and Future Prospects. Vaccines (Basel). 2022;10(8):1272.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Yu C, Li Q, Zhang Y, Wen ZF, Dong H, Mou Y. Current status and perspective of tumor immunotherapy for head and neck squamous cell carcinoma. Front Cell Dev Biol. 2022;10:941750.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Wei Z, Chen Y, Zeng Z, Peng Y, Li L, Hu N, Gao X, Cai W, Yin L, Xu Y, et al. The novel m6A writer METTL5 as prognostic biomarker probably associating with the regulation of immune microenvironment in kidney cancer. Heliyon. 2022;8(12):e12078.

    Article  Google Scholar 

  23. Han D, Liu J, Chen C, Dong L, Liu Y, Chang R, Huang X, Liu Y, Wang J, Dougherty U, et al. Anti-tumour immunity controlled through mRNA m(6)A methylation and YTHDF1 in dendritic cells. Nature. 2019;566(7743):270–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Yu J, Chen Y, Pan X, Wen W. Relationships of Ferroptosis and Pyroptosis-Related Genes with Clinical Prognosis and Tumor Immune Microenvironment in Head and Neck Squamous Cell Carcinoma. Oxid Med Cell Longev. 2022;2022:3713929.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Hu J, Yu A, Othmane B, Qiu D, Li H, Li C, Liu P, Ren W, Chen M, Gong G, et al. Siglec15 shapes a non-inflamed tumor microenvironment and predicts the molecular subtype in bladder cancer. Theranostics. 2021;11(7):3089–108.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, Yuan H, Cheng P, Li F, Long Z, et al. TIP: A Web Server for Resolving Tumor Immunophenotype Profiling. Cancer Res. 2018;78(23):6575–80.

    Article  CAS  PubMed  Google Scholar 

  29. Agarwal S, Harada J, Schreifels J, Lech P, Nikolai B, Yamaguchi T, Chanda SK, Somia NV. Isolation, characterization, and genetic complementation of a cellular mutant resistant to retroviral infection. Proc Natl Acad Sci U S A. 2006;103(43):15933–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Grundy GJ, Rulten SL, Arribas-Bosacoma R, Davidson K, Kozik Z, Oliver AW, Pearl LH, Caldecott KW. The Ku-binding motif is a conserved module for recruitment and stimulation of non-homologous end-joining proteins. Nat Commun. 2016;7:11242.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Martínez-López MJ, Alcántara S, Mascaró C, Pérez-Brangulí F, Ruiz-Lozano P, Maes T, Soriano E, Buesa C. Mouse neuron navigator 1, a novel microtubule-associated protein involved in neuronal migration. Mol Cell Neurosci. 2005;28(4):599–612.

    Article  PubMed  Google Scholar 

  32. Fei H, Shi M, Chen L, Wang Z, Suo L. MicroRNA-18 promotes apoptosis of islet β-cells via targeting NAV1. Exp Ther Med. 2019;18(1):389–96.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. Thériault S, Dina C, Messika-Zeitoun D, Le Scouarnec S, Capoulade R, Gaudreault N, Rigade S, Li Z, Simonet F, Lamontagne M, et al. Genetic Association Analyses Highlight IL6, ALPL, and NAV1 As 3 New Susceptibility Genes Underlying Calcific Aortic Valve Stenosis. Circ Genom Precis Med. 2019;12(10): e002617.

    Article  PubMed  Google Scholar 

  34. Li L, Lee KM, Han W, Choi JY, Lee JY, Kang GH, Park SK, Noh DY, Yoo KY, Kang D. Estrogen and progesterone receptor status affect genome-wide DNA methylation profile in breast cancer. Hum Mol Genet. 2010;19(21):4273–7.

    Article  CAS  PubMed  Google Scholar 

  35. Lee JJ, Kim SY, Kim SH, Choi S, Lee B, Shin JS. STING mediates nuclear PD-L1 targeting-induced senescence in cancer cells. Cell Death Dis. 2022;13(9):791.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Woo SR, Fuertes MB, Corrales L, Spranger S, Furdyna MJ, Leung MY, Duggan R, Wang Y, Barber GN, Fitzgerald KA, et al. STING-dependent cytosolic DNA sensing mediates innate immune recognition of immunogenic tumors. Immunity. 2014;41(5):830–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Yeo B, Turner NC, Jones A. An update on the medical management of breast cancer. BMJ. 2014;348: g3608.

    Article  PubMed  Google Scholar 

  38. Zwart W, Terra H, Linn SC, Schagen SB. Cognitive effects of endocrine therapy for breast cancer: keep calm and carry on? Nat Rev Clin Oncol. 2015;12(10):597–606.

    Article  PubMed  Google Scholar 

  39. Wang H, Wu D, Cai L, Li X, Zhang Z, Chen S. Aberrant methylation of WD-repeat protein 41 contributes to tumour progression in triple-negative breast cancer. J Cell Mol Med. 2020;24(12):6869–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Li H, Zhang D, Li B, Zhen H, Chen W, Men Q. PRDX6 Overexpression Promotes Proliferation, Invasion, and Migration of A549 Cells in vitro and in vivo. Cancer Manag Res. 2021;13:1245–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Zhang X, Gao F, Li N, Zhang J, Dai L, Yang H. Peroxiredoxins and Immune Infiltrations in Colon Adenocarcinoma: Their Negative Correlations and Clinical Significances, an In Silico Analysis. J Cancer. 2020;11(11):3124–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Xia KG, Wang CM, Shen DY, Song XY, Mu XY, Zhou JW, Zhu AY, Xuan Q, Tao T. LncRNA NEAT1-associated aerobic glycolysis blunts tumor immunosurveillance by T cells in prostate cancer. Neoplasma. 2022;69(3):594–602.

    Article  CAS  PubMed  Google Scholar 

  43. Liu RJ, Xu ZP, Li SY, Yu JJ, Feng NH, Xu B, Chen M. BAP1-Related ceRNA (NEAT1/miR-10a-5p/SERPINE1) Promotes Proliferation and Migration of Kidney Cancer Cells. Front Oncol. 2022;12: 852515.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Sato K, Akimoto K. Expression Levels of KMT2C and SLC20A1 Identified by Information-theoretical Analysis Are Powerful Prognostic Biomarkers in Estrogen Receptor-positive Breast Cancer. Clin Breast Cancer. 2017;17(3):e135–42.

    Article  CAS  PubMed  Google Scholar 

  45. Ning Y, Zeng Z, Deng Y, Feng W, Huang L, Liu H, Lin J, Zhang C, Fan Y, Liu L. VPS33B interacts with NESG1 to suppress cell growth and cisplatin chemoresistance in ovarian cancer. Cancer Sci. 2021;112(5):1785–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Li Y, Cao X, Li H. Identification and Validation of Novel Long Non-coding RNA Biomarkers for Early Diagnosis of Oral Squamous Cell Carcinoma. Front Bioeng Biotechnol. 2020;8:256.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Zhang X, Zhang S, Yan X, Shan Y, Liu L, Zhou J, Kuang Q, Li M, Long H, Lai W. m6A regulator-mediated RNA methylation modification patterns are involved in immune microenvironment regulation of periodontitis. J Cell Mol Med. 2021;25(7):3634–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Li HB, Tong J, Zhu S, Batista PJ, Duffy EE, Zhao J, Bailis W, Cao G, Kroehling L, Chen Y, et al. m(6)A mRNA methylation controls T cell homeostasis by targeting the IL-7/STAT5/SOCS pathways. Nature. 2017;548(7667):338–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank the authors who provided the TCGA and GEO public datasets. We thank LetPub (www.letpub.com) for linguistic assistance during the preparation of this manuscript.

Clinical Trial Number

Not applicable.

Funding

Not applicable.

Author information

Authors and Affiliations

Authors

Contributions

Wei Li, Ying Chen, Yao Zhang, Wen Wen, Yingying Lu conceived of and designed the study. Wei Li, Ying Chen, Yao Zhang, Wen Wen, Yingying Lu developed this methodology. Wei Li, Ying Chen, Yao Zhang, Wen Wen, Yingying Lu analyzed and interpreted the data. Wei Li, Ying Chen, Yao Zhang, Wen Wen, Yingying Lu wrote, reviewed, and revised the manuscript.

Corresponding author

Correspondence to Yingying Lu.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the Ethics Committee of Shengjing Hospital of China Medical University, and informed consent was obtained from all patients. All methods were performed in accordance with relevant guidelines and regulations.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Li, W., Chen, Y., Zhang, Y. et al. Comprehensive analysis of the relationship between RNA modification writers and immune microenvironment in head and neck squamous cell carcinoma. BMC Immunol 25, 76 (2024). https://doi.org/10.1186/s12865-024-00667-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12865-024-00667-3

Keywords