Identification of Expression Patterns and Potential Prognostic Significance of m5C-Related Regulators in Head and Neck Squamous Cell Carcinoma

5-Methylcytosine (m⁵C) methylation is a major epigenetic technique of RNA modification and is dynamically mediated by m⁵C “writers,” “erasers,” and “readers.” m⁵C RNA modification and its regulators are implicated in the onset and development of many tumors, but their roles in head and neck squamous cell carcinoma (HNSCC) have not yet been completely elucidated. In this study, we examined expression patterns of core m⁵C regulators in the publicly available HNSCC cohort via bioinformatic methods. The differentially expressed m⁵C regulators could divide the HNSCC cohort into four subgroups with distinct prognostic characteristics. Furthermore, a three-gene expression signature model, comprised of NSUN5, DNMT1, and DNMT3A, was established to identify individuals with a high or low risk of HNSCC. To explore the underlying mechanism in the prognosis of HNSCC, screening of differentially expressed genes, followed by the analysis of functional and pathway enrichment, from individuals with high- or low-risk HNSCC was performed. The results revealed a critical role for m⁵C RNA modification in two aspects of HNSCC: (1) dynamic m⁵C modification contributes to the regulation of HNSCC progression and (2) expression patterns of NSUN5, DNMT1, and DNMT3A help to predict the prognosis of HNSCC.

Introduction

According to the most recent report, head and neck squamous cell carcinoma (HNSCC) is a relatively lethal type of cancer, and HNSCC ranks among the top six in terms of incidence and mortality, seriously threatening public health and the quality of life of patients with HNSCC worldwide (1). Despite consecutive forms of treatment, namely, surgery, chemotherapy, and radiotherapy, and considerable advancement in the therapeutic schedule for HNSCC, the 5-year survival rate of patients with HNSCC remains far from satisfactory, owing to end-stage diagnosis, rapid development, high recurrence rate, and induction of metastasis to distant sites (2). Mounting evidence has demonstrated that molecular markers hold a promising function not only because of their prognostic value but also because of their role as molecular targets. For example, HNSCC therapy has undergone a significant change in recent years with the development of precision medicines, such as bevacizumab, against vascular endothelial growth factor (VEGF)/VEGF receptor (VEGFR) based on the special gene expression signature, and prognostic outcomes in HNSCC (3–5). More recently, increasing knowledge has brought into focus the features of several macromolecules (protein, RNA, DNA, and sugar) that are involved in tumorigenesis and progression, especially in epigenetic modifications; and the value of these features as prognostic indicators and potential therapeutic targets has been gradually recognized (6–8). To date, the molecular mechanism of HNSCC occurrence and development has not yet been completely elucidated. Therefore, it is imperative to gain a deeper insight into the molecular mechanism of carcinogenesis of HNSCC and thereby provide valuable detection and effective targets for patients with HNSCC.

To date, RNA modification methods have become prominent because of the development of detection technologies and the realization that RNA not only serves as the intermediate molecule for translation or as an auxiliary function for protein synthesis (rRNA and tRNA) but also acts as a functional regulator for the transmission of genetic signals (lncRNA and miRNA) (9, 10). Among these methods, RNA methylation is one of the most common techniques in the epigenetic modification of posttranscriptional RNA, even though extensive effort has been made on studying protein and DNA modifications (11). Generally, RNA methylation predominantly includes m⁶A, m⁵C, m¹A, and m⁷G, among which m⁶A and m⁵C modification techniques are two of the most major and most representative types of posttranscriptional RNA modification in over 170 chemical modification schemes (12–15). The m⁶A modification technique has been predominantly studied in HNSCC, but that of m⁵C does not attract much attention since the current knowledge of its functions is limited to regulation of the exportation of mRNA and the maintenance of the structure and stability of mRNA (12, 16–18). Its modulating effects have been characterized to be reversible and dynamic, similar to those of histone and DNA modification methods, owing to the involvement of m⁵C writers, erasers, and readers (19, 20). Taking these regulators further, m⁵C was established by adding the methyl group through a number of methyl transferases (NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, TRDMT1, DNMT3A, DNMT3B, and DNMT1), was removed by demethylase (TET2 and TET3), and was recognized by binding proteins (ALYREF and YBX1), which were jargonized as “writers,” “erasers,” and “readers,” respectively (12, 13, 21–23).

At present, accumulating evidence has demonstrated that the aberrant expression of m⁵C RNA regulators and specific methylated genes is involved in the pathogenesis of abnormal differentiation in progenitors, fertility damage in males, and cancer oncogenesis (24–27). For instance, DNMT1- and EZH2-mediated epigenetic silencing promotes the progression of glioblastoma and gastric cancer (28). Furthermore, NSUN6 promotes the activation of breast cancer metastasis by incorporating the adaptor proteins, LLGL2, and lncRNA MAYA, to accumulate YAP1 in the nucleus for transcriptional activation (29). In addition, one of the main m⁵C RNA writers, NSUN2, was also noted to be overexpressed in different types of cancer and was deemed an effective prognostic biomarker (25, 27). However, the gene features and prognostic values of m⁵C-related regulators in HNSCC remain obscure and need in-depth investigation.

In this study, we systematically analyzed and evaluated the expression patterns of 15 widely studied m⁵C-related regulators in 501 tumor and 44 normal control tissues and the association between clinicopathological and survival parameters from The Cancer Genome Atlas (TCGA) database.

Materials and Methods

HNSCC Dataset Acquisition and Bioinformatic Analysis

