Bioinformatics identification and validation of m6A/m1A/m5C/m7G/ac4 C-modified genes in oral squamous cell carcinoma

Abstract

Background: RNA modifications, including m6A, m1A, m5C, m7G, and ac4C, may play a role in the occurrence and development of cancer, such as proliferation. However, the effects of RNA modification-related genes (RRGs) in the development of oral squamous cell carcinoma (OSCC) have not been fully elucidated. The present study aimed to evaluate the effects and mechanisms of RRGs on OSCC development progression.

Methods: RNA-seq transcriptome data, along with clinical and prognostic information, were extracted for 328 patients with OSCC from the TCGA database. A total of 49 RRGs were analyzed for differential expression. We then performed Lasso analysis, as well as univariate and multivariate Cox regression analyses, followed by Kaplan-Meier survival analysis to identify relevant prognostic genes and establish a risk-prognosis model. Patients were categorized into high-risk and low-risk groups, and gene set enrichment analysis (GSEA) was conducted to analyze differences in gene signatures between these two groups, using data from the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) databases. RT-PCR was employed to validate the expression levels of differentially expressed genes in OSCC samples. The four most significantly differentially expressed genes were selected for further functional analysis, and small interfering RNA (siRNA) vectors targeting these genes were transfected into OSCC CAL27 cells. The Cell Counting Kit-8 (CCK-8) assay was used to evaluate cell proliferation. Additionally, a subcutaneous CAL27 xenograft model transfected with short hairpin RNA (shRNA), combined with Ki-67 immunohistochemical (IHC) staining and TUNEL assay, was used to investigate their underlying molecular mechanisms in vivo.

Results: Among the 49 RRGs, four genes (IGF2BP2, HNRNPC, NAT10, and TRMT61B) were found to be associated with the development of OSCC. Based on various methodological validations, a risk score model was constructed using these four genes. The high-risk and low-risk groups of OSCC patients exhibited significantly different survival outcomes and clinicopathological characteristics. Patients in the low-risk group had longer overall survival (OS) and lower mortality rates compared to those in the high-risk group. The nomogram and decision curve analysis (DCA) demonstrated that our risk model accurately and reliably predicted the impact of risk factors on OS at 1-, 3-, and 5-year. Additionally, risk scores correlated with the infiltration of several immune cells, particularly CD8⁺ T cells and B cells, which showed significant negative correlations. Furthermore, the results of the CCK-8 assay indicated that inhibition of NAT10 and IGF2BP2 expression using siRNA inhibited the proliferation of OSCC cell lines in vitro. Meanwhile, inhibition of NAT10 and IGF2BP2 expression using shRNA influenced proliferation of tumorigenicity in vivo.

Conclusion: In this study, we established a risk model and nomogram based on four RRGs, which can be used for risk stratification and predicting survival outcomes in patients with OSCC. This provides a reliable reference for individualized therapy in OSCC patients.

Keywords: Bioinformatics analysis; Biomarkers; Oral squamous cell carcinoma; Prognostic signature; RNA modification.

Fig. 1

Risk model based on RRGs. A LASSO coefficient path plot: The horizontal axis represents the log(λ) value, and the vertical axis shows the coefficients of candidate genes. Each colored line corresponds to one gene. B Partial likelihood deviance plot. The blue dotted line traces the partial likelihood deviance across varying log(λ). The error bars represent the standard deviation from 10-fold cross-validation. The two vertical dashed lines mark the λ values that yield the minimum deviance (left) and the simplest model (right) within one standard error of the minimum, respectively. C Forest plots for univariate (left) and multivariate (right) Cox analyses. Hazard ratios (HR) with 95% confidence intervals (CIs) are shown for each of the four selected genes (IGF2BP2, HNRNPC, NAT10, and TRMT61B). Genes with an HR > 1 indicate a higher risk of poor survival when overexpressed. p-values reflect the significance of each gene’s correlation with overall survival in oral squamous cell carcinoma (OSCC) patients

Fig. 2

