ORIGINAL RESEARCH article
Sec. Epigenomics and Epigenetics
The Prognostic Value and Immune Landscapes of a m6A/m5C/m1A-Related LncRNAs Signature in Head and Neck Squamous Cell Carcinoma
- 1Department of Otorhinolaryngology, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
- 2Department of Stomatology, Zhongshan Hospital, Fudan University, Shanghai, China
Background: N6-methyladenosine (m6A), 5-methylcytosine (m5C) and N1-methyladenosine (m1A) are the main RNA methylation modifications involved in the progression of cancer. However, it is still unclear whether m6A/m5C/m1A-related long non-coding RNAs (lncRNAs) affect the prognosis of head and neck squamous cell carcinoma (HNSCC).
Methods: We summarized 52 m6A/m5C/m1A-related genes, downloaded 44 normal samples and 501 HNSCC tumor samples with RNA-seq data and clinical information from The Cancer Genome Atlas (TCGA) database, and then searched for m6A/m5C/m1A-related genes co-expressed lncRNAs. We adopt the least absolute shrinkage and selection operator (LASSO) Cox regression to obtain m6A/m5C/m1A-related lncRNAs to construct a prognostic signature of HNSCC.
Results: This prognostic signature is based on six m6A/m5C/m1A-related lncRNAs (AL035587.1, AC009121.3, AF131215.5, FMR1-IT1, AC106820.5, PTOV1-AS2). It was found that the high-risk subgroup has worse overall survival (OS) than the low-risk subgroup. Moreover, the results showed that most immune checkpoint genes were significantly different between the two risk groups (p < 0.05). Immunity microenvironment analysis showed that the contents of NK cell resting, macrophages M2, and neutrophils in samples of low-risk group were significantly lower than those of high-risk group (p < 0.05), while the contents of B cells navie, plasma cells, and T cells regulatory (Tregs) were on the contrary (p < 0.05). In addition, patients with high tumor mutational burden (TMB) had the worse overall survival than those with low tumor mutational burden.
Conclusion: Our study elucidated how m6A/m5C/m1A-related lncRNAs are related to the prognosis, immune microenvironment, and TMB of HNSCC. In the future, these m6A/m5C/m1A-related lncRNAs may become a new choice for immunotherapy of HNSCC.
Head and neck squamous cell carcinoma (HNSCC) represents the sixth most common malignant tumor in the world (Johnson et al., 2020). These tumors originate from the mucosa of various anatomical parts of the oral cavity, pharynx, larynx and sinonasal tract and are heterogeneous in etiology and clinicopathological characteristics, and these highly aggressive malignant tumors affect more than 830,000 patients worldwide each year (Johnson et al., 2020; Qin et al., 2021). HNSCC is usually associated with exposure to tobacco-derived carcinogens, excessive alcohol consumption, or both. Tumors appearing in the oropharynx are increasingly related to human papillomavirus (HPV) carcinogenic strains, mainly HPV-16, and to a lesser extent HPV-18 and other strains of infection (Riva et al., 2021). In current treatment methods, only about 50% of HNSCC patients survive for more than 5 years. Therefore, it is urgent to analyze the heterogeneity of HNSCC to choose different treatment options. In particular, there is an urgent need for reliable biomarkers with predictive and prognostic value (Rasheduzzaman et al., 2020).
Some evidences show that epigenetic changes (including: DNA methylation, histone covalent modification, chromatin remodeling, non-coding RNA (ncRNA) activity and RNA chemical modifications) are often related to oral carcinogenesis, tumor malignancy and resistance to treatment (Castilho et al., 2017; Romanowska et al., 2020). Increasing evidence suggests that dynamic RNA modifications pathways are also misregulated in human cancers (including HNSCC, lung cancer, breast cancer, pancreatic cancer, endometrial cancer, et al.) and may be ideal targets of cancer therapy (Barbieri and Kouzarides, 2020; Zhou et al., 2020; Huang et al., 2021). Nowadays, over 100 types of chemical modifications have been identified in cellular RNAs, including messenger RNAs (mRNAs), transfer RNAs (tRNAs), ribosomal RNAs (rRNAs), and ncRNAs, most of these modifications contribute to pre-mRNA splicing, nuclear exporting, transcript stability and translation initiation of eukaryotic cellular RNAs. Among these modifications, methylation modification is the most abundant which include N6-methyladenosine (m6A), 5-methylcytosine (m5C), N1-methyladenosine (m1A) (Shi et al., 2020). M6A is the most abundant form of methylation modification in eukaryotic cellular RNAs, and it is also the most thoroughly studied type of RNA modification. The sequence near the m6A modification site on the mRNA is highly conservative, mainly occurring on the adenine of the RRACH (R: purine; A: m6A; H: non-guanine) sequence, and its function is determined by the methyltransferase (writers: methyltransferase like 3 (METTL3), methyltransferase like 14 (METTL14) and WT1 associated protein (WTAP), et al.), demethylase (erasers: fat mass and obesity associated protein (FTO) and alkB homologue 5 (ALKBH5), et al.) and binding protein (readers: YTH domain containing family, the insulin like growth factor 2 mRNA binding protein (IGF2BP) family and heterogeneous nuclear ribonucleoproteins A2/B1 (HNRNPA2B1), et al.) are jointly regulated (Yang et al., 2018; Li M. et al., 2021). M6A modification is involved in a variety of biological processes, such as stem cell differentiation, cell division, gametogenesis and biological rhythms, et al., as well as the occurrence of a variety of diseases, including tumors, obesity and infertility (Xu Y. et al., 2020; Zhang M. et al., 2020; Song et al., 2020). M5C modification is another prevalent RNA modification in multiple RNA species, including mRNAs, tRNAs, rRNAs, and ncRNAs (Chen et al., 2021). M5C modification is distributed with mRNA, enriched around 5′UTR and 3′UTR, and conserved in tRNA and rRNA. It is dynamically regulated by related enzymes, including methyltransferase (writers: NOL1/NOP2/SUN domain family member (NSUN) family and DNA methyltransferase like 2 (DNMT2), et al.), demethylases (erasers: tet methylcytosine dioxygenase 2 (TET2), et al.), and binding proteins (readers: Aly/REF export factor (ALYREF), et al.) (Wiener and Schwartz, 2021). Up to now, accumulated studies have shown that m5C is involved in various RNA metabolism, including mRNA export, RNA stability and translation. The absence of m5C modification in organisms may lead to mitochondrial dysfunction, stress response defects, gametogenesis and embryogenesis frustration, and abnormal development of nerves and brain, which is related to cell migration and tumorigenesis (Chellamuthu and Gray, 2020; Vissers et al., 2020). The majority of these m1A sites are located within the 5′-untranslated region (5′ UTR), m1A modification has been identified in tRNA, rRNA, mRNA, long non-coding RNA (lncRNA) and mitochondrial (mt) transcripts, It can also be dynamically regulated by methyltransferases, demethylases and binding proteins (Zhang and Jia, 2018).
LncRNAs are a class of RNA molecules with transcript length exceeding 200nt. it is generally believed that they do not encode proteins, but participate in the regulation of protein-encoding genes in many aspects (epigenetic regulation, transcriptional regulation and post-transcriptional regulation, etc.) in the form of RNA. More and more lncRNAs have been found to play a role in the initiation and development of HNSCC, suggesting that they can be used as novel biomarkers and therapeutic targets to provide more effective diagnosis, prognosis and treatment for HNSCC patients, such as HOTAIR, UCA1, FOXC1, AFAP1-AS1, et al. (Cossu et al., 2019; Wang S. et al., 2020; Wang et al., 2020c.). Researchers have found that methylation modification of lncRNA regulates numerous processes affecting cancer cell activity by enhancing and increasing lncRNA stability, increasing lncRNA nuclear accumulation, and promoting decay of some lncRNAs (Dai et al., 2020). However, there is currently no evidence to explain the correlation between three major RNA methylation modifications related lncRNAs and HNSCC. In this study, we integrated 44 normal samples and 501 tumor samples with clinical trials data from The Cancer Genome Atlas (TCGA). We used bioinformatics and statistical analysis methods to explore whether m6A/m5C/m1A-related lncRNAs were related to the diagnosis and prognosis of HNSCC.
2. Materials and Methods
2.1. Data Collection and Processing
RNA-seq data and clinical information of HNSCC patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/) (Wang et al., 2020d), which contained data of 501 HNSCC and 44 non-tumor tissues. Raw count data was downloaded for further analysis. In addition, based on previous publications, the expression matrix of 52 m6A/m5C/m1A-related genes was extracted from TCGA database, including expression data on writers, erasers and readers (Table 1). We obtain the gene annotation file of Genome Reference Consortium Human Build 38 (GRCh38) from the GENCODE website (https://www.gencodegenes.org/human/) for the annotation of the mRNAs in the TCGA database (Frankish et al., 2019).
We downloaded 33 types of tumor HTSeq-FPKM GDC Hub files, Phenotype GDC Hub files, survival data GDC Hub files from the UCSC Xena database (http://xena.ucsc.edu/) (Yang et al., 2019), and downloaded immune subtype files from the TCGA Pan-Cancer (PANCAN) database. We downloaded the nci 60 RNA seq files from Processed Data Set and Compound activity: DTP NCI-60 files from the CellMiner database (http://discover.nci.nih.gov/cellminer/home.do) (Reinhold et al., 2019).
2.2. Statistical Analysis and Bioinformatic Analysis
2.2.1. mRNAs–lncRNAs Co-Expression Network and Prognostic Value of m6A/m5C/m1A-Related lncRNAs
After obtaining the expression matrix of 52 m6A/m5C/m1A-related genes and expression matrix of lncRNAs from the TCGA database, we screened out the lncRNAs co-expressed with m6A/m5C/m1A-related genes by setting the absolute value of correlation coefficient > 0.5 and p value < 0.001. The univariate cox regression analysis was performed to screen out prognostic value of m6A/m5C/m1A-related lncRNAs.
2.2.2. Cluster the Samples Based on the m6A/m5C/m1A-Related lncRNAs
Based on the expression level of these 44 lncRNAs with prognostic value in tumor samples, we used the R package “ConsensusClusterPlus” to classify tumor samples into different groups and to analyze the prognosis by cluster analysis of samples in the database (Zheng et al., 2020). Principal Component Analysis (PCA) was applied to the typed samples to observe whether the samples could be intuitively distinguished (Rapa et al., 2021). In order to understand the prognosis of samples with different classifications, the R package “survival” and “Survminer” (Zhang L.-P. et al., 2020) were used to analyze the prognosis of the samples after classifications. Clinical data were included to observe whether different clinical features differed in different classifications. The NMT classification comes from the TNM classification of malignant tumours (TNM) (https://www.uicc.org/resources/tnm), which is a globally recognised standard for classifying the extent of spread of cancer. Each individual aspect of TNM is termed as a category: T category describes the primary tumor site and size; N category describes the regional lymph node involvement; M category describes the presence or otherwise of distant metastatic spread. Gene enrichment analysis (GSEA) was used to understand the differentially functional enrichment pathways of samples with different types (Tang et al., 2021). In order to explore the characteristics of immune microenvironment of different types of samples, the R package “limma” was used to filter the differences in immune cell infiltration of different types of samples (Li Z. et al., 2021). In addition, we further analyzed the expression of 29 immune checkpoint moleclues in different types of samples, including CD70, TNFSF9, FGL1, CD276, NT5E, HHLA2, VTCN1, TNFRSF18, CD274, PDCDL1LG2, IDO1, CTLA4, ICOS, HAVCR, PDCD1, LAG3, SIGLEC15, TNFSF4, TNFRSF9, ENTPD1, VSIR, ICOSLG, TNFSF14, TMIGD2, CD27, NCR3, BTLA, CD40LG, CD40, TNFRSF4 (Zhang C. et al., 2020). Moreover, we also analyzed the differences in the expression levels of the tumor suppressor gene TP53 among the subtypes.
2.2.3. A Prognostic Signature Constructed by m6A/m5C/m1A-Related lncRNAs
The least absolute shrinkage and selection operator (LASSO) cox regression was used to select the most significant prognostic lncRNAs in 44 m6A/m5C/m1A-related lncRNAs for predictive model construction (Wang S. et al., 2020). Lasso is a linear regression method based on L1-regularization. L1-regularization can reduce the complexity of the model and the risk of over-fitting of the model. The explanation in regression analysis is that if two characteristic variables are highly correlated, the regression analysis without regularization may give them approximately equal weight, so the weighted result is still small, but due to the large weight, the over-fitting problem is caused. L1-regularization greatly increases the stability of modeling and achieves the purpose of selecting characteristic variables. The most outstanding advantage of Lasso lies in the penalized regression analysis of all variable coefficients, the relatively unimportant independent variable coefficients become 0, which is excluded from modeling (Vaid et al., 2021; Yu et al., 2021; Hoq et al., 2021). However, the 44 lncRNAs are too large for constructing a prognosis model, so we used LASSO regression analysis to further compress the 44 lncRNAs to get the lncRNAs with the most prognostic value, and their correlation coefficients. Combining respective correlation coefficients of obtained m6A/m5C/m1A-related lncRNAs, the respective gene expression level of these lncRNAs in each tumor sample, and the risk score calculating formula (The risk score calculating formula is:
Considering the strength expression of markers in samples. If it is a weak expression, It might be hard to be detected in clinical diagnosis. Therefore, we searched these filtered lncRNAs included in this study in TIMER: Tumor IMmune Estimation Resource (http://timer.cistrome.org/) (Zhu Y. et al., 2021).
2.2.4. Immune Landscapes Assessment of the Prognostic Signature
In order to explore the immune landscapes of the high-risk group and low-risk group, we analyzed the correlation between risk score and immune cells by Spearman analysis. The R Package “Limma” was used to filter the immune cell infiltration of samples with different risk scores. In addition, we analyzed the expression levels of 29 immune checkpoint moleclues in samples from different risk groups, including CD70, TNFSF9, FGL1, CD276, NT5E, HHLA2, VTCN1, TNFRSF18, CD274, PDCDL1LG2, IDO1, CTLA4, ICOS, HAVCR, PDCD1, LAG3, SIGLEC15, TNFSF4, TNFRSF9, ENTPD1, VSIR, ICOSLG, TNFSF14, TMIGD2, CD27, NCR3, BTLA, CD40LG, CD40, TNFRSF4.
2.2.5. Simple Nucleotide Variations of the Prognostic Signature
Based on the tumor mutational burden (TMB) information of each sample (Bonfield et al., 2021), combined with the survival information analysis of patients, all samples were divided into the high-TMB group and the low-TMB group, and the survival of the two groups was analyzed. The survival of the high-TMB + high-risk group, the high-TMB + low-risk group, the low-TMB + high-risk group and the low-TMB + low-risk score group were jointly analyzed to observe the influence of the risk score and TMB on the survival of the patients. Genes with significant mutation rates in the high-risk group and low-risk group were further analyzed to find the ones with the most significant mutation rates for further analysis.
2.2.6. Prognostic Analysis and Immune Infiltration of Single lncRNA
Based on the expression levels of each independent lncRNA in the 6 m6A/m5C/m1A-related lncRNAs Prognostic signatures, the patient samples were divided into high expression group and low expression group. The R package “survival” was used to obtain survival curves for each independent lncRNA. There are six types of immune infiltration in human cancers, including Wound Healing (Immune C1), IFN-gamma Dominant (Immune C2), Inflammatory (Immune C3), Lymphocyte Depleted (Immune C4), Immunologically Quiet (Immune C5) and TGF-beta Dominant (Immune C6) (Soldevilla et al., 2019; Zhao et al., 2021). Based on the relevant files downloaded from the UCSC Xena database, the R package “CIBERSORT” and “ESTIMATE” (Wang S. et al., 2020) were used to obtain relevant data (Stromal Score, Immune Score, ESTIMAT Score) of tumor purity. Furthermore, the correlation between the expression levels of 6 m6A/m5C/m1A-related lncRNAs and tumor purity in the immune microenvironment of HNSCC was further analyzed.
2.2.7. Drug Sensitivity Analysis
The R package “Prophetic” was used to analyze the drug sensitivity analysis of samples from different risk groups (Ding et al., 2021). Combined with the relevant data downloaded from the CellMiner Database, the correlation analysis of drug sensitivity and the expression level of 6 m6A/m5C/m1A-related lncRNAs was further carried out (Xu X. et al., 2020).
3.1. m6A/m5C/m1A-Related lncRNAs
The flow chart of our work was shown in Figure 1.
Using the downloaded file from the “GENCODE” website, we obtained 14,086 lncRNAs expression matrix. In addition, 52 expression matrix of m6A/m5C/m1A-related genes were screened out from the expression matrix of mRNAs. Through co-expression analysis, we screened out 147 lncRNAs co-expressed with m6A/m5C/m1A-related genes (Figure 2A). We further screened out 44 lncRNAs with potential prognostic value from 147 co-expressed lncRNAs. The forest plot showed that all 44 lncRNAs were protective factors (hazard ratio < 1) in patients with HNSCC (Figure 2B). The heatmap and box figures show that the expressions of 44 lncRNAs in normal samples and tumor samples are significantly different (Figures 2C,D).
FIGURE 2. m6A/m5C/m1A-related mRNA–lncRNA co-expression network and identification of prognostic value of m6A/m5C/m1A-related lncRNAs. (A) 147 lncRNAs co-expressed with m6A/m5C/m1A-related genes. Green represents lncRNAs, while red represents m6A/m5C/m1A-related genes. (B) The forest plot of univariate Cox regression showed that 44 lncRNAs with prognostic value were protective factors of HNSCC. (C,D) The heatmap and boxplot of differential expression of 44 lncRNAs with prognostic value between 44 normal tissues and 501 tumor tissues in TCGA-HNSCC database.
3.2. Consensus Clustering of m6A/m5C/m1A-Related lncRNAs Identified Two Clusters of HNSCC Patients
By cluster analysis, the tumor samples in the TCGA database were classified into two types, consensus matrix for optimal K = 2 (Figure 3A). PCA was used to represent the characteristics of the typing samples, and it was found that the samples could be visually distinguished (Figure 3B). The Kaplan-Meier survival curves showed that the prognosis of patients in Cluster 2 was significantly better than that of patients in Cluster 1 (p = 0.03) (Figure 3C), and the two clinical characteristics of patients’ HPV status and tumor grade were significantly different among different types (p < 0.001) (Figure 3D). GSEA analysis showed that the pathways of Cluster 2 were mainly enriched in CELL_CYCLE, LYSINE_DEGRADATION, P53_SIGNALING_PATHWAY, et al. The functional enrichment pathways of Cluster1 were mainly enriched in ECM_receptor_interaction, FOCAL_ADHESION, OCAL_ADHESION, et al. (Figure 3E). According to the immune microenvironments of the two subtypes, we can found the infiltration of T cells follicular helper and T cells regulatory in Cluster 1 significantly less than Cluster 2 (p = 0.005; p = 0.009), while Macrophages M0 infiltrate significantly more than Cluster 2 (p = 0.05) (Figure 3F). We found that the expression levels of 9 immune checkpoint genes were significantly different between the two subtypes (p < 0.05). The expression levels of NT5E, ICOS, PDCD1LG2, HAVCR2, TNFSF4, TNFRSF9, EntPD1 and TNFSF14 were higher in Cluster 1, and TNFRSF18 was higher in Cluster 2 (Figure 3G). In addition, the expression level of tumor suppressor gene TP53 in Cluster 2 was significantly higher than that in Cluster 1 (p < 0.001) (Figure 3H).
FIGURE 3. Clinical characteristics and overall survival among different subgroups of HNSCC. *, **, *** represent p < 0.05, p < 0.01 and p < 0.0001, respectively. (A) Consensus matrix for optimal k = 2. (B) Principal component analysis (PCA) of cluster 1 and cluster 2. (C) Kaplan-Meier curve of overall survival (OS) time in Cluster 1 and Cluster 2 (p = 0.03). (D) Heatmap of the clinical features of the two clusters. (E) GSEA analysis displayed key pathways in in Cluster 1 and Cluster 2. (F) The violet plot showed the differential expression of immune cells between Cluster 1 and Cluster 2, in which Cluster 1 was represented in blue and Cluster 2 was represented in red. (G) The heatmap of the differential expression levels of 29 immune checkpoint genes between Cluster 1 and Cluster 2. (H) The boxplot of TP53 expression level between Cluster 1 and Cluster 2.
3.3. The Risk Prognostic Signature Built by m6A/m5C/m1A-Related lncRNAs
We got 44 lncRNAs with prognostic value, but it is too difficult to use all lncRNAs to build a prognostic model, so we compressed 44 lncRNAs with Lasso analysis method, and finally got 6 lncRNAs with the most prognostic value to build our prognostic model: AL035587.1, AC009121.3, AF131215.5, FMR1-IT1, AC106820.5, PTOV1-AS2 (Figures 4A,B). The correlation coefficient was shown in Table 2.
FIGURE 4. The m6A/m5C/m1A-related lncRNAs Prognostic signature. (A,B) 6 m6A/m5C/m1A-related lncRNAs were found by using least absolute shrinkage and selection operator (LASSO) cox regression. (C) Kaplan-Meier curve of overall survival of high-risk group and low-risk group. (D,E) Kaplan-Meier curve of overall survival of train group and test group. (F,G) The distribution of risk score, survival status and m6A/m5C/m1A-related lncRNAs expression panal.
The survival time and survival rate of patients in the high-risk subgroup were significantly lower than those in the low-risk subgroup (p = 5.262e-04) (Figure 4C), which was also verified in samples from the training set and the testing set (Train group: p = 0.008; Test group: p = 0.02) (Figures 4D,E). It was found in the training set and testing set that the prognosis of patients in the high-risk subgroup was worse than that in the low-risk subgroup, and the expression levels of 6 lncRNAs in the low-risk subgroup were significantly higher than that in the high-risk subgroup (Figures 4F,G). The univariate and multivariate regression analysis results showed that risk score, as an independent prognostic indicator, had a certain predictive value relative to individual and all clinical characteristics (age, gender, tumor grade, tumor stage) in both the training set and the testing set (p < 0.01) (Supplementary Figures S1A–D). The results of PCA and T-SNE show that the use of risk score can intuitively distinguish the high-risk group from the low-risk group in the training set and the testing set (Supplementary Figures S1E–H). As shown in Supplementary Figure S2, the area under curve (AUC) at 5 years of overall survival for the LncSig (our study) is 0.607, which was significantly higher than that of JianglncSig (AUC = 0.602) and WanglncSig (AUC = 0.589). These results indicate that LncSig has better prognostic performance in predicting survival than the two recently published lncRNA signatures (Wang et al., 2018; Jiang et al., 2020).
In combination with clinical features, analysis showed that patients in the tumor grade G3-4, N0, N1-3, stage I–II, stage III–IV, T1-2, T3-4, the survival of high-risk subgroup was significantly worse than low-risk subgroup (p < 0.05) (Figures 5A–H). The heatmap showed the expression of 6 selected lncRNAs in patients with high and low risk. We observed the statistically significant differences between the high-risk and low-risk groups in gender, HPV-status, tumor stage and tumor grade (p < 0.05) (Figure 5I). The results of sankey diagram showed that most of the patients in Cluster 2 belonged to the low-risk group and had a higher survival rate than Cluster 1, which also supported the results of survival status analysis of the above classifications (Figure 5J). GO analysis showed that cell functions enriched by different genes between the high-risk group and the low-risk group were as follows: B cell mediated immunity, complement activation, complement activation classical pathway, humoral immune response mediated by circulating immunoglobulin, immunoglobulin mediated immune response (Figure 5K). KEGG analysis showed that the signal pathways enriched by differential genes between the high-risk group and the low-risk group were: drug metabolism-cytochrome P450, estrogen signaling pathway, leukocyte transendothelial migration, retinol metabolism, signaling pathway regulating pluripotency of stem cells (Figure 5L).
FIGURE 5. Correlation of risk model subgroups with overall survival and different clinical characteristics in TCGA-HNSCC database. (A–H) Kaplan–Meier curve of the overall survival of patients with different clinical traits in high-risk group and low-risk group. (I) Heatmap of association between expression levels of six prognosis lncRNAs and clinicopathological features in TCGA database. (J) The sankey diagram of the relationship among the cluster, risk score and survival state. (K,L) GO and KEGG analysis.
The gene expression box plots (Supplementary Figure S3) of FMR1-IT1 and PTOV1-AS2 were found in the Module of exploring associations between gene expression and tumor features in TCGA from TIMER, which showed that the expression of FMR1-IT1 and PTOV1-AS2 in HSNCC tissues is significantly higher than that in normal tissues.
3.4. Immune Landscapes and Tumor Mutational Burden (TMB) of the Prognostic Signature
We found that 12 immune checkpoint genes were significantly different between the two risk groups (p < 0.05). CD70, ICOS, HAVCR2, PDCD1, SIGLEC15, TNFSF4, TNFRSF9, ENTPD1, TNFSF14, CD27, BTLA, CD40LG. The expression of BTLA was higher in low-risk group. CD70, ICOS, HAVCR2, PDCD1, SIGLEC15, TNFSF4, TNFRSF9, ENTPD1, TNFSF14, CD27 and CD40LG were higher expressed in high-risk group (Figure 6A). Among them, the overlapping immune checkpoint genes ICOS, HAVCR2, TNFSF4, TNFRSF9, ENTPD1 and TNFSF14 were all increased in Cluster 1 and high-risk subgroups.
FIGURE 6. Differences in immune microenvironment between high-risk and low-risk subgroups. (A) Violet plot of the expression of different checkpoint genes between high-risk subgroup and low-risk subgroup. (B) Heatmap of infiltration of different immune cells between high-risk subgroup and low-risk subgroup. (C–L) The correlation between immune cell infiltration and risk score.
Immunity microenvironment analysis showed that the contents of NK cell resting, macrophages M2, and neutrophils in samples of low-risk group were significantly lower than those of high-risk group (p < 0.05), while the contents of B cells navie, plasma cells, and T cells regulatory (Tregs) were on the contrary (p < 0.05) (Figure 6B). The correlation analysis between risk score and immune cell infiltration showed that neutrophils, dendritic cells activated, NK cells resting, dendritic cells resting, macrophages M2, mast cells activated were positively correlated with risk score. B cells naive, plasma cells, T cells helper and T cells regulatory (Tregs) were negatively correlated with risk score (p < 0.05) (Figures 6C–L).
TMB analysis showed that the survival of the high-TMB group was significantly lower than that of the low-TMB group (p = 0.002) (Figure 7A). Combined analysis of risk score showed that patients in the low-TMB + low-risk score group had a significantly better prognosis than those in the other three groups (p < 0.001), suggesting that patients with low TMB tend to have a better prognosis (Figure 7B). Looking for the significant mutant genes in the high-risk group and the low-risk group, it was found that the mutation rate of TP53 was the highest in both groups, and the expression of TP53 in the high-risk group was significantly lower than that in the low-risk group (p = 1.3e-07) (Figures 7C–E). The combined analysis of TP53 and risk score showed that the expression level of TP53 in the high-risk group was higher than that in the low-risk group (p = 1.3e-07). Patients in the low-risk subgroup of TP53 wild-type had significantly better outcomes than those in the other three groups (p = 0.002) (Figure 7F).
FIGURE 7. The impact of differential TMB level on patients’ prognosis and TP53 somatic mutation. (A) Comparison of prognosis between high TMB subgroup and low TMB subgroup. (B) Comparison of prognosis of TMB combined with risk score. (C,D) Top 20 genes with the most significant mutation rates in high-risk subgroup and low-risk subgroup. (E) The differential expression level TP53 between high-risk group and low-risk group. (F) Kaplan–Meier curve analysis of overall survival was shown for patients classified according to TP53 mutation status and risk score. TP53-wt: wild type of TP53 sequence; TP53-mt: mutation type of TP53 sequence.
3.5. Prognostic Analysis and Immune Infiltration of m6A/m5C/m1A-Related lncRNAs
The Kaplan-Meier survival curves showed that there was no significant difference between high-expression and low-expression of AF131215.5. However, the high-expression of AL035587.1, AC009121.3, FMR1-IT1, AC106820.5, PTOV1-AS2 had the better overall survival than low-expression in TCGA database (p < 0.05) (Figures 8A–F).
FIGURE 8. Kaplan-Meier curves of overall survival for differential expression level of 6 m6A/m5C/m1A-related lncRNAs Prognostic signatures respectively. (A) The association of AL035587.1 expression level with patients’ prognosis of TCGA database. (B) The association of AC009121.3 expression level with patients’ prognosis of TCGA database. (C) The association of AF131215.5 expression level with patients’ prognosis of TCGA database. (D) The association of FMR1-IT1 expression level with patients’ prognosis of TCGA database. (E) The association of AC106820.5 expression level with patients’ prognosis of TCGA database. (F) The association of PTOV1-AS2 expression level with patients’ prognosis of TCGA database.
According to the results of the immunoinfiltration analysis, the high expression of AL035587.1 was considered to be related to the type of C3 and C6 infiltration according to the results of Figure 9A. The high expression of AC106820.5 was considered to be correlated with the type of C3 and C6 infiltration (p < 0.01). Moreover, the tumor microenvironment (TME) was assessed by estimation analysis. The expression levels of PTOV1-AS2, AC009121.3, FMR1-IT1 and AC106820.5 were negatively correlated with the stromalscore (Figure 9B). In addition, AL035587.1 and AF131215.5 were positively correlated with immunescore, while AC106820.5 was negatively correlated. Estimatescore (combining stromalscore and immunescore) showed that the expression levels of AC009121.3, AC106820.5 and FMR1-IT1 were negatively correlated with estimatescore, while AL035587.1 and AF131215.5 were positively correlated (Figure 9B) (p < 0.05).
FIGURE 9. Association of expression of prognostic lncRNAs with with immune response, TME. *, **, *** represent p < 0.05, p < 0.01 and p < 0.0001, respectively. (A) Association of prognostic lncRNAs expression with immune subtypes in HNSCC patients. (B) Correlation of expression of prognostic lncRNAs with TME (Stromal score, Immune score, and ESTIMATE score).
3.6. Drug Sensitivity
Drug sensitivity analysis showed that samples from the low-risk group had higher sensitivity to Paclitaxel, Docetaxel and Gefitinib and lower sensitivity to Methotrexate than those from the high-risk group (Figures 10A–D) (p < 0.001). Two m6A/m5C/m1A-related lncRNAs: PTOV1-AS2 and FMR1-IT1 were extracted from the CellMiner Database NCI 60 RNA seq and Compound activity: DTP NCI-60. It was found that the sensitivity of drug Nelarabine had the highest correlation with the expression level of PTOV1-AS2 (Correlation = 0.418, p < 0.001), and the sensitivity of Nelarabine was also the highest associated with FMR1-IT1 expression (Correlation = 0.337, p = 0.008) (Figure 10E).
FIGURE 10. Drug sensitivity analysis. (A–D) The boxplot of Sensitivity of Paclitaxel, Docetaxel, Gefitinib, Methotrexate in patients with high-risk subgroup and low-risk subgroup. (E) The scatter plot showed the correlation between prognostic lncRNAs expression and drug sensitivity (Z score of CellMiner interface). Top 10 associations were displayed, sorted by p value.
Studies have shown that dynamic RNA methylation, the modification events (such as m6A, m5C, and m1A) are involved in the progression of tumors, including promoting cancer cell proliferation or regulating invasiveness and metastatic potential (Barbieri and Kouzarides, 2020). m6A can also be involved in tumor proliferation, migration, drug resistance, or as a prognostic marker of tumor (Han et al., 2019; Wang M. et al., 2020; Wang H. et al., 2021; Paramasivam and Priyadharsini, 2021; Xu et al., 2021). Some studies have shown that there is an interaction between RNA methylation and lncRNA in tumors, and m6A reader YTHDF2-mediated degradation of lncRNA FENDRR promoted cell proliferation in endometrioid endometrial carcinoma (Shen et al., 2021). m6A-mediated upregulation of LINC00857 promoted pancreatic cancer tumorigenesis (Meng L. et al., 2021). LncRNA KB-1980E6.3 maintained breast cancer stem cell stemness via interacting with IGF2BP1 (Zhu P. et al., 2021). Abnormal m5C modification of lncRNA H19 mediated by NSUN2 was associated with poor differentiation of hepatocellular carcinoma (Sun et al., 2020). Conversely, lncRNAs can also affect the function of genes related to RNA methylation modifications. LncRNA MALAT1 contributed to thyroid cancer progression through the upregulation of IGF2BP2 by binding to miR-204 (Ye et al., 2021). LncRNA LINRIS can stabilize IGF2BP2 and promote aerobic glycolysis in colorectal cancer (Wang et al., 2019). The above evidences indicate the interaction between RNA methylation modifications and lncRNAs in the tumorigenesis of a variety of cancers. However, the evidence of how the RNA methylation modification involved in the development of HNSCC by influencing the lncRNA is still insufficient. Therefore, this study aims to explore whether the relevant lncRNAs are involved in the development of the HNSCC.
We downloaded the data of 44 Normal samples and 501 tumor samples from the TCGA database to search for m6A/m5C/m1A-related lncRNAs and we found 44 lncRNAs with prognostic value. Tumor samples were divided into Cluster 1 and Cluster 2 according to the expression levels of these 44 lncRNAs in all tumor samples, and the survival of Cluster 2 was significantly better than that of Cluster 1. We further used the lasso cox regression to construct the prognostic signature using 6 m6A/m5C/m1A-related lncRNAs to predict survival in patients with HNSCC. Based on the correlation coefficients of the 6 m6A/m5C/m1A-related lncRNAs and their expression levels in each tumor sample, we obtained risk scores for each tumor sample and constructed a prognostic signature. According to Figure 5J the sankey diagram of the relationship among the cluster, risk score and survival state, it can be found that most patients in high risk group come from cluster 1, while most patients in cluster 2 are assigned to low risk group. The survival time of patients in the high risk group was significantly shorter than that in the low risk group (p < 0.05), which also proved that both the clustering type composed of 44 lncRNAs and the prognosis model composed of 6 lncRNAs can well predict the prognosis of tumor patients. Differences in survival were also demonstrated in the training and validation sets of randomly assigned tumor samples. The univariate and multivariate cox regression analysis showed that risk score could be used as a prognostic marker for HNSCC independent of age, sex, tumor grade, and tumor stage. This suggested that the prognostic signature constructed by 6 lncRNAs had certain prognostic value. These results supply a reference for further functions and mechanisms of prognostic signatures in HNSCC.
Six m6A/m5C/m1A-related lncRNAs were obtained, which were AL035587.1, AC009121.3, AF131215.5, FMR1-IT1, AC106820.5 and PTOV1-AS2. The results in Supplementary Figure S3 showed that FMR1-IT1 and PTOV1-AS2 can be detected in clinical diagnosis. AF131215.5 is one of lncRNAs in lung adenocarcinoma (LUAD)-related competing endogenous RNA networks, and showed an independent prognostic value of overall survival for patients with LUAD (Hou and Yao, 2021). FMR1 intronic transcript 1 (FMR1-IT1) is also known as FMR1 intronic transcript 1 (non-protein coding) or Fragile X Mental Retardation Syndrome Protein (FMRP) or FMR1, there were evidence that FMRP regulated tumor invasiveness-related pathways in melanoma cells by impacting cell migration, invasion and adhesion (Zalfa et al., 2017), and its overexpression was associated with lung metastasis of murine breast cancer (Lucá et al., 2013). However, there is still not enough evidence to clarify the relationship between tumorigenesis and prognosis of these 6 lncRNAs and HNSCC. We hope that our work can explain the role of these 6 lncRNAs in the carcinogenesis and development of HNSCC and clarify their prognostic value.
Based on the expression levels of these 6 m6A/m5C/m1A-related lncRNAs in the samples, we further studied the tumor microenvironment of the samples and analyzed the expression levels of 29 immune checkpoint genes. There were significant differences in the expression of various immune cells and immune checkpoint genes between the high-risk subgroup and the low-risk subgroup. Among them, the overlapping immune checkpoint genes of Cluster 1 and high-risk subgroups were ICOS, HAVCR2, TNFSF4, TNFRSF9, ENTPD1 and TNFSF14. Their expression levels were all elevated, revealing the characteristics of tumor immune microenvironment between different subgroups. Currently, only a small number of HNSCC patients benefit from immunotherapy, and the need to discover new biomarkers to optimize treatment strategies has become increasingly important in clinical practice. With the rise of the application of immune checkpoint inhibitor in tumor treatment, such as Programmed Cell death-1 (PD-1) inhibitors (Kitamura et al., 2020), there are many ongoing trials in HNSCC that focus on identifying new biomarkers (Gavrielatou et al., 2020). As a member of the B7 superfamily, ICOS is a co-stimulatory receptor for T-cell enhancement. In preclinical studies, ICOS agonist monoclonal antibodies have been shown to enhance inhibitory checkpoint blockade. On the contrary, the antagonistic monoclonal antibody against ICOS can not only inhibit the lymphoid tumor cells expressing ICOS, but also inhibit the immunosuppressive Treg (Amatore et al., 2018). ENTPD1 (CD39) converts extracellular ATP into adenosine, thereby inhibiting T cell effector functions through the adenosine receptor A2A. The CD39 inhibitor POM-1 inhibits adenosine production and reduces T cell suppression. The adenosine pathway can be targeted to multiple myeloma, and blocking this pathway can replace PD1/PDL1 to inhibit multiple myeloma and other blood system cancers (Yang et al., 2020). At present, TNFSF14-based therapy has shown great efficacy in cancer immunotherapy, which can modify tumor microenvironment by normalizing tumor vessels and significantly improve the infiltration of effector infiltrating lymphocytes, thus reducing tumor burden and forming lasting anti-tumor memory (Skeate et al., 2020). These evidences also suggest the possibility of selecting differential immune checkpoint genes as therapeutic targets for HNSCC. In addition, we have analyzed the joint tumor mutation burden and risk score model. It was found that the prognosis of patients with high burden of tumor mutations was worse than that of patients with low burden, while the prognosis of patients with high burden of mutations in the high-risk subgroup was the worst, and the prognosis of patients with low burden of mutations in the low-risk subgroup was the best. This undoubtedly indicated the importance of tumor mutation burden in predicting prognosis. Moreover, we observed that the tumor suppressor gene TP53 was highly expressed in Cluster 2 and the low-risk subgroup. It is generally believed that TP53 gene can be divided into wild type and mutant type. The wild-type TP53 gene, namely tumor suppressor gene, has a control and negative regulation effect on cell proliferation. However, the mutant TP53 gene lost its tumor suppressor effect and promoted the transformation of cells to malignancy (Duffy et al., 2020; Zhu et al., 2020; Meng X. et al., 2021). Our results showed that TP53 mutant patients in the high-risk subgroup had the worst prognosis, and TP53 wild-type patients in the low-risk subgroup had the best prognosis, which also indicated that higher mutation rate of TP53 was associated with poorer prognosis. In addition, some studies have revealed the connection between lncRNA and TP53. Huarte, M, et al. found that LincRNA-p21 is the first lncRNAs identified as being induced by wt-p53 transcription. It acts as a transcriptional repressor in the p53 pathway (Huarte et al., 2010). Chang et al. (2011) found that lncRNA XIST and the p53-inducible lncRNA NEAT1 act as ceRNA for miR-34a, which leads to the activation of epithelial-mesenchymal transition. The research of Olivero et al. (2020) showed that TP53 activates the lncRNA pvt1b to inhibit myelocytomasis and suppress tumorigenesis. It is worth noting that targeting mut-TP53/ncRNA signals may have a great potential impact on the treatment of cancer patients. Interestingly, lncRNA and immune checkpoint are also closely related, lncRNAs are involved in the regulation of various cellular pathways and immune checkpoints (Di Martino et al., 2021). Tang et al. (2017) found that lncRNA AFAP1-AS1 was co-expressed with PD-1 lymphocytes associated with nasopharyngeal carcinoma, which resulted in invasive and poor prognosis phenotype. In hepatocellular carcinoma, lncRNA Tim3 mediates T cell exhaustion and inhibits T cell immune response, which leads to tumor immunosuppression through TIM-3 specific binding (Ji et al., 2018). In the future, the potential combination of lncRNAs targeted drugs and tumor immune checkpoint inhibitor antibodies may be used for experimental treatment of human cancer. Considering that the six m6A/m5C/m1A-related lncRNAs were all protective factors of HNSCC (HR < 1), we analyzed the survival of the high-risk and low-risk subgroups of these 6 lncRNAs, and found that the prognosis of the high-expression subgroup of 5 lncRNAs was significantly better than that of the low-expression subgroup, which also proved their role as protective factors. The high expression of AL035587.1 and AC106820.5 is thought to correlate with the type of Immune C3 and C6 invasion, indicating their tumor suppressive effect. Estimatedscores (a combination of stromalscore and immunescore) showed that high expression levels of AC009121.3, AC106820.5 and FMR1-IT1 were significantly associated with lower tumor purity, while AL035587.1 and AF131215.5 had the opposite effect (Sawatani et al., 2020).
Drug sensitivity analysis showed that Paclitaxel, Docetaxel, and Gefitinib had higher sensitivity in the low-risk group than in the high-risk group, while Methotrexate had higher sensitivity in the high-risk group. These four drugs are commonly used in the clinical treatment of HNSCC (Yin et al., 2019; Ham et al., 2020; Guigay et al., 2021), and our results also demonstrate their value in the application of HNSCC. For the therapeutic potential of these 6 lncRNAs, we analyzed their sensitivity to different small molecule drugs. The results showed that PTOV1-AS2 was sensitive to four small molecular agents (Nelarabine, Decitabine, Acrichine, Dasatinib) (p < 0.05), and FMR1-IT1 was sensitive to five small molecular agents (Nelarabine, Chelerythrine, Tamoxifen, Tanespimycin, Sunitinib (p < 0.05). These eight drugs have certain anticancer activity and can be considered as candidates for the targeted therapy of HNSCC. Nelarabine is a purine analogue approved for the treatment of patients with T-cell lymphoblastic lymphoma (Teachey and O’Connor, 2020). Decitabine is a DNA demethylation agent, which plays an unmethylated role in inducing apoptosis and activating autophagy in myeloid leukemia (Li L. et al., 2021). Dasatinib is an antagonist of the SRC family of kinases, and may potentially induce T-cell response leading to improved survival of malignant pleural mesothelioma (Chen et al., 2020). Chelerythrine, Tamoxifen, Tanespimycin and Sunitinib have also been shown to have anticancer effects in a variety of cancers (Heng and Cheah, 2020; Wang and Liu, 2020; Ocadlikova et al., 2021; Su et al., 2021). These evidences also indicate the prospect of targeting these two lncRNAs for the treatment of HNSCC.
The poor prognosis of the HNSCC affects the health of tens of millions of people every year, and improving the prognosis of the HNSCC is particularly important. Studies have shown that m6A/m5C/m1A RNA methylation modifications was involved in the progression of cancer. Currently, many studies have chosen to use lncRNAs to establish prognostic markers, in order to provide new targets for diagnosis and treatment of cancers. Our study elucidates how m6A/m5C/m1A -related lncRNAs are related to the prognosis, immune microenvironment, and tumor mutation load of HNSCC. The prognostic signature can be used as an independent factor to better predict the prognosis of HNSCC, which may potentially become a new option for immunotherapy of HNSCC in the future.
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.
Conceptualization, EW, YL, SZ, HX; Data Collection, EW, YL, RM, JW, PD, PZ; Methodology, EW, YL, SZ, HX, RM, JW, PD, PZ; Writing–original draft, EW, YL, SZ, RM, JW, PD, PZ; Writing–review and editing, EW, YL, SZ, HX. All authors have read and agreed to the published version of the manuscript.
This research was supported by the Natural Science Foundation of China (81771002, 82071057).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.718974/full#supplementary-material
Amatore, F., Gorvel, L., and Olive, D. (2018). Inducible Co-stimulator (ICOS) as a Potential Therapeutic Target for Anti-cancer Therapy. Expert Opin. Ther. Targets 22 (4), 343–351. doi:10.1080/14728222.2018.1444753
Bonfield, J. K., Marshall, J., Danecek, P., Li, H., Ohan, V., Whitwham, A., et al. (2021). HTSlib: C Library for reading/writing High-Throughput Sequencing Data. Gigascience 10 (2), giab007. doi:10.1093/gigascience/giab007
Castilho, R., Squarize, C., and Almeida, L. (2017). Epigenetic Modifications and Head and Neck Cancer: Implications for Tumor Progression and Resistance to Therapy. Ijms 18 (7), 1506. doi:10.3390/ijms18071506
Chang, C.-J., Chao, C.-H., Xia, W., Yang, J.-Y., Xiong, Y., Li, C.-W., et al. (2011). p53 Regulates Epithelial-Mesenchymal Transition and Stem Cell Properties through Modulating miRNAs. Nat. Cel Biol 13 (3), 317–323. doi:10.1038/ncb2173
Chen, R., Lee, W.-C., Fujimoto, J., Li, J., Hu, X., Mehran, R., et al. (2020). Evolution of Genomic and T-Cell Repertoire Heterogeneity of Malignant Pleural Mesothelioma under Dasatinib Treatment. Clin. Cancer Res. 26 (20), 5477–5486. doi:10.1158/1078-0432.Ccr-20-1767
Cossu, A. M., Mosca, L., Zappavigna, S., Misso, G., Bocchetti, M., De Micco, F., et al. (2019). Long Non-coding RNAs as Important Biomarkers in Laryngeal Cancer and Other Head and Neck Tumours. Ijms 20 (14), 3444. doi:10.3390/ijms20143444
Dai, F., Wu, Y., Lu, Y., An, C., Zheng, X., Dai, L., et al. (2020). Crosstalk between RNA m6A Modification and Non-coding RNA Contributes to Cancer Growth and Progression. Mol. Ther. - Nucleic Acids 22, 62–71. doi:10.1016/j.omtn.2020.08.004
Di Martino, M. T., Riillo, C., Scionti, F., Grillone, K., Polerà, N., Caracciolo, D., et al. (2021). miRNAs and lncRNAs as Novel Therapeutic Targets to Improve Cancer Immunotherapy. Cancers 13 (7), 1587. doi:10.3390/cancers13071587
Ding, C., Shan, Z., Li, M., Chen, H., Li, X., and Jin, Z. (2021). Characterization of the Fatty Acid Metabolism in Colorectal Cancer to Guide Clinical Therapy. Mol. Ther. - Oncolytics 20, 532–544. doi:10.1016/j.omto.2021.02.010
Frankish, A., Diekhans, M., Ferreira, A.-M., Johnson, R., Jungreis, I., Loveland, J., et al. (2019). GENCODE Reference Annotation for the Human and Mouse Genomes. Nucleic Acids Res. 47 (D1), D766–d773. doi:10.1093/nar/gky955
Gavrielatou, N., Doumas, S., Economopoulou, P., Foukas, P. G., and Psyrri, A. (2020). Biomarkers for Immunotherapy Response in Head and Neck Cancer. Cancer Treat. Rev. 84, 101977. doi:10.1016/j.ctrv.2020.101977
Gong, Z., Chen, J., Cheng, J.-N., Sun, C., Jia, Q., Diao, X., et al. (2018). Tumor Microenvironment Properties Are Associated with Low CD68-Positive Cell Infiltration and Favorable Disease-free Survival in EGFR-Mutant Lung Adenocarcinoma. Clin. Lung Cancer 19 (5), e551–e558. doi:10.1016/j.cllc.2018.03.011
Guigay, J., Aupérin, A., Fayette, J., Saada-Bouzid, E., Lafond, C., Taberna, M., et al. (2021). Cetuximab, Docetaxel, and Cisplatin versus Platinum, Fluorouracil, and Cetuximab as First-Line Treatment in Patients with Recurrent or Metastatic Head and Neck Squamous-Cell Carcinoma (GORTEC 2014-01 TPExtreme): a Multicentre, Open-Label, Randomised, Phase 2 Trial. Lancet Oncol. 22, 463–475. doi:10.1016/s1470-2045(20)30755-5
Ham, J. C., Meerten, E., Fiets, W. E., Beerepoot, L. V., Jeurissen, F. J. F., Slingerland, M., et al. (2020). Methotrexate Plus or Minus Cetuximab as First‐line Treatment in a Recurrent or Metastatic (R/M) Squamous Cell Carcinoma Population of the Head and Neck (SCCHN), Unfit for Cisplatin Combination Treatment, a Phase Ib‐randomized Phase II Study Commence. Head & Neck 42 (5), 828–838. doi:10.1002/hed.26053
Han, J., Wang, J.-z., Yang, X., Yu, H., Zhou, R., Lu, H.-C., et al. (2019). METTL3 Promote Tumor Proliferation of Bladder Cancer by Accelerating Pri-miR221/222 Maturation in m6A-dependent Manner. Mol. Cancer 18 (1), 110. doi:10.1186/s12943-019-1036-9
Heng, W. S., and Cheah, S.-C. (2020). Chelerythrine Chloride Downregulates β-Catenin and Inhibits Stem Cell Properties of Non-small Cell Lung Carcinoma. Molecules 25 (1), 224. doi:10.3390/molecules25010224
Huang, J., Chen, Z., Chen, X., Chen, J., Cheng, Z., and Wang, Z. (2021). The Role of RNA N6-Methyladenosine Methyltransferase in Cancers. Mol. Ther. - Nucleic Acids 23, 887–896. doi:10.1016/j.omtn.2020.12.021
Huarte, M., Guttman, M., Feldser, D., Garber, M., Koziol, M. J., Kenzelmann-Broz, D., et al. (2010). A Large Intergenic Noncoding RNA Induced by P53 Mediates Global Gene Repression in the P53 Response. Cell 142 (3), 409–419. doi:10.1016/j.cell.2010.06.040
Ji, J., Yin, Y., Ju, H., Xu, X., Liu, W., Fu, Q., et al. (2018). Long Non-coding RNA Lnc-Tim3 Exacerbates CD8 T Cell Exhaustion via Binding to Tim-3 and Inducing Nuclear Translocation of Bat3 in HCC. Cell Death Dis 9 (5), 478. doi:10.1038/s41419-018-0528-7
Jiang, H., Ma, B., Xu, W., Luo, Y., Wang, X., Wen, S., et al. (2020). A Novel Three-lncRNA Signature Predicts the Overall Survival of HNSCC Patients. Ann. Surg. Oncol. 28, 3396–3406. doi:10.1245/s10434-020-09210-1
Johnson, D. E., Burtness, B., Leemans, C. R., Lui, V. W. Y., Bauman, J. E., and Grandis, J. R. (2020). Head and Neck Squamous Cell Carcinoma. Nat. Rev. Dis. Primers 6 (1), 92. doi:10.1038/s41572-020-00224-3
Kitamura, N., Sento, S., Yoshizawa, Y., Sasabe, E., Kudo, Y., and Yamamoto, T. (2020). Current Trends and Future Prospects of Molecular Targeted Therapy in Head and Neck Squamous Cell Carcinoma. Ijms 22 (1), 240. doi:10.3390/ijms22010240
Li, L., Liu, W., Sun, Q., Zhu, H., Hong, M., and Qian, S. (2021a). Decitabine Downregulates TIGAR to Induce Apoptosis and Autophagy in Myeloid Leukemia Cells. Oxidative Med. Cell Longevity 2021, 1–15. doi:10.1155/2021/8877460
Li, Z., Lin, Y., Cheng, B., Zhang, Q., and Cai, Y. (2021c). Identification and Analysis of Potential Key Genes Associated with Hepatocellular Carcinoma Based on Integrated Bioinformatics Methods. Front. Genet. 12, 571231. doi:10.3389/fgene.2021.571231
Longlong, Z., Xinxiang, L., Yaqiong, G., and Wei, W. (2020). Predictive Value of the Texture Analysis of Enhanced Computed Tomographic Images for Preoperative Pancreatic Carcinoma Differentiation. Front. Bioeng. Biotechnol. 8, 719. doi:10.3389/fbioe.2020.00719
Lucá, R., Averna, M., Zalfa, F., Vecchi, M., Bianchi, F., Fata, G. L., et al. (2013). The Fragile X Protein Binds M RNA S Involved in Cancer Progression and Modulates Metastasis Formation. EMBO Mol. Med. 5 (10), 1523–1536. doi:10.1002/emmm.201302847
Meng, L., He, X., Hong, Q., Qiao, B., Zhang, X., Wu, B., et al. (2021a). CCR4, CCR8, and P2RY14 as Prognostic Factors in Head and Neck Squamous Cell Carcinoma Are Involved in the Remodeling of the Tumor Microenvironment. Front. Oncol. 11, 618187. doi:10.3389/fonc.2021.618187
Meng, X., Deng, Y., He, S., Niu, L., and Zhu, H. (2021b). m6A-Mediated Upregulation of LINC00857 Promotes Pancreatic Cancer Tumorigenesis by Regulating the miR-150-5p/E2F3 Axis. Front. Oncol. 11, 629947. doi:10.3389/fonc.2021.629947
Ocadlikova, D., Lecciso, M., Broto, J. M., Scotlandi, K., Cavo, M., Curti, A., et al. (2021). Sunitinib Exerts In Vitro Immunomodulatory Activity on Sarcomas via Dendritic Cells and Synergizes with PD-1 Blockade. Front. Immunol. 12, 577766. doi:10.3389/fimmu.2021.577766
Olivero, C. E., Martínez-Terroba, E., Zimmer, J., Liao, C., Tesfaye, E., Hooshdaran, N., et al. (2020). p53 Activates the Long Noncoding RNA Pvt1b to Inhibit Myc and Suppress Tumorigenesis. Mol. Cel 77 (4), 761–774. doi:10.1016/j.molcel.2019.12.014
Qin, Y., Zheng, X., Gao, W., Wang, B., and Wu, Y. (2021). Tumor Microenvironment and Immune-Related Therapies of Head and Neck Squamous Cell Carcinoma. Mol. Ther. - Oncolytics 20, 342–351. doi:10.1016/j.omto.2021.01.011
Rapa, M., Ciano, S., Ruggieri, R., and Vinci, G. (2021). Bioactive Compounds in Cherry Tomatoes (Solanum Lycopersicum Var. Cerasiforme): Cultivation Techniques Classification by Multivariate Analysis. Food Chem. 355, 129630. doi:10.1016/j.foodchem.2021.129630
Rasheduzzaman, M., Kulasinghe, A., Dolcetti, R., Kenny, L., Johnson, N. W., Kolarich, D., et al. (2020). Protein Glycosylation in Head and Neck Cancers: From Diagnosis to Treatment. Biochim. Biophys. Acta (Bba) - Rev. Cancer 1874 (2), 188422. doi:10.1016/j.bbcan.2020.188422
Reinhold, W. C., Varma, S., Sunshine, M., Elloumi, F., Ofori-Atta, K., Lee, S., et al. (2019). RNA Sequencing of the NCI-60: Integration into CellMiner and CellMiner CDB. Cancer Res. 79 (13), 3514–3524. doi:10.1158/0008-5472.Can-18-2047
Romanowska, K., Sobecka, A., Rawłuszko-Wieczorek, A. A., Suchorska, W. M., and Golusiński, W. (2020). Head and Neck Squamous Cell Carcinoma: Epigenetic Landscape. Diagnostics 11 (1), 34. doi:10.3390/diagnostics11010034
Sawatani, Y., Komiyama, Y., Nakashiro, K.-i., Uchida, D., Fukumoto, C., Shimura, M., et al. (2020). Paclitaxel Potentiates the Anticancer Effect of Cetuximab by Enhancing Antibody-dependent Cellular Cytotoxicity on Oral Squamous Cell Carcinoma Cells In Vitro. Ijms 21 (17), 6292. doi:10.3390/ijms21176292
Shen, J., Feng, X.-p., Hu, R.-b., Wang, H., Wang, Y.-l., Qian, J.-h., et al. (2021). N-methyladenosine Reader YTHDF2-Mediated Long Noncoding RNA FENDRR Degradation Promotes Cell Proliferation in Endometrioid Endometrial Carcinoma. Lab. Invest. 101, 775–784. doi:10.1038/s41374-021-00543-3
Shi, H., Chai, P., Jia, R., and Fan, X. (2020). Novel Insight into the Regulatory Roles of Diverse RNA Modifications: Re-defining the Bridge between Transcription and Translation. Mol. Cancer 19 (1), 78. doi:10.1186/s12943-020-01194-6
Skeate, J. G., Otsmaa, M. E., Prins, R., Fernandez, D. J., Da Silva, D. M., and Kast, W. M. (2020). TNFSF14: LIGHTing the Way for Effective Cancer Immunotherapy. Front. Immunol. 11, 922. doi:10.3389/fimmu.2020.00922
Soldevilla, B., Carretero-Puche, C., Gomez-Lopez, G., Al-Shahrour, F., Riesco, M. C., Gil-Calderon, B., et al. (2019). The Correlation between Immune Subtypes and Consensus Molecular Subtypes in Colorectal Cancer Identifies Novel Tumour Microenvironment Profiles, with Prognostic and Therapeutic Implications. Eur. J. Cancer 123, 118–129. doi:10.1016/j.ejca.2019.09.008
Song, T., Yang, Y., Jiang, S., and Peng, J. (2020). Novel Insights into Adipogenesis from the Perspective of Transcriptional and RNA N6‐Methyladenosine‐Mediated Post‐Transcriptional Regulation. Adv. Sci. 7 (21), 2001563. doi:10.1002/advs.202001563
Su, Y., Zhang, Y., Hua, X., Huang, J., Bi, X., Xia, W., et al. (2021). High-dose Tamoxifen in High-Hormone-Receptor-Expressing Advanced Breast Cancer Patients: a Phase II Pilot Study. Ther. Adv. Med. Oncol. 13, 175883592199343. doi:10.1177/1758835921993436
Sun, Z., Xue, S., Zhang, M., Xu, H., Hu, X., Chen, S., et al. (2020). Aberrant NSUN2-Mediated m5C Modification of H19 lncRNA Is Associated with Poor Differentiation of Hepatocellular Carcinoma. Oncogene 39 (45), 6906–6919. doi:10.1038/s41388-020-01475-w
Tang, J., Liu, J., Li, J., Liang, Z., Zeng, K., Li, H., et al. (2021). Upregulation of SKA3 Enhances Cell Proliferation and Correlates with Poor Prognosis in Hepatocellular Carcinoma. Oncol. Rep. 45 (4), 48. doi:10.3892/or.2021.7999
Tang, Y., He, Y., Shi, L., Yang, L., Wang, J., Lian, Y., et al. (2017). Co-expression of AFAP1-AS1 and PD-1 Predicts Poor Prognosis in Nasopharyngeal Carcinoma. Oncotarget 8 (24), 39001–39011. doi:10.18632/oncotarget.16545
Teachey, D. T., and O’Connor, D. (2020). How I Treat Newly Diagnosed T-Cell Acute Lymphoblastic Leukemia and T-Cell Lymphoblastic Lymphoma in Children. Blood 135 (3), 159–166. doi:10.1182/blood.2019001557
Vaid, A., Jaladanki, S. K., Xu, J., Teng, S., Kumar, A., Lee, S., et al. (2021). Federated Learning of Electronic Health Records to Improve Mortality Prediction in Hospitalized Patients with COVID-19: Machine Learning Approach. JMIR Med. Inform. 9 (1), e24207. doi:10.2196/24207
Wang, M., Liu, J., Zhao, Y., He, R., Xu, X., Guo, X., et al. (2020a). Upregulation of METTL14 Mediates the Elevation of PERP mRNA N6 Adenosine Methylation Promoting the Growth and Metastasis of Pancreatic Cancer. Mol. Cancer 19 (1), 130. doi:10.1186/s12943-020-01249-8
Wang, P., Jin, M., Sun, C.-h., Yang, L., Li, Y.-s., Wang, X., et al. (2018). A Three-lncRNA Expression Signature Predicts Survival in Head and Neck Squamous Cell Carcinoma (HNSCC). Biosci. Rep. 38 (6), bsr20181528. doi:10.1042/bsr20181528
Wang, Q., and Liu, X. (2020). VDAC Upregulation and αTAT1-mediated α-tubulin A-cetylation C-ontribute to T-anespimycin-induced A-poptosis in Calu-1 C-ells. Oncol. Rep. 44 (6), 2725–2734. doi:10.3892/or.2020.7789
Wang, S., Su, W., Zhong, C., Yang, T., Chen, W., Chen, G., et al. (2020b). An Eight-CircRNA Assessment Model for Predicting Biochemical Recurrence in Prostate Cancer. Front. Cel Dev. Biol. 8, 599494. doi:10.3389/fcell.2020.599494
Wang, W., Zhu, D., Zhao, Z., Sun, M., Wang, F., Li, W., et al. (2021b). RNA Sequencing Reveals the Expression Profiles of circRNA and Identifies a Four-circRNA Signature Acts as a Prognostic Marker in Esophageal Squamous Cell Carcinoma. Cancer Cel Int 21 (1), 151. doi:10.1186/s12935-021-01852-9
Wang, Y., Hui, J., Li, R., Fu, Q., Yang, P., Xiao, Y., et al. (2020c). GBX2, as a Tumor Promoter in Lung Adenocarcinoma, Enhances Cells Viability, Invasion and Migration by Regulating the AKT/ERK Signaling Pathway. J. Gene Med. 22 (2), e3147. doi:10.1002/jgm.3147
Wang, Y., Lu, J.-H., Wu, Q.-N., Jin, Y., Wang, D.-S., Chen, Y.-X., et al. (2019). LncRNA LINRIS Stabilizes IGF2BP2 and Promotes the Aerobic Glycolysis in Colorectal Cancer. Mol. Cancer 18 (1), 174. doi:10.1186/s12943-019-1105-0
Wang, Y., Wang, S., Ren, Y., and Zhou, X. (2020d). The Role of lncRNA Crosstalk in Leading Cancer Metastasis of Head and Neck Squamous Cell Carcinoma. Front. Oncol. 10, 561833. doi:10.3389/fonc.2020.561833
Xu, X., Zhou, E., Zheng, J., Zhang, C., Zou, Y., Lin, J., et al. (2021). Prognostic and Predictive Value of m6A "Eraser" Related Gene Signature in Gastric Cancer. Front. Oncol. 11, 631803. doi:10.3389/fonc.2021.631803
Xu, Y., Liu, J., Chen, W.-J., Ye, Q.-Q., Chen, W.-T., Li, C.-L., et al. (2020b). Regulation of N6-Methyladenosine in the Differentiation of Cancer Stem Cells and Their Fate. Front. Cel Dev. Biol. 8, 561703. doi:10.3389/fcell.2020.561703
Yang, D., Ma, Y., Zhao, P., Ma, J., and He, C. (2019). Systematic Screening of Protein-Coding Gene Expression Identified HMMR as a Potential Independent Indicator of Unfavorable Survival in Patients with Papillary Muscle-Invasive Bladder Cancer. Biomed. Pharmacother. 120, 109433. doi:10.1016/j.biopha.2019.109433
Yang, R., Elsaadi, S., Misund, K., Abdollahi, P., Vandsemb, E. N., Moen, S. H., et al. (2020). Conversion of ATP to Adenosine by CD39 and CD73 in Multiple Myeloma Can Be Successfully Targeted Together with Adenosine Receptor A2A Blockade. J. Immunother. Cancer 8 (1), e000610. doi:10.1136/jitc-2020-000610
Yang, Y., Hsu, P. J., Chen, Y.-S., and Yang, Y.-G. (2018). Dynamic Transcriptomic m6A Decoration: Writers, Erasers, Readers and Functions in RNA Metabolism. Cell Res 28 (6), 616–624. doi:10.1038/s41422-018-0040-8
Ye, M., Dong, S., Hou, H., Zhang, T., and Shen, M. (2021). Oncogenic Role of Long Noncoding RNAMALAT1 in Thyroid Cancer Progression through Regulation of the miR-204/IGF2BP2/m6A-MYC Signaling. Mol. Ther. - Nucleic Acids 23, 1–12. doi:10.1016/j.omtn.2020.09.023
Yin, X., Han, S., Song, C., Zou, H., Wei, Z., Xu, W., et al. (2019). Metformin Enhances Gefitinib Efficacy by Interfering with Interactions between Tumor-Associated Macrophages and Head and Neck Squamous Cell Carcinoma Cells. Cell Oncol. 42 (4), 459–475. doi:10.1007/s13402-019-00446-y
Yu, K., Xie, W., Wang, L., and Li, W. (2021). ILRC: a Hybrid Biomarker Discovery Algorithm Based on Improved L1 Regularization and Clustering in Microarray Data. BMC Bioinformatics 22 (1), 514. doi:10.1186/s12859-021-04443-7
Zalfa, F., Panasiti, V., Carotti, S., Zingariello, M., Perrone, G., Sancillo, L., et al. (2017). The Fragile X Mental Retardation Protein Regulates Tumor Invasiveness-Related Pathways in Melanoma Cells. Cel Death Dis 8 (11), e3169. doi:10.1038/cddis.2017.521
Zhang, C., Zhang, Z., Zhang, G., Zhang, Z., Luo, Y., Wang, F., et al. (2020a). Clinical Significance and Inflammatory Landscapes of a Novel Recurrence-Associated Immune Signature in Early-Stage Lung Adenocarcinoma. Cancer Lett. 479, 31–41. doi:10.1016/j.canlet.2020.03.016
Zhang, L.-P., Ren, H., Du, Y.-X., and Wang, C.-F. (2020b). Prognostic Value of the Preoperative Fibrinogen-To-Albumin Ratio in Pancreatic Ductal Adenocarcinoma Patients Undergoing R0 Resection. Wjg 26 (46), 7382–7404. doi:10.3748/wjg.v26.i46.7382
Zhang, M., Zhai, Y., Zhang, S., Dai, X., and Li, Z. (2020c). Roles of N6-Methyladenosine (m6A) in Stem Cell Fate Decisions and Early Embryonic Development in Mammals. Front. Cel Dev. Biol. 8, 782. doi:10.3389/fcell.2020.00782
Zhao, Y., Xu, L., Wang, X., Niu, S., Chen, H., and Li, C. (2021). A Novel Prognostic mRNA/miRNA Signature for Esophageal Cancer and its Immune Landscape in Cancer Progression. Mol. Oncol. 15, 1088–1109. doi:10.1002/1878-0261.12902
Zheng, J., Zhang, T., Guo, W., Zhou, C., Cui, X., Gao, L., et al. (2020). Integrative Analysis of Multi-Omics Identified the Prognostic Biomarkers in Acute Myelogenous Leukemia. Front. Oncol. 10, 591937. doi:10.3389/fonc.2020.591937
Zhou, X., Han, J., Zhen, X., Liu, Y., Cui, Z., Yue, Z., et al. (2020). Analysis of Genetic Alteration Signatures and Prognostic Values of m6A Regulatory Genes in Head and Neck Squamous Cell Carcinoma. Front. Oncol. 10, 718. doi:10.3389/fonc.2020.00718
Zhu, P., He, F., Hou, Y., Tu, G., Li, Q., Jin, T., et al. (2021). A Novel Hypoxic Long Noncoding RNA KB-1980E6.3 Maintains Breast Cancer Stem Cell Stemness via Interacting with IGF2BP1 to Facilitate C-Myc mRNA Stability. Oncogene 40 (9), 1609–1627. doi:10.1038/s41388-020-01638-9
Zhu, Y., Huang, Y., Chen, L., Guo, L., Wang, L., Li, M., et al. (2021). Up-Regulation of SLC26A6 in Hepatocellular Carcinoma and its Diagnostic and Prognostic Significance. Crit. Rev. Eukaryot. Gene Expr. 31 (5), 79–94. doi:10.1615/CritRevEukaryotGeneExpr.2021039703
Keywords: RNA methylation, long non-coding RNA, prognostic signature, head and neck squamous cell carcinoma, epigenetic change
Citation: Wang E, Li Y, Ming R, Wei J, Du P, Zhou P, Zong S and Xiao H (2021) The Prognostic Value and Immune Landscapes of a m6A/m5C/m1A-Related LncRNAs Signature in Head and Neck Squamous Cell Carcinoma. Front. Cell Dev. Biol. 9:718974. doi: 10.3389/fcell.2021.718974
Received: 01 June 2021; Accepted: 05 November 2021;
Published: 30 November 2021.
Edited by:Nejat Dalay, Istanbul University, Turkey
Reviewed by:Leong-Perng Chan, Kaohsiung Medical University Hospital, Taiwan
Narendra Singh, Stowers Institute for Medical Research, United States
Copyright © 2021 Wang, Li, Ming, Wei, Du, Zhou, Zong and Xiao. 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.
‡These authors have contributed equally to this work and share first authorship
†Present Address: Yang Li, Department of Prosthodontics, Peking University School and Hospital of Stomatology, National Engineering Laboratory for Digital and Material Technology of Stomatology, National Clinical Research Center for Oral Diseases, Beijing Key Laboratory of Digital Stomatology, Beijing, China