The tumor samples (primary OSCC tissue specimens and adjacent matched normal mucosa tissues) were obtained from 38 OSCC patients with or without lymph node metastasis who underwent surgery at the Department of Oral and Maxillofacial Surgery of the Second Xiangya Hospital. These patients underwent surgery between April 2011 and July 2013, and the follow-up period used for survival analyses ended in May 2018. No patients involved in this investigation received chemotherapy before surgery. Tumor histology and grading were classified according to the WHO guidelines. The details of the patients used in this study are present in Table 1 .

Immunohistochemistry (IHC) was performed to evaluate HNRNPA2B1 and LINE-1 expression using the standard protocol. The tissue sections were dewaxed with xylenes, hydrated with graded alcohols. Endogenous peroxidase was inhibited with 3% hydrogen peroxide. Antigen retrieval was performed using citrate buffer (pH 6.0) in a microwave oven. Anti-rabbit HNRNPA2B1 (1:50, ProteinTech, USA) and anti-mouse ORF1p (1:500, Sigma-Aldrich, Germany) primary antibodies were used to the sections overnight at 4°C, then incubated with a secondary antibody for an hour at room temperature. Signals were visualized by using a 3,30-diaminobenzidine (DAB) detection kit (DakoCytomation, Denmark). Sections were counterstained with hematoxylin, dehydrated, and mounted. The expressions of HNRNPA2B1 and ORF-1p are scored by two researchers. For each specimen, one score was assigned according to the percentage of positive cells: <25%: one point; 26%–50%: two points: 51%–75%, three points; and 76%–100%: four points. A second score was assigned according to the intensity of the staining, with negative staining equaling 0, weak staining equaling one point, moderate staining equaling two points, and strong staining equaling three points. IHC score was then calculated by multiplying the two scores described above. If the expression score was ≥5, the tissue was considered as positive expression.