The transcriptome TCGA-HNSCC datasets and the corresponding clinical features applied in this study were obtained from the TCGA database via the GDC Data Portal, as described earlier (https://portal.gdc.cancer.gov/) (30). In total, 501 tumor and 44 normal control datasets from 527 patients were available for further experimental procedure. In addition, the mutation data and the expression values in the pan-cancer analysis of three selected risk genes were obtained from the cBioPortal database (https://www.cbioportal.org/) and the Gene Expression Profiling Interactive Analysis 2 (GEPIA2) database (http://gepia2.cancer-pku.cn/) (31–33). The differentially expressed genes (DEGs) were generated by the R (3.6.0) package “limma,” the heat map by “pheatmap,” the volcano plots and the bubble plots by “ggplot2,” and the chord plots by “GOplot.”

Landscape of m⁵C RNA Methylation Regulators

In total, 15 m⁵C-associated regulators composed of 11 writers (NOP2, NSUN2, NUSN3, NSUN4, NSUN5, NSUN6, NSUN7, DNMT1, DNMT3A, DNMT3B, and TRDMT1), two erasers (TET2 and TET3), and two readers (ALYREF and YBX1) were retrieved from the published literature. The expression profile of these 15 regulators, accompanied by clinicopathological parameters, was then systematically extracted and analyzed in patients with HNSCC.

Consensus Clustering of m⁵C-Related Regulators

To better investigate and construct the distinct classification model, the 13 selected m⁵C-related regulators were screened out to primarily conduct the consensus clustering analysis using the “ConsensusClusterPlus” package of R (3.6.0) (34). Furthermore, the survival analysis of different clusters was performed to determine the best clustering in HNSCC samples.

Construction of Prognostic Prediction Model

To obtain a better prediction of m⁵C RNA methylation regulators in HNSCC, first, we calculated the hazard ratio (HR), which indicates the result of comparing the hazard function between individuals who are exposed to the hazard function and those who are not, and 95% confidence interval (CI) of the m⁵C-related regulators to identify the appropriate candidate genes by univariate Cox regression analysis. Second, the appropriate candidate m⁵C-related regulators were assigned and built for potential HNSCC prognostic signatures using LASSO Cox regression, which was calculated with the formula described below:

$\begin{array}{l}RiskScore={\displaystyle \sum _{i=1}^n}coefi\times x_i\end{array}$

where n represents the number of module RNAs, coefi is the coefficient, and x_i denotes the z-score-transformed relative expression level (log₂ (FPKM + 1)) for each gene. Then, the HNSCC cohort from TCGA was divided into two subgroups, high risk and low risk, on account of the median risk score. The Kaplan–Meier survival analysis/risk prediction model was applied to estimate the prognostic values of the candidate risk genes.

Protein-Protein Interactions and Functional Annotations Analysis

Protein-protein interactions (PPIs) of the 15 m⁵C-related regulators with a combined confidence score > 0.4 and the DEGs between the two risk subgroups were evaluated via the STRING database (https://string-db.org/) and visualized using Cytoscape (3.7.1) (35, 36). The Reactome (https://reactome.org/), the Kyoto Encyclopedia of Genes and Genomes (KEGG), and the Annotation, Visualization, and Integrated Discovery (DAVID) databases (https://david.ncifcrf.gov/) were used to evaluate the enriched functional annotations (37–41).

Immunofluorescence Analysis

Following the guidelines set by the Research Ethics Committee of the Shanghai Ninth People's Hospital, which is affiliated to the Shanghai Jiao Tong University School of Medicine, 10 pairs of tumors, in which the pathological results revealed HNSCC, were selected for this study. Normal human oral mucosal epithelial tissues obtained from the oral mucosa during the surgical resection of HNSCC were used as a control. Paraffin-embedded samples corresponding to the most representative tumor area on H&E-stained slides were selected for performing the IF assay. Briefly, the original fresh-frozen IF sections with a thickness of 6 μm were acquired through cryosection, air-dried for 10 min, fixed with acetone, washed with phosphate buffered saline (PBS) (1 ×), and incubated at room temperature for 2 h with NSUN5 (Proteintech Group, Rosemont, IL, USA, 1:50), DNMT1 (Bioworld Technology, Inc., St. Louis Park, MN, USA, 1:50), and DNMT3A (Novus Biologicals, Littleton, CO, USA, 1:170) primary antibodies. The sections were rinsed with sterile PBS (1 ×) and then incubated with a Goat Anti-Rabbit IgG secondary antibody (Abcam, Cambridge, UK, 1:200) protected from light. Next, the nuclei of the sections were counterstained with DAPI (Beyotime, Shanghai, China). Primary antibodies were replaced with PBS as a negative control.

Statistical Analysis

Statistical analysis was conducted using R software (3.6.0). Moreover, the DEGs between the two groups were analyzed by the Student's t-test. p < 0.05 was deemed statistically significant.

Results

Expression Patterns of m⁵C RNA Methylation Regulators in HNSCC

The overall flowchart of the procedure that was applied to in this study and that elucidated the risk score for investigating the prognostic values of HNSCC is summarized in Figure 1. To decipher the expression of essential biological functions of m⁵C-related regulators in HNSCC, we first downloaded and extracted the efficacious gene expression data from the TCGA database. Expression levels of the individual m⁵C RNA methylation catalase in HNSCC and control samples are presented in the heat map (Figure 2A). Among the 545 cases, 14 of the 15 m⁵C methylation regulatory genes were differentially expressed between tumors and healthy samples, with p < 0.05. Specifically, NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, DNMT1, DNMT3A, DNMT3B, TET2, TET3, ALYREF, and YBX1 exhibited different expression patterns. However, no distinct discrepancies in the TRDMT1 (p = 0.83) expression were analyzed in HNSCC tissues when compared with normal tissues (Figure 2B). Of these 14 genes, the majority of writers (NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, DNMT3A, DNMT3B, and DNMT1) were more substantially upregulated in the HNSCC samples compared with normal tissues, with the exception of NSUN7, which showed an opposite expression trend with the other 9 m⁵C writers (p < 0.05). In addition, in HNSCC tumor tissues, the expression level of the m⁵C eraser TET3 was significantly elevated, while that of TET2 was downregulated. Similarly, the readers ALYREF and YBX1 showed a higher expression level in HNSCC compared to normal tissues. In summary, we concluded that m⁵C RNA methylation regulators had distinct expression changes in HNSCC and corresponding normal tissues.

FIGURE 1

Figure 1. Workflow of the study design and different analyses conducted in the study.

FIGURE 2

Figure 2. Expressional and interactive landscape of the m⁵C modification signature in head and neck squamous cell carcinoma (HNSCC). (A) The heat map shows the expression levels of 15 m⁵C RNA decoration regulators in each clinical sample (N, normal; T, tumor; *p < 0.05; **p < 0.01; ***p < 0.001). (B) A violin plot was applied to demonstrate the significantly differentially expressed m⁵C RNA modification regulators between tumor tissues and those of normal control. (C) The protein–protein interaction (PPI) network of 15 m⁵C-related regulators was constructed to visualize the interaction (the size of the node is applied to reflect the degree of regulator; the size of the line is applied to denote the combined_score). (D) Pearson correlation analysis was delineated to determine the correlation among 15 selected m⁵C-related regulators in The Cancer Genome Atlas (TCGA) HNSCC cohort.

Interaction and Correlation Patterns Among the m⁵C RNA Methylation Regulators in HNSCC

To investigate the associations between the main m⁵C RNA methylation regulators, we built a PPI network. The interactions among the 15 m⁵C-related regulators are shown in Figure 2C. TRDMT1 appeared to be the hub gene of the interaction network and was predominantly associated with most of the m⁵C RNA methylation regulators, except for ALYREF and YBX1. Since the PPI network did not provide details of correlation, we performed a further correlation analysis on HNSCC, as shown in Figure 2D. There was a close correlation between TET2 and TET3, two members of the TET gene family. Except for NSUN5 and DNMT3B, NSUN6 was correlated with the other 12 m⁵C RNA methylation regulators. Moreover, the expression of DNMT1 was positively related to the other m⁵C-related methylation genes but not to TRDMT1 and NSUN5. Similarly, NSUN4 also shows a positive relationship with other m⁵C-related methylation regulators except for NSUN5 and DNMT3B. Furthermore, it is worthwhile to note that both erasers, TET2 and TET3, were most negatively correlated with NSUN5 among all the correlations of m⁵C RNA methylation regulators with the Pearson's correlation coefficient of −0.34 and −0.19, respectively.

Identification of Four Clusters of HNSCC Samples With Different Clinical Outcomes and Characteristics

Transcriptome data of 545 HNSCC samples from the TCGA database were used for consensus clustering analysis. From the differentially expressed m⁵C-related methylation regulators described above, we used 13 m⁵C-related genes (NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, DNMT1, DNMT3A, DNMT3B, TET2, TET3, ALYREF, and YBX1) for further research. Based on the expression similarity profiling of the 13 m⁵C-related regulators combined with consensus clustering cumulative distribution function (CDF) and relative change in the area under the CDF curve, as shown in Figures 3A,B, k = 4 was deemed as the appropriate number of clusters when the clustering stability datasets varied from k = 2 to k = 10 (Figures 3C,D). Then, the HNSCC samples with survival parameters were classified accordingly into four groups. A noticeably shorter overall survival (OS) was observed in HNSCC cases in cluster 1 compared with the other clusters. We highly postulate that the expression of these 13 m⁵C-related genes can divide the HNSCC samples into four groups with distinct prognoses (Figure 3E). Additionally, to better predict the clinicopathological characteristics of HNSCC with these four subgroups, a heat map was applied to present the significant difference in grade and survival condition with both p < 0.01, while no huge difference was witnessed with other features, such as TNM classification, gender, age, and clinical stage (Figure 3F). Overall, we can conclude that the expression features of the 13m⁵C-related genes were associated with the grade and survival condition of patients with HNSCC.

FIGURE 3

Figure 3. Differential patterns of expressional and clinical features of patients with TCGA HNSCC in cluster 1–4 subgroups. (A) Consensus clustering cumulative distribution function (CDF) and (B) relative change in the area under the CDF curve from k = 2 to k = 9. (C) Consensus clustering matrix was assessed by k = 4. (D) The tracking plot of the HNSCC samples. (E) The survival analysis for the four clusters by the Kaplan–Meier method. (F) The clinicopathological characteristics of the four clusters were determined by coexpression patterns of the consensus expression of 13 selected m⁵C-related regulators (*p < 0.05; **p < 0.01; ***p < 0.001).

Construction of a Three-Gene Risk Signature With Distinct Prognostic Value

To investigate the prognostic role of m⁵C RNA methylation regulators in HNSCC, univariate Cox regression analysis was applied to the 13 m⁵C-related gene expression profiles. According to the details contained in these results (Figure 4A), three (NSUN5, DNMT1, and DNMT3A) of the 13 genes that presented a significant prognostic value (p < 0.1) were specifically selected to establish the risk signature. Among these three selected genes, DNMT1 and DNMT3A were protective genes with HR <1, while NSUN5 served as a risk factor with HR > 1.

FIGURE 4

Figure 4. The prognostic risk signature was created with the three m⁵C RNA methylation regulators. (A) The process of analyzing the signature regulators by the univariate Cox regression model. (B,C) LASSO Cox regression analysis shows the coefficients of the three selected regulators (NSUN5, DNMT1, and DNMT3A). (D) The Kaplan–Meier overall survival analysis of the two subgroups of TCGA HNSCC cohorts assigned to high- and low-risk groups followed by risk scores. (E) Significant differences were observed in the area of M stage and grade across the high- and low-risk subgroups. (F) Univariate and (G) multivariate Cox regression showed the relationship between clinicopathological features and the risk score.

Then the three screened genes with prognostic values were applied to build the survival risk model using LASSO Cox regression. The coefficients of individual candidate genes were generated based on the minimum criteria (Figure 4B). Subsequently, the risk score of each patient with HNSCC from the TCGA database was calculated as follows: risk score = (−0.029572) *expression of DNMT1 + (−0.038603) *expression of DNMT3A + (0.021496) *expression of NSUN5 (Figure 4C). Afterward, patients with HNSCC were divided into two subgroups, low risk and high risk, according to the median bound. It was observed that patients in the high-risk subgroup had a significantly shorter OS than patients in the low-risk subgroup (Figure 4D, p < 0.05). Taken together, these results suggest that the three risk genes can be used as a predictor for the clinical outcomes of HNSCC.

The Prognostic Patterns of Signature-Based Risk Scores Were Associated With Clinical Features in HNSCC

To evaluate the clinical parameters with the three selected m⁵C methylation regulators in the two subgroups, a heat map was generated to visualize the relationship pattern. The expression of DNMT1 and DNMT3A was high in the low-risk subgroup, whereas, the expression of NSUN5 mainly emerged in the high-risk subgroup. Furthermore, significant differences were disclosed in this heat map concerning the M stage and grade (Figure 4E, p < 0.01).

Next, we further performed univariate and multivariate Cox regression analyses to confirm whether the risk signature was an independent prognostic element. Remarkably, age (p = 0.002, HR = 1.021, 95% CI = 1.008–1.034), gender (p = 0.036, HR = 0.726, 95% CI = 0.538–0.979), T stage (p = 0.013, HR = 1.164, 95% CI = 1.033–1.312), N stage (p < 0.001, HR = 1.278, 95% CI = 1.129–1.447), and the risk score (p < 0.001, HR = 2.545, 95% CI = 1.548–4.184) were determined to be independent prognostic factors from univariate analysis (Figure 4F), while these parameters were calculated by the multivariate Cox regression model. The valuable factors was reduced to age (p = 0.008, HR = 1.018, 95% CI = 1.005–1.032), N stage (p = 0.012, HR = 1.203, 95% CI = 1.042–1.390), and the risk score (p = 0.004, HR = 2.048, 95% CI = 1.256–3.340) (Figure 4G).

External Confirmation of the Expression Patterns and Genetic Alterations of Three Risk Signature Genes in HNSCC

To better explore the expression profiles of NSUN5, DNMT1, and DNMT3A in human organs and corresponding tumors, the GEPIA2 database was utilized to investigate the distinct organ expression features of these three prognostic markers. As shown in Supplementary Figure 1 and Supplementary Table 1, these three markers were almost overexpressed in tumors compared with normal organs in 26 pairs of comparison. Next, further comparison of the three genes at the transcriptional level among the 31 human pan-cancer samples was conducted, and an almost identical result was obtained from these contrasts. This validation indicates that NSUN5, DNMT1, and DNMT3A were nearly overexpressed in human tumors compared with normal samples (Figures 5A–C).

FIGURE 5

Figure 5. The expression and genetic alteration patterns of three risk genes of various types of cancer in human. The expression values of (A) NSUN5, (B) DNMT1, and (C) DNMT3A in pan-cancers and corresponding normal tissues. (D) Expression levels of three risk genes in HNSCC and normal controls based on the GEPIA2 database. (E) The profiling of the expression alterations of the three genes in the HNSCC sample (n = 1,332) from TCGA datasets.

Considering the complexity of the posttranscriptional regulation and the protein expression, we further employed the GEPIA2 database and IF at the protein level to evaluate and verify the expression of three m⁵C signature writers. The GEPIA2 database contains TCGA and Genotype-Tissue Expression (GTEx) data with a larger sample population, which may increase the statistical confidence in a precise estimate. The mRNA expression profiles illustrate that all the three risk signature genes were upregulated in GEPIA2, and similar results were obtained from IF (Figures 5D, 6). Among the 1,330 patients with HNSCC, 92 patients (7.1%) had genetic mutations of NSUN5, DNMT1, or DNMT3A (Figure 5E).

FIGURE 6

Figure 6. Immunofluorescence analysis of NSUN5, DNMT1, and DNMT3A in tissues of HNSCC tumors and normal oral epithelium (×100). Compared with oral epithelium controls (A'–C'), a higher expression of (A) NSUN5, (B) DNMT1, and (C) DNMT3A was detected in HNSCC tumor tissues.

Functional Annotation and Pathways of Two Risk Subgroups Determined by Three Prognostic Genes

The results mentioned above suggest that the two risk subgroups may be closely related to the prognostic capacity of patients with HNSCC. Next, we sought to explore the potential markers targeting the two subgroups and identified the involved biological functions. In total, 725 DEGs, with 359 being upregulated (p < 0.05, log₂FC > 0.5) and 366 being downregulated (p < 0.05, log₂FC < −0.5), were screened in the high-risk subgroup by the limma package (p < 0.05). Among them, we ranked the DEGs followed by log₂FC. From a total of 20 polar DEGs, 10 upregulated mRNAs (KRT83, BPIFA1, KRT79, FABP4, LOR, LINC01214, FLG2, IGFL2, DSC1, and S100A12) and 10 downregulated mRNAs (PRR27, NPS, STATH, DEFA6, FGF3, TKTL1, BPIFA2, IGKV2D-28, PIP, and ZG16B), were found to be at the top of rankings with the highest |log₂FC| (Figure 7A). To better understand the interactions between the upregulated and downregulated genes, we also evaluated the PPIs via the STRING database and visualized them by Cytoscape (Supplementary Figures 2, 3). In addition, to better summarize the function of DEGs, Reactome, GO, and KEGG analyses were performed to illustrate the functional annotations of DEGs using the GOplot and ggplot2 packages.

FIGURE 7

Figure 7. The profiling of differentially expressed genes (DEGs) and their functional annotations from the two risk subgroups. (A) The volcano plot shows the DEGs in the two subgroups. The top 10 upregulated and downregulated genes with the highest |log₂FC| value are labeled with blue dots. Significantly enriched the Reactome pathway analysis of (B) 359 upregulated genes and (C) 366 downregulated genes was visualized using the package GOplot in R. The plot includes eight Reactome pathway terms. The genes with at least two terms and the terms that included at least eight genes are demonstrated in this plot. log₂FC represents the difference between the high- and low-risk subgroups. The bubble chart shows the top 10 significant enrichment terms of (D) 359 upregulated genes and (E) 366 downregulated genes in the Reactome database.

The GOplot and ggplot2 data of the relationship between the listed 359 upregulated and 366 downregulated genes and their corresponding metabolic pathways in the Reactome database, together with the log₂FC of the two-part genes, are presented in Figures 7B–E (p < 0.05). Furthermore, we also applied ggplot2 to demonstrate the functional annotations in the GO database. The leading highly enriched GO terms of biological process (BP), cellular component (CC), and molecular function (MF) in the high-risk subgroup were “translation,” “cytoplasm,” and “protein binding,” respectively (Figures 8A–E, p < 0.05). On the other hand, in the low-risk subgroup, the principal terms of these three aspects degenerated into “positive regulation of transcription from RNA polymerase II promoter,” “plasma membrane,” and “ATP binding” (Figures 8B–F, p < 0.05). Moreover, the enriched signaling cascade for the 359 upregulated and 366 downregulated genes identified by the analysis of the KEGG pathway was selected according to log₂FC (p < 0.05). In the high-risk subgroup, the top-ranking terms were associated with “ribosome,” “oxidative phosphorylation,” “non-alcoholic fatty liver disease,” and “synaptic vesicle cycle,” and so on (p < 0.05); among them, “oxidative phosphorylation” was testified to be related to larger tumor size of HNSCC (Figure 8G) (42). Verification of the low-risk subgroup with KEGG analysis also revealed that “cell cycle” and “HTLV-infection” were linked with these genes (Figure 8H).

FIGURE 8

Figure 8. Bubble chart of the enriched terms in GO and KEGG. Functional annotations of upregulated genes in the high-risk subgroup and low-risk subgroup are analyzed by (A,B) GO biological processes, (C,D) cellular component, (E,F) molecular function, and (G,H) KEGG pathways.

Discussion

Recently, increasing evidence has demonstrated that HNSCC is a complex and heterogeneous disease that is attributed to the combination of virus infection, environmental risk factors, and genetic predisposition. Of note, tobacco smoking and alcohol abuse are considered to be the leading carcinogenic factors for HNSCC (43). In addition, m⁵C RNA modification in HNSCC has garnered substantial attention among researchers worldwide. Its functions should include numerous BPs, such as mRNA export, RNA stability, translation, and alternative splicing (12, 17, 18, 44–47). m⁵C RNA modification can be detected in most types of RNA and is associated with a wide range of disorders (13, 48). In particular, abnormal m⁵C methylation has been implicated in the development of many malignant tumors, namely human skin squamous cell carcinomas and breast cancer (8). Although, Xue et al. (49) reported that the gene significance of m⁵C regulators can predict the prognosis of patients with HNSCC, the role of m⁵C modification in HNSCC is still obscure, and in-depth investigations in the field are urgent. In this study, we analyzed the expression patterns of 15 m⁵C regulators in HNSCC and constructed a three-gene risk signature to predict the prognosis of patients with HNSCC.

Because of the advancements in sensitive, quantitative, and specific technologies, the identification of modification techniques on the low abundance RNA m⁵C regulators has been brought into focus. Mounting evidence has shown that these regulators, namely m⁵C “writers,” “erasers,” and “readers,” modulate the occurrence and progression of tumors primarily through their methylation function. For example, the well-identified “writer” NSUN5 exhibits tumor-suppressing characteristics in gliomas. DNA methylation-associated epigenetic silencing of NSUN5 is observed in human gliomas, and it helps glioma cells overcome hostile stress conditions (50). Chen et al. (21) revealed that YBX1 is an m⁵C “reader” that recognizes m⁵C-modified mRNAs in human urothelial carcinoma of the bladder (UCB) and maintains the stability of its target mRNAs. Moreover, YBX1 targets the m⁵C methylation site in the HDGF 3′ untranslated region to drive UCB pathogenesis (21). Notably, most m⁵C regulators have been reported to participate in cancer pathogenesis via non-methylated pathways. Xu et al. revealed the critical role of YBX1 in modulating abnormal ubiquitination in hepatocellular carcinoma (HCC). Protection of YBX1 from PRP19-mediated ubiquitination degradation by circRNA-SORE, a newly discovered circRNA highly expressed in HCC, increases sorafenib resistance in patients with HCC (51).

In this study, we intend to investigate the expression profile of 15 m⁵C-related regulators in HNSCC. Analysis of the TCGA HNSCC cohort revealed that 14, namely, NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, DNMT1, DNMT3A, DNMT3B, TET2, TET3, ALYREF, and YBX1, out of 15 m⁵C-related RNA regulators exhibited different expression hallmarks among tumors and normal controls. Considering that m⁵C regulators are differentially expressed in other tumors and involved in the regulation of their pathogenesis, interpretation of the analyzed results somewhat indicates that the differentially expressed regulators may affect HNSCC development and therapy. For example, the levels of DNMT3B, NSUN2, DNMT3A, NOP2, DNMT1, NSUN4, NSUN5, and ALYREF were upregulated in lung adenocarcinoma (52). Zhang et al. (53) have found that DNMT1 could enhance the radiosensitivity of HPV-positive HNSCC through suppressor of morphogenesis in genitalia 1 (SMG1). The NSUN family member, NSUN2, is found to be implicated in regulating cell cycles and accumulates in a variety of tumor lesions compared with normal samples (25, 54). Therefore, it seems valuable to further investigate the role of m⁵C regulators in HNSCC.

This study attempted to uncover the prognostic effects of m⁵C RNA methylation regulators in HNSCC. We identified four subgroups of HNSCC based on m⁵C RNA methylation regulators by means of consensus clustering and found that the classification was related to OS and tumor grade, implying that the expression pattern of m⁵C-related genes was positively correlated with the malignant process and prognosis of HNSCC. A previous study divided the TCGA HNSCC cohort into two subgroups depending on the 13 m⁶A RNA methylation regulators and applied consensus clustering (55). The OS and tumor grade of patients were also found to be strongly different between two identified subgroups. Taken together, these results indicate that RNA methylation regulators (m⁶A, m⁵C) might be associated with the prognosis of HNSCC.

Afterward, the prognostic three-gene risk signature, comprised of NSUN5, DNMT1, and DNMT3A, was built to effectively distinguish between patients with high risk and patients with low risk and robustly predict OS in the subgroups of HNSCC. Specifically, we observed that low NSUN5, together with high DNMT1 and DNMT3A levels, were positively associated with favorable functional outcomes in patients with HNSCC. Besides, both univariate and multivariate Cox analyses revealed that the risk score could act as an independent prognostic factor in HNSCC, implying that NSUN5, DNMT1, and DNMT3A could be involved in tumor oncogenes and suppressors. It is noteworthy that all three risk genes were m⁶A “writers,” though they showed the opposite effect on survival in HNSCC, which may hint that the NSUN and DNMT family proteins affect the OS in HNSCC with diverse impacts.

Currently, the roles of NSUN5, DNMT1, and DNMT3A in tumors have been widely explored. NSUN5 was reported to be responsible for modifying the second m⁵C position in eukaryotic rRNA and maintaining global protein synthesis and normal growth in mice (56, 57). Another study demonstrated that NSUN5 epigenetic inactivation is a hallmark of long-term survival for patients with glioma (50). Additionally, emerging research has indicated that DNMT1 and DNMT3A share oncogene affections (58, 59). These results agree well with the constructed risk model using the three risk genes. Taken together, the risk model offers a basis for further studies of pathogenesis, and for the determination of the novel classification and construction of the prognosis model of HNSCC.

Moreover, the m⁵C-related risk model was found to be associated with the signaling pathways and biological functions of HNSCC. The role of m⁵C RNA regulators was discovered not so long ago. Thus, we identified several functional annotations and signaling pathways related to the two risk subgroups. In the high-risk subgroup, the translation, protein binding, and oxidative phosphorylation terms were enriched by the Reactome, the GO, and the KEGG databases, which were ascertained to be correlated with tumorigenesis (60–62). On the other hand, the top 20 DEGs identified in the two risk groups are probably associated with a high likelihood of underregulation by the three m⁵C writers (NSUN5, DNMT1, and DNMT3A). In detail, the results suggest an m⁵C-regulated mechanism in HNSCC, by which NSUN5 might target KRT83, KRT79, and BPIFA1, while DNMT1 and DNMT3A methylate PRR27, NPS, and STATH. Consistent with these m⁵C “writers,” m⁵C “readers” ALYREF and YBX1 might be implicated in HNSCC by functioning as mediators of the mRNA output from the nucleus and by maintaining the stability of their target mRNAs (12, 21). Based on this, the hypotheses on the role of m⁵C in HNSCC initiation and progression through its internal interactions and signaling pathways have been proposed (Figure 9). Besides, the IF analysis revealed that NSUN5, DNMT1, and DNMT3A were similarly overexpressed at the protein level, which was concordant with the transcriptomic level. This is the first characteristic for detecting RNA methylation regulators using IF. These findings, if verified in a larger cohort of prospective clinical cases combined with prognostic data, might be precise for the forecasting and management of patients with HNSCC. In addition to this shortcoming, we also need to acknowledge another limitation that we came across. Whether m⁵C decoration on RNA was associated with HNSCC prognosis was not directly proven.

FIGURE 9

Figure 9. The presumptive outline of m⁵C RNA modification in HNSCC. Methyltransferases of the NSUN family and DNMT family are the main m⁵C “writers.” Among them, NSUN5 is predicted to methylate KRT83, KRT79, and BPIFA1; whereas DNMT1 and DNMT3A target PRR27, NPS, and STATH. TET2 and TET3 serve as m⁵C “erasers,” which regulate RNA m⁵C dynamics together with “writers” in the nucleus. m⁵C readers bind m⁵C-modified RNA and modulate the processes of mRNA translation, alternative splicing, and RNA decay. Specifically, ALYREF functions in the nucleus and regulates mRNA export, while YBX1 maintains mRNA stability in the cytoplasm.

Conclusion

In summary, we systematically illustrated the expression profile, biological function, and clinical prognostic value of m⁵C regulators in HNSCC. The association between m⁵C-related genes and HNSCC progression has been identified. Furthermore, a 3-m⁵C-related gene-based risk score model was built using NSUN5, DNMT1, and DNMT3A, hinting at a prognostic value in HNSCC.

Data Availability Statement

The publicly available datasets were analyzed in this study. All of the raw data can be found in The Cancer Genome Atlas (TCGA), the Gene Expression Profiling Interactive Analysis 2 (GEPIA2), and the cBioPortal databases.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee of Shanghai Ninth People's Hospital affiliated to Shanghai Jiao Tong University, School of Medicine. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

ZT conceived and developed the outline of the study. ZH and BY downloaded the data and performed the statistical study. ZH and XZ generated figures and tables and wrote the manuscript. YW revised the paper. All of the authors read and approved the final manuscript.

Funding

This work was supported by the Natural Science Foundation of Shanghai (Grant No. 18ZR1422200).

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.

Acknowledgments

The authors would like to thank Yuanming Lei and Qin Wang for their valuable assistance and instructions to improve this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.592107/full#supplementary-material

References

1. Zhou C, Parsons JL. The radiobiology of HPV-positive and HPV-negative head and neck squamous cell carcinoma. Expert Rev Mol Med. (2020) 22:e3. doi: 10.1017/erm.2020.4

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Huang F, Xin C, Lei K, Bai H, Li J, Chen Q. Non-coding RNAs in oral premalignant disorders and oral squamous cell carcinoma. Cell Oncol (Dordr). (2020) 43:763–77. doi: 10.1007/s13402-020-00521-9

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Qi Y, Li W, Kang S, Chen L, Hao M, Wang W, et al. Expression of BDNF, TrkB, VEGF and CD105 is associated with pelvic lymph node metastasis and prognosis in IB2-stage squamous cell carcinoma. Exp Ther Med. (2019) 18:4221–30. doi: 10.3892/etm.2019.8100

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ahn PH, Machtay M, Anne PR, Cognetti D, Keane WM, Wuthrick E, et al. Phase I trial using induction ciplatin, docetaxel, 5-FU and erlotinib followed by cisplatin, bevacizumab and erlotinib with concurrent radiotherapy for advanced head and neck cancer. Am J Clin Oncol. (2018) 41:441–6. doi: 10.1097/COC.0000000000000317

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Argiris A, Li S, Savvides P, Ohr JP, Gilbert J, Levine MA, et al. Phase III Randomized Trial of Chemotherapy With or Without Bevacizumab in Patients With Recurrent or Metastatic Head and Neck Cancer. J Clin Oncol. (2019) 37:3266–74. doi: 10.1200/JCO.19.00555

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Zhou X, Han J, Zhen X, Liu Y, Cui Z, Yue Z, et al. Analysis of genetic alteration signatures and prognostic values of m6A regulatory genes in head and neck squamous cell carcinoma. Front Oncol. (2020) 10:718. doi: 10.3389/fonc.2020.00718

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Deng X, Jiang Q, Liu Z, Chen W. Clinical significance of an m6A reader gene, IGF2BP2, in head and neck squamous cell carcinoma. Front Mol Biosci. (2020) 7:68. doi: 10.3389/fmolb.2020.00068

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Barbieri I, Kouzarides T. Role of RNA modifications in cancer. Nat Rev Cancer. (2020) 20:303–22. doi: 10.1038/s41568-020-0253-2

CrossRef Full Text | Google Scholar

9. Yao RW, Wang Y, Chen LL. Cellular functions of long non-coding RNAs. Nat Cell Biol. (2019) 21:542–51. doi: 10.1038/s41556-019-0311-8

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Gebert LFR, MacRae IJ. Regulation of microRNA function in animals. Nat Rev Mol Cell Biol. (2019) 20:21–37. doi: 10.1038/s41580-018-0045-7

CrossRef Full Text | Google Scholar

11. Zhang J, Yang PL, Gray NS. Targeting cancer with small molecule kinase inhibitors. Nat Rev Cancer. (2009) 9:28–39. doi: 10.1038/nrc2559

CrossRef Full Text | Google Scholar

12. Yang X, Yang Y, Sun BF, Chen YS, Xu JW, Lai WY, et al. 5-methylcytosine promotes mRNA export - NSUN2 as the methyltransferase and ALYREF as an m(5)C reader. Cell Res. (2017) 27:606–25. doi: 10.1038/cr.2017.55

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Bohnsack KE, Höbartner C, Bohnsack MT. Eukaryotic 5-methylcytosine (m5C) RNA Methyltransferases: Mechanisms, Cellular Functions, and Links to Disease. Genes. (2019) 10:102. doi: 10.3390/genes10020102

CrossRef Full Text | Google Scholar

14. Boccaletto P, Machnicka MA, Purta E, Piatkowski P, Baginski B, Wirecki TK, et al. MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res. (2018) 46:D303–7. doi: 10.1093/nar/gkx1030

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Andersen NM, Douthwaite S. YebU is a m5C methyltransferase specific for 16 S rRNA nucleotide 1407. J Mol Biol. (2006) 359:777–86. doi: 10.1016/j.jmb.2006.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ban Y, Tan P, Cai J, Li J, Hu M, Zhou Y, et al. LNCAROD is stabilized by m6A methylation and promotes cancer progression via forming a ternary complex with HSPA1A and YBX1 in head and neck squamous cell carcinoma. Mol Oncol. (2020) 14:1282–96. doi: 10.1002/1878-0261.12676

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Trixl L, Lusser A. The dynamic RNA modification 5-methylcytosine and its emerging role as an epitranscriptomic mark. Wiley Interdiscip Rev RNA. (2019) 10:e1510. doi: 10.1002/wrna.1510

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Zhang X, Liu Z, Yi J, Tang H, Xing J, Yu M, et al. The tRNA methyltransferase NSun2 stabilizes p16INK4 mRNA by methylating the 3'-untranslated region of p16. Nat Commun. (2012) 3:712. doi: 10.1038/ncomms1692

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Dor Y, Cedar H. Principles of DNA methylation and their implications for biology and medicine. Lancet. (2018) 392:777–86. doi: 10.1016/S0140-6736(18)31268-6

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Kouzarides T. Chromatin modifications and their function. Cell. (2007) 128:693–705. doi: 10.1016/j.cell.2007.02.005

CrossRef Full Text | Google Scholar

21. Chen X, Li A, Sun BF, Yang Y, Han YN, Yuan X, et al. 5-methylcytosine promotes pathogenesis of bladder cancer through stabilizing mRNAs. Nat Cell Biol. (2019) 21:978–90. doi: 10.1038/s41556-019-0361-y

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Fu L, Guerrero CR, Zhong N, Amato NJ, Liu Y, Liu S, et al. Tet-mediated formation of 5-hydroxymethylcytosine in RNA. J Am Chem Soc. (2014) 136:11582–5. doi: 10.1021/ja505305z

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Yang Y, Wang L, Han X, Yang WL, Zhang M, Ma HL, et al. RNA 5-Methylcytosine Facilitates the Maternal-to-Zygotic Transition by Preventing Maternal mRNA Decay. Mol Cell. (2019) 75:1188–202.e11. doi: 10.1016/j.molcel.2019.06.033

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Flores JV, Cordero-Espinoza L, Oeztuerk-Winder F, Andersson-Rolf A, Selmi T, Blanco S, et al. Cytosine-5 RNA methylation regulates neural stem cell differentiation and motility. Stem Cell Reports. (2017) 8:112–24. doi: 10.1016/j.stemcr.2016.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Frye M, Watt FM. The RNA methyltransferase Misu (NSun2) mediates Myc-induced proliferation and is upregulated in tumors. Curr Biol. (2006) 16:971–81. doi: 10.1016/j.cub.2006.04.027

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Harris T, Marquez B, Suarez S, Schimenti J. Sperm motility defects and infertility in male mice with a mutation in Nsun7, a member of the Sun domain-containing family of putative RNA methyltransferases. Biol Reprod. (2007) 77:376–82. doi: 10.1095/biolreprod.106.058669

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Jhiang SM, Yaneva M, Busch H. Expression of human proliferation-associated nucleolar antigen p120. Cell Growth Differ. (1990) 1:319–24.

PubMed Abstract | Google Scholar

28. Ning X, Shi Z, Liu X, Zhang A, Han L, Jiang K, et al. DNMT1 and EZH2 mediated methylation silences the microRNA-200b/a/429 gene and promotes tumor progression. Cancer Lett. (2015) 359:198–205. doi: 10.1016/j.canlet.2015.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Li C, Wang S, Xing Z, Lin A, Liang K, Song J, et al. A ROR1-HER3-lncRNA signalling axis modulates the Hippo-YAP pathway to regulate bone metastasis. Nat Cell Biol. (2017) 19:106–19. doi: 10.1038/ncb3464

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Mounir M, Lucchetta M, Silva TC, Olsen C, Bontempi G, Chen X, et al. New functionalities in the TCGAbiolinks package for the study and integration of cancer data from GDC and GTEx. PLoS Comput Biol. (2019) 15:e1006701. doi: 10.1371/journal.pcbi.1006701

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. (2012) 2:401–4. doi: 10.1158/2159-8290.CD-12-0095

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. (2019) 47:W556–60. doi: 10.1093/nar/gkz430

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. (2013) 6:pl1. doi: 10.1126/scisignal.2004088

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. (2019) 47:D607–13. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Jassal B, Matthews L, Viteri G, Gong C, Lorente P, Fabregat A, et al. The reactome pathway knowledgebase. Nucleic Acids Res. (2020) 48:D498–503. doi: 10.1093/nar/gkz1031

CrossRef Full Text | Google Scholar

38. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. (2000) 25:25–9. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Huang da W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. (2009) 37:1–13. doi: 10.1093/nar/gkn923

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. (2009) 4:44–57. doi: 10.1038/nprot.2008.211

CrossRef Full Text | Google Scholar

41. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Qin T, Koneva LA, Liu Y, Zhang Y, Arthur AE, Zarins KR, et al. Significant association between host transcriptome-derived HPV oncogene E6^* influence score and carcinogenic pathways, tumor size, and survival in head and neck cancer. Head Neck. (2020) 42:2375–89. doi: 10.1002/hed.26244

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Plzák J, Bouček J, Bandúrová V, Kolár M, Hradilová M, Szabo P, et al. The head and neck squamous cell carcinoma microenvironment as a potential target for cancer therapy. Cancers. (2019) 11:440. doi: 10.3390/cancers11040440

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Chernyakov I, Whipple JM, Kotelawala L, Grayhack EJ, Phizicky EM. Degradation of several hypomodified mature tRNA species in Saccharomyces cerevisiae is mediated by Met22 and the 5'-3' exonucleases Rat1 and Xrn1. Genes Dev. (2008) 22:1369–80. doi: 10.1101/gad.1654308

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Courtney DG, Tsai K, Bogerd HP, Kennedy EM, Law BA, Emery A, et al. Epitranscriptomic addition of m(5)C to HIV-1 transcripts regulates viral gene expression. Cell Host Microbe. (2019) 26:217–27.e6. doi: 10.1016/j.chom.2019.07.005

CrossRef Full Text | Google Scholar

46. Hoernes TP, Clementi N, Faserl K, Glasner H, Breuker K, Lindner H, et al. Nucleotide modifications within bacterial messenger RNAs regulate their translation and are able to rewire the genetic code. Nucleic Acids Res. (2016) 44:852–62. doi: 10.1093/nar/gkv1182

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Sharma S, Lafontaine DLJ. ‘View from a bridge': a new perspective on eukaryotic rRNA base modification. Trends Biochem Sci. (2015) 40:560–75. doi: 10.1016/j.tibs.2015.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Guo G, Wang H, Shi X, Ye L, Yan K, Chen Z, et al. Disease activity-associated alteration of mRNA m(5) C methylation in CD4(+) T cells of systemic lupus erythematosus. Front Cell Dev Biol. (2020) 8:430. doi: 10.3389/fcell.2020.00430

CrossRef Full Text | Google Scholar

49. Xue M, Shi Q, Zheng L, Li Q, Yang L, Zhang Y. Gene signatures of m5C regulators may predict prognoses of patients with head and neck squamous cell carcinoma. Am J Transl Res. (2020) 12:6841–52.

PubMed Abstract | Google Scholar

50. Janin M, Ortiz-Barahona V, de Moura MC, Martínez-Cardús A, Llinàs-Arias P, Soler M, et al. Epigenetic loss of RNA-methyltransferase NSUN5 in glioma targets ribosomes to drive a stress adaptive translational program. Acta Neuropathol. (2019) 138:1053–74. doi: 10.1007/s00401-019-02062-4

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Xu J, Ji L, Liang Y, Wan Z, Zheng W, Song X, et al. CircRNA-SORE mediates sorafenib resistance in hepatocellular carcinoma by stabilizing YBX1. Signal Transduct Target Ther. (2020) 5:298. doi: 10.1038/s41392-020-00375-5

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Sun L, Liu WK, Du XW, Liu XL, Li G, Yao Y, et al. Large-scale transcriptome analysis identified RNA methylation regulators as novel prognostic signatures for lung adenocarcinoma. Ann Transl Med. (2020) 8:751. doi: 10.21037/atm-20-3744

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Zhang C, Mi J, Deng Y, Deng Z, Long D, Liu Z. DNMT1 enhances the radiosensitivity of HPV-positive head and neck squamous cell carcinomas via downregulating SMG1. Onco Targets Ther. (2020) 13:4201–11. doi: 10.2147/OTT.S227395

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Hussain S, Benavente SB, Nascimento E, Dragoni I, Kurowski A, Gillich A, et al. The nucleolar RNA methyltransferase Misu (NSun2) is required for mitotic spindle stability. J Cell Biol. (2009) 186:27–40. doi: 10.1083/jcb.200810180

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Zhao X, Cui L. Development and validation of a m(6)A RNA methylation regulators-based signature for predicting the prognosis of head and neck squamous cell carcinoma. Am J Cancer Res. (2019) 9:2156–69.

Google Scholar

56. Heissenberger C, Liendl L, Nagelreiter F, Gonskikh Y, Yang G, Stelzer EM, et al. Loss of the ribosomal RNA methyltransferase NSUN5 impairs global protein synthesis and normal growth. Nucleic Acids Res. (2019) 47:11807–25. doi: 10.1093/nar/gkz1043

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Schosserer M, Minois N, Angerer TB, Amring M, Dellago H, Harreither E, et al. Methylation of ribosomal RNA by NSUN5 is a conserved mechanism modulating organismal lifespan. Nat Commun. (2015) 6:6158. doi: 10.1038/ncomms7158

CrossRef Full Text | Google Scholar

58. Xu SF, Zheng Y, Zhang L, Wang P, Niu CM, Wu T, et al. Long non-coding RNA LINC00628 interacts epigenetically with the LAMA3 promoter and contributes to lung adenocarcinoma. Mol Ther Nucleic Acids. (2019) 18:166–82. doi: 10.1016/j.omtn.2019.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Zhang J, Yang C, Wu C, Cui W, Wang L. DNA methyltransferases in cancer: biology, paradox, aberrations, and targeted therapy. Cancers. (2020) 12:2123. doi: 10.3390/cancers12082123

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Diamantopoulou Z, Castro-Giner F, Aceto N. Circulating tumor cells: Ready for translation? J Exp Med. (2020) 217:e20200356. doi: 10.1084/jem.20200356

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Ghosh P, Guo Y, Ashrafi A, Chen J, Dey S, Zhong S, et al. Oxygen-enhanced optoacoustic tomography reveals the effectiveness of targeting heme and oxidative phosphorylation at normalizing tumor vascular oxygenation. Cancer Res. (2020) 80:3542–55. doi: 10.1158/0008-5472.CAN-19-3247

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Kang YT, Purcell E, Palacios-Rolston C, Lo TW, Ramnath N, Jolly S, et al. Isolation and profiling of circulating tumor-associated exosomes using extracellular vesicular lipid-protein binding affinity based microfluidic device. Small. (2019) 15:e1903600. doi: 10.1002/smll.201903600

PubMed Abstract | CrossRef Full Text | Google Scholar