Clinical characteristics of risk scores. A OncoPrint of IGF2BP2, HNRNPC, TRMT61B, and NAT10 using the cBioPortal database. Each column represents one patient sample, and each row corresponds to a gene. Different colors or bars indicate various genetic alterations (e.g., mutation, amplification, deep deletion). The bottom legend summarizes the frequency of each alteration type across the cohort. B, C Unpaired (B) and paired (C) differential analyses for the 4 genes between paraneoplastic and tumor tissues. The expression level of the 4 genes was higher in tumor tissues than in normal tissues. D Immunohistochemistry (IHC) of IGF2BP2, HNRNPC, TRMT61B, and NAT10 in normal versus tumor tissues. Representative tissue microarray images (from the Human Protein Atlas or in-house staining) show the protein localization and intensity. Brown staining indicates positive protein expression, while blue (hematoxylin) represents the nuclear counterstain. The upper panels depict normal oral mucosa, and the lower panels show corresponding OSCC tissues, confirming higher protein expression in tumor samples for most genes. Scale bars (if included) indicate the magnification used

Fig. 3

KM survival analysis, risk score assessment by the RRGs signature and time-dependent ROC curves in the TCGA OSCC cohort. A KM survival analysis of high- and low-risk samples of the TCGA OSCC training set. B ROC curves for OS of the training set. The AUC at 1-, 3-, and 5-year. C Risk score distribution in the training set, (D) survival status of the training set, and (E) 4 RRGs expression patterns for patients in high- and low- risk groups as determined by the four-RRGs signature in the training set. F KM survival analysis of high- and low-risk samples of the GSE85446 validation cohorts. G ROC curves for OS of the validation cohorts. The AUC at 1-, 3-, and 5-year. H Risk score distribution in the validation cohorts, (I) survival status of the validation cohorts, and (J) 4 RRGs expression patterns for patients in high- and low-risk groups based on the four-RRGs signature in the validation cohorts

Fig. 4

Univariate and multivariate Cox regression analysis of the 4-gene risk model, along with the establishment and assessment of a nomogram. A Univariate cox regression and multivariate cox regression analysis based the risk score and clinical characteristics. B The nomograms for predicting probabilities of OSCC patients OS. C The calibration curves. D The DCA plot

Fig. 5

Differences in survival rates between high- and low-risk group patients stratified by clinicopathological parameters in the validated OSCC cohort. A Male, Female. B Age ≤60 y, Age >60 y. C T1&2, T3, T4. D N0, N1, N2&3. The RRGs retained their prognosis prediction value in multiple subgroups analysis of OSCC patients

Fig. 6

Spearman's coefficient for the relationship between risk scores and the infiltration level of 6 immune cell types. A B cell, (B) Myeloid dendritic cell, (C) neutrophil, (D) CD8⁺T cell, (E) macrophage, and (F) CD4⁺T cell. The blue line in each graph corresponds to the linear model and represents the proportionate trend of tumor-infiltrating immune cell score (TICS) and risk score. The gray shading surrounding the blue line represent the 95% CI

Fig. 7

GSEA and KEGG/GO enrichment analysis showing the potential biological functions and signaling pathways of the 4 genes (IGF2BP2, HNRNPC, TRMT61B, and NAT10). A GESA analysis revealed that the 4 genes were significantly in critical biological functions and signal pathways which were correlated with cell circulation, drug metabolism, and glutathione metabolism. B The result of KEGG/GO analyses indicated that the 4 genes were primarily involved in critical biological functions and signaling pathways involved in drug metabolism, skin development, epidermal growth, and keratinocyte differentiation

Fig. 8

Inhibition of NAT10 and IGF2BP2 expression using siRNA inhibits proliferation in vitro. A The mRNA level of IGF2BP2, HNRNPC, NAT10, and TRMT61B in OSCC cell lines CAL27, WSU-HN30, WSU-HN6, SCC9, and HaCaT detected by RT-PCR. B IGF2BP2, HNRNPC, NAT10, and TRMT61B expression was decreased by IGF2BP2, HNRNPC, NAT10, and TRMT61B siRNAs in CAL27. C The proliferation potential of cells was assayed by CCK8 in CAL27 and WSU-HN30, and inhibition of NAT10 and IGF2BP2 suppressed the proliferation of CAL27 and WSU-HN30 cells

Fig. 9

Inhibition of NAT10 and IGF2BP2 expression using shRNA inhibits proliferation in vivo. A The sh-NAT10 and sh-IGF2BP2 CAL27 subcutaneous tumors showed smaller size compared with control tumors. B IHC staining in the mice subcutaneous tumors. sh-NAT10 and sh-IGF2BP2 CAL27 tumor showed lower staining of Ki-67 compared with control CAL27 tumors. Scale bar: 50 μm. C Analysis of TUNEL-positive cells in the mice subcutaneous tumors. sh-NAT10 and sh-IGF2BP2 CAL27 tumor showed higher staining of TUNEL compared with control CAL27 tumors. Scale bar: 50 μm