The transcriptome profiles of 330 OSCC and 32 normal cases were downloaded from TCGA (https://tcga-data.nci.nih.gov/tcga/). TCGA is a publicly available dataset. No ethics approval was needed. Based on previous literatures, we have identified 19 m⁶A-related genes, including METTL3, METTL14, METTL16, WTAP, KIAA1429, RBM15, FTO, ALKBH5, ZC3H13, YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, HNRNPC, HNRNPA2B1, LRPPRC, IGF2BP1, and FMR1. Differential expression analysis was performed in R software, with the cutoff value of |log2(fold change [FC]|>2 and p <0.05. The heatmap and boxplot were constructed using the ggplot2 package in R software.

Cluster analysis was performed according to the m⁶A-related gene transcriptome profiles. Meanwhile, the result of cluster analysis was used to make a principal component analysis (PCA). The patient’s clinical data of survival time and survival status were extracted from 330 OSCC samples. The correlation analysis between clinical features and clustering results was performed through R. Finally, heatmaps were constructed through “pheatmap,” survival, and “survminer” packages.

In this study, we performed three survival analyses. For the first analysis, we divided all OSCC patients from the TCGA database into cluster 1 and cluster 2 based on the results of cluster analysis. For the second analysis, we divided OSCC patients into high-risk groups and low-risk groups according to the Cox regression model. For the third analysis, we divided the OSCC patient samples from the department or TCGA database into positive expression and negative expression, or high expression (with counts higher than the median) and low expression (with counts lower than the median). Then, we used the Kaplan–Meier method to analyze the candidate genes of significant prognostic value with p < 0.05.

First, the original expression profile data of m⁶A-related genes were normalized by log2(x + 1), and prognosis-associated factors were selected by univariate Cox regression. Next, we performed Cox regression analysis combined with LASSO regression to establish a risk regression and the penalty regularization parameter lambda (λ) was chosen through the cross-validation routine with an n-fold equal to 10 by using R package “glmnet” (32). Meanwhile, “lambda.min” was identified to pick out the variables. Finally, HNRNPA2B1, METTL3, YTHDF3, and HNRNPC were enrolled in risk Cox regression and survival analysis. Moreover, univariate and multivariate Cox regressions are performed through survival and “forestplot” packages to identify the independent prognostic factor of OSCC. Moreover, we evaluate the reliability of the risk regression through “survivalROC” package in R software.

To identify the key gene of m⁶A, we used the caret package to construct three different machine learning models to investigate features of importance involving 19 m⁶A-related genes. Neural network (nn), random forest (rf), and gradient boosting machine (gbm) were used to rank the 19 genes by their expression levels in two disease states (tumor and normal). The most important feature of the three models was selected as the key gene; results were visualized by ggplot2 package.

To provide a quantitative estimate of EMT, we constructed a model by multinominal logistic regression; each sample is given EMT scores, which range from 0 (pure epithelial) to 2 (pure mesenchymal), with a score of 1 indicating a maximal hybrid epithelial/mesenchymal (E/M) phenotype (33). The transcriptome profiles of OSCC used for GSEA were downloaded from TCGA (https://tcga-data.nci.nih.gov/tcga/) and GEO data (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138206). GSEA (34, 35) was performed with KEGG gene set (c2) or oncogenic signature gene set (c6) collections of the Molecular Signature Database v7.0(http://www.gsea-msigdb.org/gsea/downloads.jsp).

CAL27 and SCC4 cells were cultured in Dulbecco’s modified Eagle medium (DMEM) supplemented with 10% fetal bovine serum (Invitrogen, Germany), 2 mM L-glutamine (Invitrogen, Germany), and antibiotics (50 μg/ml penicillin–streptomycin, Invitrogen, Germany). Mouse anti-α-tubulin, rabbit anti-E-cadherin (24E10), and anti-N-cadherin (D4R1H) were from Cell Signaling Technology (Beverly, MA). Rabbit anti-Snail and anti-TGF-β1 were from Abcam (Cambridge, MA). Rabbit anti-HNRNPA2B1 was from ProteinTech (Rosemont, IL, USA), and mouse anti-ORF1p was purchased from Sigma-Aldrich (clone 4H1, Germany).

Commercial shRNA plasmid vectors (U6–MCS–Ubiquitin–Cherry–IRES–puromycin) carrying target-specific ENAuences against human HNRNPA2B1 and non-target scrambled control were purchased from GeneChem (Shanghai, China). Lentiviral particles were produced in HEK293T packaging cells with shRNA, psPAX2 (#12260, Addgene), and pMD2.G (#12259, Addgene) plasmids. Then, CAL27 and SCC4 cells were infected with 0.45-μm-pore filtered viral supernatants. Infected cells of CAL27 and SCC4 were selected for 2 weeks using 1 μg/ml puromycin. The expression of HNRNPA2B1 was confirmed by RT-PCR and Western blot. Oligonucleotide sequences for shRNAs are the following: HNRNPA2B1 shRNA#1 and shRNA#2 are 5′-TGACAACTATGGAGGAGGAAA-3′ and 5′-AGAAGCTGTTTGTTGGCGGAA-3′, respectively. The sequence of scrambled shRNA is 5′-TTCTCCGAACGTGTCACGT-3′.

The HNRNPA2B1 overexpression lentiviral vector-plasmid (GV208, Ubi-MCS-HNRNPA2B1-EGFP) was purchased from GeneChem (Shanghai, China) and used to transfect the HEK293T packaging cells. CAL27 and SCC4 cells were transfected with lentivirus containing HNRNPA2B1 cDNA for subsequent analysis.

Total RNA was isolated using TRIzol Reagent (Invitrogen, Thermo Fisher Scientific, USA) and reserved transcribed with TaqMan Reverse Transcription Reagents Kit (Applied Biosystems, Foster City CA). The resulting cDNAs (5 μl) were used as templates for qRT-PCR with specific primers including HNRNPA2B1 (forward: 5′-TAGTACTACAGGTTGCTGCC-3′, reverse: 5′-CTGTGAGGCTAGACTAGTAGT-3′) and GAPDH (forward: 5′-ACAACTTTGGTATCGTGGAAGG-3′, reverse: 5′-GCCATCACGCCACAGTTTC-3). qRT-PCR was carried out with Power SYBR Green PCR Master Mix. The qRT-PCR cycling conditions were as follows: 95°C for 10 s; 40 cycles of 95°C for 5 s, 60°C for 20 s, 95°C for 60 s, 55°C for 30 s, and 95°C for 30 s. Fold changes were determined by comparing the ΔCT value of each product normalized to GAPDH as an internal control.

The whole-cell lysates were prepared using RIPA lysis buffer according to the standard protocol. For Western blot, 15 μg protein for each group was resolved by SDS-PAGE and transferred onto polyvinyldinefluoride (PVDF) membranes (Millipore, Billerica, MA). The membranes were blocked with appropriate primary antibodies indicated above. After incubating with secondary antibody, the membranes were visualized with enhanced chemiluminescence and the bands were analyzed with ImageJ software.

Cell proliferation was measured using Cell Counting Kit-8 (CCK8, Bimake, Houston, TX, USA) assay. CAL27 and SCC4 cells were seeded in 96-well plates, three wells for each group, with each well containing 2*10³ cells. At 24, 48, 72, 96, and 120 h, 10 μl of CCK8 reagent was added to each well and the cells were incubated for 1 h at 37°C. Then the absorbance at 450 nm was measured using the microplate reader.

CAL27 and SCC4 shRNA stable cells were plated in six-well plates at a density of 2.5 × 10⁶ cells/well. After the cells reached 100% confluency, straight wounds were made by using 200-μl tips and the cells were incubated in medium with 1% FBS. The healing of the gap was captured at 0, 24, 48, and 72 h. Experiments were repeated at least three times.

Cell invasion assays were performed in 24-well plates using chambers with 8-μm pores (Corning Incorporated, Kennebunk, ME, USA) coated with an indicator layer of growth factor reduced Matrigel. 5 × 10⁴ CAL27 and SCC4 stable cells were seeded in the upper chamber under serum-free conditions, and DMEM with 10% FBS was placed in the lower chamber. After 24 h, cells on the upper chamber were scraped and washed away; cells on the lower chambers were fixed in 4% formaldehyde and stained with 0.5% crystal violet for 20 min. Next, the invasion cell numbers of each field were counted via microscope.

TRIzol (Thermo Fisher Scientific, Waltham, MA, USA) was used to extract total RNA from OSCC patients’ tissue samples. m⁶A methylation was quantified according to the guideline of m⁶A RNA Methylation Assay Kit (ab185912, Abcam, Cambridge, MA).

RIP was performed with Magna RIP RNA-binding protein immunoprecipitation kits (Millipore, Billerica, MA). Antibodies against HNRNPA2B1 were used. Total RNA from cells was extracted and depleted of ribosomal RNA. The RNA protein complexes were washed and mixed with 900 μl RIP Immunoprecipitation Buffer. Next, RNA was purified with 150 μl proteinase K buffer. Finally, RNAs were extracted and evaluated by qRT-PCR, which was normalized to input and IgG.

The primers of LINE-1 are forward: 5′-GGGCTGCACCTGTCAAGATA-3′, reverse: 5′-ACCTTGCCATCTTTTCCCGT-3′.

All analyses were performed with the SPSS 25.0 software (SPSS Inc., USA), R software (Version 4.0.3), and GraphPad Prism software version 7.0 (GraphPad Software, Inc., La Jolla, CA, USA). If there are no special instructions, all bioinformatics analysis and visualization in this study are performed using R software. Statistical significance was determined by Student’s t-test. The analysis of overall survival and progression-free survival time was calculated using the Kaplan–Meier method, and the differences in survival between the groups were compared using a log-rank test. The effect of clinicopathological factors on survival was determined with univariate and multivariate Cox proportional hazard models. Correlations between HNRNPA2B1 and EMT marker expression levels were computed using Pearson correlation analysis. Data from three independent experiments were presented as mean ± SD. A p-value of <0.05 was considered statistically significant. *, **, *** and **** indicated p < 0.05, p < 0.01, p < 0.001 and p < 0.001, respectively.

More and more evidence showed that m⁶A RNA methylation promotes tumor initiation and progression. However, there are fewer studies to explore the role of m⁶A in OSCC. To study the role of m⁶A regulatory genes in OSCC, we detect the m⁶A level in OSCC and normal adjacent tissues. The result showed that the m⁶A level upregulated in tumor tissues compared with normal adjacent tissues ( Figure 1B ). Further, we evaluated all 19 m⁶A-related genes in 330 OSCC samples and 32 normal samples from the TCGA database, including METTL3, METTL14, METTL16, WTAP, KIAA1429, RBM15, FTO, ALKBH5, ZC3H13, YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, HNRNPC, HNRNPA2B1, LRPPRC, IGF2BP1, and FMR1. The expression profiles of these genes were extracted from the transcriptome data, and it was found that 13 genes were abnormally expressed in tumor tissues compared with control ( Figures 1A, C ).

Cluster analysis was performed to analyze the 330 OSCC samples from the TCGA database. Based on the expression similarity of m⁶A-related genes, k = 2 was demonstrated to be the most appropriated selection to divide the OSCC patient cohort into two clusters, namely, cluster 1 and cluster 2 ( Figures 2A, B ). Subsequently, the overall survival rate of the two clusters estimated by Kaplan–Meier was significantly different ( Figure 2C ). On the basis of correlation analysis of clinical characteristics, an obvious difference was found between cluster 1 and cluster 2 for the metastasis, stage, and T classification, while no significant difference was observed for other parameters such as age, gender, and grade ( Figure 2D ).

Univariate and multivariate Cox regression analyses were used to identify the m⁶A-related genes that are associated with OS in the TCGA OSCC cohort. The results show that METTL3, YTHDF3, HNRNPC, and HNRNPA2B1 might be independent biomarkers in OSCC ( Figures 3A, B ). Furthermore, LASSO Cox regression was performed to determine variables ( Figures 3C, D ). Finally, four variables including METTL3, YTHDF3, HNRNPC, and HNRNPA2B1 were selected in Cox regression. The risk score for each patient was calculated with the following formula: risk scores = (0.1378*HNRNPA2B1) + (0.0036*METTL3) + (0.0021*YTHDF3) + (0.0095*HNRNPC). Then, all patients were grouped into the high-risk group and low-risk group, respectively. The survival analysis showed that the OSCC patients in the high-risk group had a significantly shorter overall survival than those in the low-risk group ( Figure 3F ). The prognostic signature model showed good prediction efficiency with the area under the ROC curve (AUC) value equal to 0.632 ( Figure 3E ). There were significant differences between the high-risk and low-risk groups with respect to the metastasis and T classification ( Figure 3I ). Here, we found that risk score may be an independent prognostic risk factor for OSCC ( Figures 3G, H ).

To evaluate the important characteristics of the 19 m⁶A-related genes, we used the neural network (nn), random forest (rf), and gradient boosting machine (gbm) to build predictive models, respectively ( Figures 4A–C ). We found that the HNRNPA2B1 gene was consistently ranked the highest among the 19 genes in the three models, and it was selected as the most prognostically important locus of the 19 m⁶A regulatory genes in OSCC. Combining with the results of the previous univariate and multivariate Cox regression analyses, HNRNPA2B1 was selected to conduct further research. After analyzing the relationship between HNRNPA2B1 and each clinical parameter from the TCGA database, we found that HNRNPA2B1 was overexpressed in tumors ( Figure 4D ) and HNRNPA2B1 expression was significantly correlated with tumor stage, T-class, and metastasis ( Figures 4E–G ) but without grade ( Figure 4H ). Moreover, the mRNA expression of HNRNPA2B1 was examined in eight patients’ fresh cervical lymph nodes samples (four lymph node metastasis and four non-lymph node metastasis) by real-time PCR. The results showed that HNRNPA2B1 mRNA levels in lymph node metastasis patients were generally much higher than those in non-lymph node metastasis patients (the mean CT values of tumor samples with or without lymph node metastasis were 15.2572 and 16.5325, respectively. The mean ΔCT values were -0.675 and 0.1675, respectively, and the mean fold changes were 1.794 and 0.829, respectively) ( Figure 4I ). Moreover, to evaluate the expression of HNRNPA2B1 in OSCC, we examined HNRNPA2B1 expression of 38 OSCC patient tissues from our department by immunohistochemistry (IHC) staining. The expression status of HNRNPA2B1 was observed mainly in the nucleus of cells. After the IHC was scored, the results indicated lower HNRNPA2B1 expression in normal cases and higher expression in tumor tissues and metastatic lymph node ( Figures 4J, K ). These results were consistent with the HNRNPA2B1 mRNA expression results from the TCGA database ( Figure 4D ). Fold changes were determined by comparing the ΔCT value of each product normalized to GAPDH as an internal control.

According to the follow-up data of the 38 OSCC patients from our department, we found that patients with high HNRNPA2B1 expression had a significantly poor 5-year OS compared with those with low expression of HNRNPA2B1 (p < 0.05, Figure 5A ). The survival analysis from TCGA further validated this result (p < 0.05, Figure 5B ). The univariate and multivariate Cox proportional hazards regression analyses revealed that HNRNPA2B1 expression was an independent prognostic factor for poor OS (p < 0.001, Table 1 ). Interestingly, the relationship between HNRNPA2B1 expression and the clinicopathological characteristics was statistically analyzed, and the results showed a direct association between HNRNPA2B1 expression and histological grade and metastasis ( Table 2 ). Here, the definition of smoking is as follows: participants who smoked more than 100 cigarettes were classified as smokers (36). The definition of drinking is as follows: the participant had drunk alcohol regularly (i.e., drank at least once a week on a regular basis) during the past 12 months, otherwise were considered non-drinkers (37).

CAL27 cell and SCC4 cell were selected to establish HNRNPA2B1 knockdown and overexpress stable cell lines by lentiviral vectors. Western blot was used to measure the knockdown efficacy ( Figure 6A ).

The CCK8 assay indicated that downregulation of HNRNPA2B1 significantly inhibited cell growth compared with the control (**p < 0.01) ( Figures 6B, C ). Wound healing assay was performed to determine the role of HNRNPA2B1 in regulating OSCC cell migration. After knocking down HNRNPA2B1#1 and HNRNPA2B1#2, the cell migration ability was remarkably suppressed ( Figures 6D, E ). In addition, Transwell assay showed a significant decrease in the HNRNPA2B1 knockdown group ( Figure 6F ). We also found that overexpression of HNRNPA2B1 could enhance the proliferation, migration, and invasion ability of OSCC cells ( Figures 6B–F ). Together, we found that HNRNPA2B1 was positively correlated with OSCC progression.

As a tumor-related gene, HNRNPA2B1 is associated with tumor progression (30). To investigate the mechanisms, we divided the OSCC patients into two groups (HNRNPA2B1 low expression group and HNRNPA2B1 high expression group). Through GSEA, we found that the top upregulated gene sets in the HNRNPA2B1 high expression group were related to EMT, such as tight junction, focal adhesion, adherens junction, and TGF-β signaling pathway ( Figures 7B, C ). Moreover, the HNRNPA2B1 high expression group had a higher EMT score ( Figure 7A ) and multiple EMT markers were differential expression in the HNRNPA2B1 high and low expression group ( Figures 7D–G ). Epithelial–mesenchymal transition (EMT) is always observed at invasive OSCC and is significantly correlated with metastasis in tumor progression (38). In this study, we showed that reduced expression of HNRNPA2B1 in OSCC cells leads to decreased proliferation, migration, and invasion. These results illuminate that HNRNPA2B1 may promote OSCC progression by regulates EMT.

The TGF-β signaling pathway is a classic pathway that promotes EMT (39), GSEA indicated that the TGF-β signaling pathway in the HNRNPA2B1 high expression group was stimulated ( Figures 7B, C ). Our previous research demonstrated that LINE-1 may mediate the EMT process through the TGF-β1/Smad2/Snail signaling pathway in OSCC. HNRNPA2B1, as an RNA-binding protein, participates in the translation of many proteins (40–42). To explore whether HNRNPA2B1 plays a role in LINE-1 expression and whether HNRNPA2B1 mediates the EMT process by regulating the expression of LINE-1, we initially analyzed the transcription and translation of LINE-1 and HNRNPA2B1. The mRNA expression profile data of HNRNPA2B1 and LINE-1 from the TCGA database and show that their mRNA expression levels were not correlated ( Figure 7H ). Through the immunohistochemistry results of patient tumor tissues of 38 OSCC patients, we found that the protein levels of HNRNPA2B1 and LINE-1 are significantly positively correlated ( Figure 4L ).

Then, we examined the transcription and translation of LINE-1 upon HNRNPA2B1 loss. As expected, HNRNPA2B1 overexpression increased the protein abundance of LINE-1 but did not affect its RNA level in both CAL27 and SCC4 cells ( Figures 7I, J ). Furthermore, we detected the EMT-related markers as indicated by Western blot after downregulating and upregulating the expression of HNRNPA2B1 in CAL27 and SCC4 cells ( Figure 7J ). The results were consistent with the results from TCGA.

To find out the precise mechanisms underpinning the observed HNRNPA2B1-dependent phenotypes, we conducted RIP assay using the HNRNPA2B1 antibody in OSCC cells to identify HNRNPA2B1 as m⁶A “reader” to combine with LINE-1 mRNA and promote the translation of LINE-1. We found that HNRNPA2B1 could enrich LINE-1 mRNA, implying that LINE-1 may be regulated in protein level upon interaction with HNRNPA2B1 ( Figure 7K ).

These results explain that HNRNPA2B1 may, as an m⁶A “reader,” regulate the proliferation, migration, and invasion of OSCC cells by targeting EMT via the LINE-1/TGF-β1/Sail signaling pathway.