Molecular Characteristics, Prognostic Value, and Immune Characteristics of m6A Regulators Identified in Head and Neck Squamous Cell Carcinoma.

N6-methyladenosine (m6A) plays crucial roles in a diverse range of physiological and pathological processes, and it is believed that it tremendously promotes neoplasia and progression. However, knowledge of the molecular characteristics of m6A modification, its prognostic value, and the infiltration of immune cell populations in head and neck squamous cell carcinoma (HNSCC) is still insufficient. Therefore, a pan-cancer genomic analysis was systematically performed here by examining m6A regulators at the molecular level within 33 multiple cancer types, and the correlations between the expression of m6A molecules were researched using datasets from The Cancer Genome Atlas (TCGA). Based on the above analysis, insulin-like growth factor 2 mRNA-binding protein 2 (IGF2BP2) is upregulated in HNSCC and may serve as an independent prognostic factor of overall survival, thus showing potential as a prognostic biomarker in HNSCC. Genetic alteration analyses elucidated the reasons for the abnormal upregulation of IGF2BP2 in HNSCC. As a result, IGF2BP2 was selected for further univariate and multivariate analyses. The functions of the related genes were annotated through gene set enrichment analysis, and the activation states of multiple biological pathways were shown by gene set variation analysis. We found that LRRC59 and STIP1 may act as IGF2BP2-associated genes to have a regulatory function in the m6A modification. In addition, we found that the status of immune cell infiltration was correlated with the level of IGF2BP2 gene expression. Our results provide supplementation at the molecular level for epigenetic regulation in HNSCC and insight into effective immunotherapy targets and strategies.

Introduction

The classic epigenetic modifications are DNA methylation, histone modification, and chromatin remodeling, which take part in many basic biological activities of carcinogenesis, including malignant development, and are associated with the prognosis of various kinds of cancer, including head and neck squamous cell carcinoma (HNSC, HNSCC) (1). In the past, reports on tumor methylation modifications have concentrated mostly on DNA methylation and offered potential biomarkers for cancer diagnosis (2). In addition to these traditional epigenetic modifications, mRNA modification is known as another epigenetic regulation of gene expression. Its dynamic regulation and disorder are closely related to translation control, RNA splicing defects, and the occurrence, maintenance, and progression of various tumors (3, 4). The degree of RNA modification is highly sensitive to changes in the cell microenvironment or the transformation of physiological states, and RNA modification changes will in turn influence cell regulation and adaptation (5). Among the types of RNA modifications, N6-methyladenosine (m6A) is recognized as the first and most common modification in eukaryotic mRNA (6).

m6A modification has been studied for 40 years. As a posttranscriptional modification with dynamic and reversible characteristics, m6A methylation has a profound impact on gene regulation (7), and its biological importance has also been emphasized by recent research progress. m6A modification participates in numerous basic biological processes, such as mRNA stability, mRNA translation, RNA splicing, and the phase separation potential of mRNA (8–10). Many physiological processes of the human body are closely related to m6A, such as spermatogenesis (11) and spermatogonial differentiation (12), hematopoietic stem cell development (13), and antiviral innate immunity (14). Moreover, a disturbance in m6A modification is associated with the pathogenesis of diverse cellular processes, including reduced cell proliferation, impaired self-renewal ability, developmental defects, and cell death (15). m6A modification has crucial functions in many human diseases, including infertility (16), nervous system diseases (17), early embryonic retardation (18), immune diseases, and multiple cancers (19). Currently, m6A modification, especially the role of its regulatory enzymes (writers, erasers, and readers), in cancer biology has been recognized as a prominent hot spot in terms of tumorigenesis, malignant progression, and potential biological target screening. Writers form a multisubunit methyltransferase complex to upregulate the level of m6A, while erasers are m6A demethylases, making this event reversible. Readers can decode m6A methylation information and convert it into effectors of functional signals. Recently, increasing evidence has shown that changes in m6A regulators can promote the development of several cancers (20), including glioblastoma, breast and cervical cancers, liver and gallbladder cancers, and bladder and prostate cancers (21–23). However, despite some progress in revealing their role in multiple cancers, the characteristics of m6A modification in HNSCC are still very inadequate.

HNSCC is a highly invasive malignant tumor that originates from the mucous membrane of the hypopharynx, larynx, oral cavity, and oropharynx, accounting for 90% of all head and neck cancers (24). HNSCC has more than 650000 confirmed cases worldwide, causing more than 300000 deaths every year (25) and ranking as the sixth most common and fatal cancer throughout the world (26). Although much progress has been made in the treatment methods (surgery, chemotherapy, and radiotherapy), the clinical prognosis of HNSCC patients is still poor. The 5-year survival rate of patients is still approximately 50%, which remains to be improved, and suggests that approximately 30-50% of HNSCC patients will experience local recurrence and distant metastasis (24). Therefore, the search for specific molecular markers of HNSCC is of great significance for better understanding of its progress and for developing new therapeutic methods. At present, an increasing number of studies have concentrated on the mechanism of pathogenesis, especially genetic and epigenetic events, at the molecular level (27). Currently, m6A modification has been reported to play a vital function in cell infiltration in the tumor microenvironment, and the m6A modification pattern can be used to help identify different immune phenotypes in gastric cancer (28). In addition, HNSCC was one of the first diseases for which immunotherapy was developed, which has greatly changed the therapeutic prospects of cancer patients. However, until now, there have been few reports about the relationship between m6A methyltransferase and clinicopathological features, its prognostic value, and the infiltration of immune cells in HNSCC. Accordingly, the function of m6A modification in tumor diagnosis, prognosis, and immune cell infiltration in HNSCC remains to be explored and could provide treatment targets for HNSCC.

As a result, we investigated the expression modes of 24 m6A regulatory genes in 33 cancers and explored the differential expression of IGF2BP2 in HNSCC and its clinical significance. Based on the clinical and sequencing data of The Cancer Genome Atlas (TCGA)-HNSCC cohort, the relationships among differential expression, genetic alterations, and clinicopathological features, including survival, were analyzed. Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) were used to analyze the molecular mechanism of IGF2BP2 abnormalities in HNSCC. On the other hand, the characteristics of m6A genes, including IGF2BP2, were analyzed to examine their relationship with immune cell infiltration in HNSCC. The purpose of this study was to provide supplementation at the molecular level for epigenetic regulation in HNSCC and ideas for effective immunotherapy targets and strategies.

Materials and Methods

Acquisition of Datasets and Data Preprocessing

For datasets in the TCGA, the R package TCGAbiolinks (29) was used to download RNA sequencing data and full clinical annotation information from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) and was specialized for an integrated analysis of GDC data (29). Patients for whom survival information was lacking were excluded from further evaluation. Information on somatic mutation data for HNSCC patients was obtained from the TCGA database as well. The expression of m6A regulatory genes in 33 cancers was analyzed with the TIMER2.0 database (http://timer.comp-genomics.org/) (30), in which the Wilcoxon test (also named the Mann–Whitney test) was applied to evaluate the P value. Sequencing and copy number variation (CNV) data provided by cBioPortal (http://www.cbioportal.org/) for cancer genomics were used to detect the genetic alterations in m6A regulatory proteins within 504 HNSCC patients/samples. IGF2BP2 mRNA expression, somatic mutations, CNVs, and the level of methylation in HNSCC were analyzed with UCSC Xena (https://xenabrowser.net/heatmap/). Twenty-four interactions between m6A regulators were elucidated with the help of the STRING database (http://www.string-db.org/). Immune cell infiltration was analyzed using the TIMER2.0 database.

Prognostic Value and Survival Analyses

Univariate and multivariate Cox regression analyses were used to assess the prognostic value of m6A regulators through their expression pattern in the TCGA. Both the “survival” and “survminer” (https://CRAN.R-project.org/package=survminer) R packages were used to investigate every gene. Four hundred and twenty-nine HNSCC patients were divided into different subgroups (high and low groups) by their characteristics and median gene expression values.

Differential Expression Analyses

The TCGAbiolinks R package was used to download TCGA-HNSC read count data. Patients with HNSCC were assigned to high and low IGF2BP2 groups according to the median cutoff value. A total of 11,297 long noncoding RNAs (lncRNAs) were ultimately screened for differential expression analyses. The DEseq2 R package (version 1.26.0) in R software (https://www.r-project.org/) was used to perform data standardization and calculate differences in the expression data. An absolute log2-fold change (FC) > 1 and p-value < 0.05 were considered the criteria for the differential expression of lncRNAs. The ggplot2 R package was used to draw a volcano plot.

Gene Set Cancer Analysis (GSCA) Lite Analysis

The differential expression profiles and survival outcomes of 24 m6A regulator genes among 33 cancer types were assessed using the GSCALite web server (31). Only 14 cancer types (THCA, KIRP, BLCA, LIHC, HNSC, BRCA, LUAD, PRAD, ESCA, KICH, LUSC, KIRC, STAD, and COAD) had paired samples (paired tumor-normal tissue) with available expression data. GSCALite was also used to study the effect (activation or inhibition) of m6A molecules on pathways related to cancer. In GSCALite, the scores of 7876 samples, 10 cancer-related pathways, and 32 cancer types (acute myeloid leukemia (AML) was excluded) were obtained from reversed-phase protein array (RPPA) data from The Cancer Proteome Atlas (TCPA). A total of ten cancer-related pathways, including apoptosis, PI3K/AKT, and TSC/mTOR, were identified, and m6A regulators that impacts (activate or inhibit) at least five different cancers were revealed by GSCALite.

Gene Set Enrichment Analysis (GSEA)

GSEA was performed with the JAVA program of the MSigDB v6.1 acquired from the website of the Broad Institute (http://software.broadinstitute.org/gsea/) (32). Then, we applied GSEA with standard settings as implemented in the R package clusterProfiler to assign pathway activity estimates to each sample. Finally, 18419 genes were included in the GSEA program, where “c2.cp.kegg.v7.1.symbols.gmt” was regarded as the hallmark gene set (33). A normalized p-value < 0.05 and a false discovery rate < 0.25 of the hallmark gene set were considered significantly enriched.

Gene Set Variation Analysis (GSVA) and Functional Annotation

GSVA was performed with the “GSVA” R package to determine the differences in physiological processes based on m6A modification patterns. GSVA is widely employed as a method to estimate variation in signaling pathways and physiological processes within a nonparametric and unsupervised method for samples of an expression dataset (34). Thus, the gene set “c2.cp.kegg.v7.1.symbols.gmt” was accessed through the MSigDB for GSVA. An adjusted p-value < 0.05 was considered statistically significant.

The Human Protein Atlas

The Human Protein Atlas (https://www.proteinatlas.org/) provides pathological specimens for the immunohistochemical analysis of prognostic RNAs, such as m6A methyltransferases. A survival analysis of associated genes in HNSCC was also performed. Information on the patients, staining intensity, and quantity can be obtained online.

Statistical Analysis

SPSS 20.0 (SPSS, Chicago, IL), R software V3.6.2 and GraphPad Prism 8 (GraphPad Software, Inc., La Jolla, CA, USA) were applied for statistical analysis. Continuous variables are presented as the mean ± standard deviation. Univariate and multivariate Cox regression analyses were used to illustrate their prognostic value and associations with clinical characteristics at the molecular/pathological level. Kaplan–Meier curves with a two-sided log-rank test were used to compare the overall survival (OS) of the patients in the high and low expression groups. p-values < 0.05 were considered statistically significant.

Results

Gene Expression Landscape of m6A Regulators

The genome-wide omics data, consisting of more than 10,000 human samples of 33 cancer types, were acquired from the TCGA for analyses ( ). The literature was reviewed, and a catalog of 24 genes that function mainly as regulators of m6A modifications was curated here: 8 writers (WTAP, METTL3, METTL14, ZC3H13, RBM15, RBM15B, KIAA1429 (VIRMA), and CBLL1), 2 erasers (ALKBH5, FTO) and 14 readers (HNRNPA2B1, HNRNPC, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, RBMX, FMR1, ELAVL1, and LRPPRC). The proportions of readers, writers, and erasers among m6A regulators are listed in . The correlations between the expression of 24 m6A regulators in the TCGA-HNSC cohort were determined by Pearson correlation analysis, and the results are displayed in . The TIMER database was used for analysis, and we found that m6A regulators were abnormally expressed in cancers in a cancer-specific pattern. Compared with their paired normal controls, RBM15, METTL3, and IGF2BP1/2/3 were obviously upregulated in the TCGA-HNSC dataset, while WTAP, METTL14, ZC3H13, YTHDF2, YTHDF3, ALKBH5, FTO, CBLL1, and LRPPRC were not upregulated ( ). To better visualize the interactions between the 24 m6A regulators, their relationships were determined with the STRING database ( ). The GSCALite web server was used, and we found twenty-four m6A regulators to be differentially expressed in 14 cancer types ( ). We found that IGF2BP1/2/3 were upregulated in HNSCC, consistent with the data from TIMER.

Figure 1

Gene expression landscape of m6A regulators in pan-cancer. (A) The proportions of readers, writers, and erasers among m6A regulators. (B) Correlations between the expression of 24 m6A regulators. The correlation between RBMX and HNRNPA2B1 is illustrated as a scatter plot (Pearson correlation). A positive correlation is indicated in red, and a negative correlation is indicated in blue. The color intensity and size of the circle are proportional to the correlation coefficient. (C) The expression profiles of the 24 m6A regulatory genes in the TCGA pan-cancer dataset originating from the TIMER database. (D) The interactions among 24 m6A regulatory genes (STRING). (E) GSCALite was used to analyze the gene expression profiles of 24 m6A regulatory genes in 14 types of cancer.

Relationships Between m6A Regulatory Genes and Pan-Cancer

Among the different types of cancer, the prevalent genetic and expression alterations in m6A regulators could largely give rise to the breakthrough of translational medicine, as m6A regulatory genes were closely related to patient OS across 33 cancer types in our study. It was also implied by our results that 24 m6A regulatory genes were related to the OS of patients with at least one cancer type, where some of the genes, including IGF2BP1/2/3, showed oncogenic features, and their increased expression was related to worse survival in patients with various cancers. In particular, upregulation of the IGF2BP2 gene was associated with poor patient survival in seven cancer types, including HNSCC ( ).

To further explore the underlying mechanisms by which m6A regulatory genes participate in multiple cancers at the molecular level, the effects (activation or inhibition) of m6A regulatory genes on pathways accompanying cancers were analyzed with GSCALite. It was found that m6A regulatory genes are functionally expressed and are related to oncogenic activation or inhibitory pathways ( ). A large number of activated pathways, namely, the cell cycle, epithelial-mesenchymal transition (EMT), DNA damage response, PI3K-AKT, RTK, and RAS/MAPK, were found to be associated with the expression of IGF2BP2, RBMX, and HNRNPC. These results confirmed that m6A regulators are functionally related to multiple cancers.

Single Nucleotide Polymorphism (SNP) Mutations in m6A Regulatory Genes

We found that single nucleotide polymorphism (SNP) occupied the first place in variant type ( ). SNP analysis showed that missense mutations were the most frequent variants in the m6A-regulatory genes ( ). The most frequent SNV type was C > T ( ).

In addition, we summarized the incidence of single nucleotide variants (SNVs) in 24 m6A regulators in 33 pan-cancer samples with GSCALite. Among 1443 samples, 1002 harbored mutations in m6A regulators, with a frequency of 69.44%, among which ZC3H13 showed the highest mutation frequency, followed by KIAA1429 ( ).

The Prognostic Role of m6A Regulators and Their Correlations With Clinical Characteristics in HNSCC

Considering the relationship between m6A methyltransferases and the malignant progression of HNSCC, univariate and multivariate survival analyses were performed to further determine the prognostic value of m6A methyltransferases in HNSCC using Cox proportional hazards models of their expression levels in the TCGA dataset. A total of 429 HNSCC patients were assigned to two groups (high and low) based on the median expression value of each gene. Stage IV disease and high LRPPRC and IGF2BP2 expression levels were shown to be independent prognostic factors for OS. Multivariate analysis showed that patients aged ≥60 years, with stage IV disease, of African American race, and high ELVAL1 and IGF2BP2 expression were risk factors for poor OS (hazard ratio (HR)> 1) ( ). As shown in , IGF2BP2 was correlated with T-stage based on the median cutoff of IGF2BP2 (P = 0.023). However, IGF2BP2 did not correlate with clinical characteristics such as age, race, stage, or sex. ELVAL1 was correlated with the clinical characteristic race (P = 0.014). LRPPRC was correlated with the clinical characteristics sex (P = 0.016), T-stage (P = 0.044), and race (P = 0.040) ( ).

Figure 2

Univariate Cox regression analysis of m6A regulators in HNSCC patients in the TCGA dataset. Tumor node metastasis classification: T (tumor), N (lymph node), and M (metastasis).

The Expression and Prognostic Value of IGF2BP2

Based on the above analysis, we next selected IGF2BP2 for the multivariate analysis of 24 m6A regulatory genes. The expression of the IGF2BP2 gene in pan-cancer was examined with TIMER ( ), and we verified the expression and prognostic value of the IGF2BP2 gene in HNSCC using GEPIA2 ( ). Kaplan–Meier curve analysis also indicated that patients with higher IGF2BP2 expression experienced significantly worse OS than their counterparts ( ).

Figure 3

The expression and prognostic value of IGF2BP2 in HNSCC. (A) The expression of the IGF2BP2 gene in pan-cancer was analyzed with TIMER. (*p < 0.05, **p < 0.01, ***p < 0.001). (B, C) The expression and prognostic value of the IGF2BP2 gene in HNSCC was analyzed with GEPIA2. (HR, hazard ratio) (D) Kaplan–Meier survival curves illustrate the difference in OS based on the expression level of IGF2BP2.

Genetic Alterations in m6A Regulatory Genes and Mutation, CNV, and Methylation Analyses of IGF2BP2

The incidence of genetic alterations in the 24 m6A regulators in HNSCC (TCGA, Firehose Legacy) was analyzed by using the cBioPortal tool. Twenty-four m6A methylation regulators showed genetic alterations in 227 (45%) queried patients/samples (504 patients/samples). The most common genetic alterations were amplification and deep deletion mutations. In our analysis, IGF2BP2 exhibited the highest mutation frequency (an alteration rate of 20%), followed by KIAA1429 (VIRMA) (an alteration rate of 8%), while both METTL14 and YTHDF2 showed few mutations in HNSCC samples ( ). Furthermore, the potential mechanism of dysregulated IGF2BP2 was studied because it was identified as an independent prognostic factor for OS. In the UCSC Xena database, 604 samples from the TCGA-HNSC dataset were analyzed, which revealed that IGF2BP2 expression was associated with CNV and partial DNA methylation sites but not with somatic mutations ( ). In conclusion, CNV and DNA methylation may give rise to the abnormal upregulation of IGF2BP2 in HNSCC.

Prediction of IGF2BP2 Transcription Factors

To identify the mechanisms controlling the tumorigenesis and poor prognosis of HNSCC at the transcriptional level, PROMO and Cistrome Data Browser were used to predict the transcription factor binding sites for IGF2BP2 ( ). The transcription factor CCAAT/enhancer-binding protein beta (CEBPB; C/EBP-β) was obtained after taking the intersections. IGF2BP2 (chr3:185643129-185821141) was found to be a putative target (score=2.843) of the transcription factor CEBPB by searching the ChIPseq results (GSM1010802) in Cistrome DB: 46089 (35) ( ). The same result was validated in the hTFtarget database (http://bioinfo.life.hust.edu.cn/hTFtarget/) (36). The motif of the transcription factor CEBPB is shown in (ID: MC00121, hits: 3187, cutoff: 6.012, z score: -107.2, -10log(pval): 690.776). A statistical plot showed the top 20 factors that regulate IGF2BP2 queried from Cistrome ( ), and CEBPB ranked first. In summary, the transcription factor CEBPB is most likely the transcription factor that regulates IGF2BP2.

Figure 4

Transcription factor prediction of IGF2BP2. (A) Visual presentation of the high binding peak of the transcription factor CEBPB in the promoter region of IGF2BP2 (B) Representative sequence logo of the binding specificity of the transcription factor CEBPB queried from Cistrome. (C) Static plot of the top 20 factors that regulate IGF2BP2. Y axis shows the RP score. X axis shows different factors. Dots on the x axis line indicate the same factor.

IGF2BP2-Associated Genes Mediate m6A Modification

Initially, we first collated an “A” list of dysregulated genes in HNSCC (n=5481) (log2FC cutoff = 0.58, p-value cutoff = 0.01), a “B” list of the top 1000 similar gene correlations of IGF2BP2 (n=1000), and then a “C” list of genes related to the survival of HNSCC patients (n=500). The number of genes in each list is displayed as a Venn diagram ( ). Twelve genes were common to all three gene sets. Therefore, they may act as IGF2BP2-associated genes and play a regulatory role in m6A modification. We analyzed their relationships with IGF2BP2 again. For example, we found that leucine-rich repeat-containing protein 59 (LRRC59) and stress-induced phosphoprotein 1 (STIP1) were correlated with IGF2BP2 (R=0.53 and R=0.53, respectively) ( ). LRRC59 and STIP1 were upregulated in HNSCC according to GEPIA2 (log2FC cutoff = 1, p-value cutoff = 0.01) ( ). High LRRC59 and STIP1 expression was correlated with poor OS in HNSCC patients ( ). It was predicted that IGF2BP2 facilitates the recruitment of RBPs to m6A-modified target mRNAs. By analyzing m6A Var, we analyzed the RBPs of potential associated genes, and their relationships were visualized with Cytoscape (37) ( ).

Figure 5

IGF2BP2-associated genes that mediate m6A modification. (A) Venn chart of IGF2BP2-related genes in HNSCC. (B) The correlation of LRRC59, STIP1, and IGF2BP2. (C) The expression of LRRC59 and STIP1 was confirmed with GEPIA2. T, tumor tissue. N, normal tissue. * means p-value < 0.01. (D) Kaplan-Meier survival curves illustrate the difference in OS according to the expression of LRRC59 and STIP1. (E) The relationships of RBPs and IGF2BP2-associated genes with the most potential.

Potential Regulatory Mechanisms of IGF2BP2

Based on the above analysis, a total of 301 upregulated and 149 downregulated lncRNAs are significantly related to IGF2BP2 expression ( ). Among them, some lncRNAs are associated with HNSCC. The lncRNA SOX2-OT was reported to inhibit PTEN expression to facilitate laryngeal squamous cell carcinoma (LSCC) development through EZH2-mediated H3K27me3 (38). The potential target miRNAs (at least four of seven platforms predicted) and corresponding lncRNAs were predicted with starBase v2.0 (http://starbase.sysu.edu.cn/) (39). Finally, a Sankey plot consisting of lncRNAs, microRNAs (miRNAs), and two mRNAs (LRRC59 and STIP1) was constructed to illustrate a competing endogenous RNA (ceRNA) regulatory network ( ).

Figure 6

Potential regulatory mechanisms of IGF2BP2. (A) The dysregulated lncRNAs associated with IGF2BP2 expression. (B) Sankey plot illustrating the ceRNA regulatory network.

Unsupervised Consensus Analysis of m6A Regulatory Genes

With the help of GSEA, the malignant hallmarks of tumors were investigated. The results showed that the tumor hallmarks primary immunodeficiency (NES=-2.21, normalized P < 0.0001), intestinal immune network for IgA production (NES=-1.95, normalized P < 0.0001), PPAR signaling pathway (NES=-1.80, normalized P< 0.001), basal cell carcinoma (NES=1.53, normalized P =0.013), Wnt signaling pathway (NES=1. 30, normalized P < 0.001), and Hedgehog signaling pathway (NES=1. 26, normalized P =0.080) were enriched in the IGF2BP2 group ( ).

Figure 7

Unsupervised consensus analysis of m6A regulatory genes. (A) GSEA revealed that the IGF2BP2 group was enriched in the hallmarks malignant progression and tumor immunology. (B) GSVA revealed the state of biological pathways in the high and low IGF2BP2 groups. The column represents the TCGA-HNSCC samples.

Furthermore, GSVA was applied to illustrate the biological behaviors among distinct m6A modification patterns. As shown in , the high IGF2BP2 group was markedly enriched in cell adhesion molecules, neuroactive ligand receptor interaction pathways, and calcium signaling pathways. In our study, it was concluded that low IGF2BP2 expression was enriched in the cytokine receptor interaction pathway.

Based on the above analyses, we speculate that the IGF2BP2-mediated m6A methylation modification has significant associations with the malignant progression and tumor immunology of HNSCC.

The Protein Expression of m6A Target Genes

To evaluate the expression differences in target genes at the protein level, immunohistochemistry data were obtained from pathological specimens from The Human Protein Atlas. We found that LRRC59 and STIP1 were dysregulated between HNSCC tumor tissues and normal tissues ( ).

Figure 8

The protein expression of m6A target genes. A Comparison of the protein expression of LRRC59 and STIP1 in HNSCC tumor tissues and normal tissues.

Relationship Between the Tumor-Immune Microenvironment and mRNA Expression of IGF2BP2

After determining the malignant hallmarks of IGF2BP2, the relationship between IGF2BP2 expression and tumor immune cell infiltration in cancers was explored. First, the correlation of IGF2BP2 expression and immune infiltration levels was evaluated in diverse cancer types. The relationship between IGF2BP2 and 21 kinds of immune cells indicated a negative correlation between IGF2BP2 and B cells, CD4+ T cells, follicular helper T cells, regulatory T cells (Tregs), myeloid-derived suppressor cells, and cancer-associated fibroblast cells in HNSCC ( ). For example, the expression level of IGF2BP2 had a negative correlation with the immune infiltration level of B cells (Spearman r = -0.624, p=6.39e-11). The mutation frequency of the IGF2BP2 gene in each TCGA cancer type is shown as a bar plot ( ). We found that HNSCC has the lowest IGF2BP2 mutation frequency among multiple cancer types. Therefore, IGF2BP2 gene mutations may have a weak relationship with immune cell infiltration across 22 multiple cancer types and immune cell types. The relative proportions of IGF2BP2 in different sCNA states in all tumor types are shown as a stacked bar chart ( ). We found that high amplification, arm-level gain, and diploidy were the most common states. It was also found that the distribution of immune cell infiltration in different sCNA states was different. For example, the differential T cell CD8+ (TIMER) and monocyte (XCELL) infiltration levels were observed in these sCNA states in HNSCC ( ).

Figure 9

Relationship between the tumor-immune microenvironment and IGF2BP2 expression. (A) The relationship between IGF2BP2 and immune cell infiltration in HNSCC. (B) Bar plot of IGF2BP2 mutation frequency in each TCGA cancer type. (C) The stacked bar chart shows the relative proportions of IGF2BP2 in different sCNA states in all tumor types. (D) Differential T cell CD8+ infiltration level in HNSCC (TIMER). (E) Differential monocyte infiltration level in HNSCC (XCELL).

Discussion

The abnormal epigenetic modification of RNA is important to tumorigenesis and cancer patient prognosis. As the most common modification type in eukaryotic mRNA, m6A is characterized by the N6 methylation of adenosine. Its role in cancer has gradually been recognized (6). It was recently found that m6A methylation has an indispensable function in innate immunity and the tumor microenvironment. This study focused on the abnormal gene expression of m6A regulators in pan-cancer and the molecular characteristics, prognostic value, and immune characteristics of m6A regulators in HNSCC.

In this study, we generated a catalog of 24 genes that act as regulators of m6A (8 writers, 2 erasers and 14 readers). After exploring their expression landscape in pan-cancer, it was noted that most m6A regulatory genes were abnormally expressed in tumors in a cancer-specific manner, and all the m6A regulatory genes were correlated with patient OS in at least one cancer type. These results were consistent with the findings of previous studies. For instance, METTL3 was found to be upregulated in colorectal cancer (CRC), and high METTL3 expression was associated with poor OS in CRC patients (40). Similarly, IGF2BP2 was upregulated and found to be related to a poor prognosis in pancreatic cancer (41) and esophagogastric junction adenocarcinoma (42). Moreover, it was found that combined levels of METTL3 and YTHDF1 could reflect the malignancy state and prognosis of hepatocellular carcinoma (HCC) (43). Our findings showed that most m6A regulatory genes were upregulated or downregulated in HNSCC, suggesting that these genes might be related to the carcinogenic activity of tumor cells and/or the prognosis of HNSCC patients. Among these m6A regulatory genes, IGF2BP2 was upregulated in HNSCC tumor tissues, and high IGF2BP2 expression was associated with poor patient OS in seven cancer types, including HNSCC. Similar studies in other cancers, such as colorectal, pancreatic, and breast cancers, have reported similar results (44–46). We found that the expression of IGF2BP2, RBMX, and HNRNPC was correlated with various activated cancer-related pathways, providing insight into the molecular mechanisms related to tumor progression. Furthermore, m6A regulators are useful as biomarkers for prognostic stratification (47). Our results showed that IGF2BP2 expression was negatively correlated with the prognosis of HNSCC, so IGF2BP2 is very likely an oncogene in HNSCC. In addition, CNV and DNA methylation, but not somatic mutations, give rise to the abnormal upregulation of IGF2BP2 in HNSCC based on data from TCGA-HNSCC patients. This result further shows that epigenetic modifications are also pivotal to the aberrant expression of IGF2BP2, an epigenetic modifier (41). This research suggests that the abnormal expression of m6A regulatory genes such as IGF2BP2 is common in HNSCC, but relatively little information on the function of IGF2BP2 in tumorigenesis is currently available, and its precise expression and mechanisms require further investigation.

Regarding the mechanisms controlling the tumorigenesis and poor prognosis of HNSCC at the transcriptional level, we found that the transcription factor CEBPB was the most likely to bind to the promoter region of IGF2BP2. Dysregulated IGF2BP2 could also affect the associated downstream genes via the m6A modification to accelerate cancer progression (41). Our study showed that 12 genes, including LRRC59 and STIP1, may act as IGF2BP2-associated genes to play a regulatory role in the m6A modification. LRRC59 and STIP1 have been reported to be closely associated with the progression of human malignancies. For instance, high LRRC59 expression is closely related to poor survival in LUAD patients and an aggressive phenotype of breast cancer (48). STIP1 is also upregulated in LUAD and may be associated with the enhanced proliferation, adhesion, and migration and inhibition of apoptosis in LUAD cells via the JAK2/STAT3 signaling pathway and EMT (49). Similar findings have been reported in HCC and ovarian, pancreatic, renal, and gastric cancers (50, 51). Noncoding RNAs are associated with the regulation of the m6A modification and can further affect the malignant progression of cancer (52, 53). Interestingly, a series of abnormally expressed m6A-modified miRNAs or lncRNAs could perform their function by acting as ceRNAs (54, 55), suggesting potential mechanisms in the dysregulation of IGF2BP2-associated genes in HNSCC. For example, IGF2BP2 serves as a reader of m6A-modified lncRNA-DANCR and can stabilize DANCR, which reversely contributes to tumor stemness-like properties and the pathogenesis of pancreatic cancer (41). So, in order to predict the key factors that exist in this regulatory relationship, we derived a regulatory network of IGF2BP2-associated genes, miRNAs, and lncRNAs.

In recent years, there has been increasing interest in m6A modification in multiple tumor biological processes and signaling pathways not limited to the growth of tumor stem cells, neoplasia (56–58), RNA metabolism (59), and the DNA damage response after radiotherapy and chemotherapy (60). By performing GSEA, the malignant hallmarks of IGF2BP2 in HNSCC were analyzed and may include basal cell carcinoma, the Wnt signaling pathway, the Hedgehog signaling pathway, and some signaling pathways related to tumor immunology, such as primary immunodeficiency and the intestinal immune network for IgA production. GSVA showed that the high IGF2BP2 group was markedly enriched in the cell adhesion molecule CAMS, neuroactive ligand receptor interaction pathways, and the calcium signaling pathway. This study uncovered that IGF2BP2-mediated m6A methylation was strongly correlated with the malignant progression of HNSCC and tumor immunology signaling pathways. The subsequent assays revealed a negative correlation between IGF2BP2 and multiple tumor-infiltrating immune cells, such as B cells, CD4+ T cells, follicular helper T cells, and Tregs. Our findings may provide novel insight into immunotherapy targets in HNSCC.

However, this study also has some limitations. For example, the potential mechanisms by which IGF2BP2 is upregulated are only descriptive. In the future, more basic functional mechanistic studies will be needed to validate our findings, and the clinical application of m6A methyltransferase in HNSCC needs further exploration.

Conclusions

In conclusion, by systematically analyzing the gene expression landscape of 24 m6A regulators in pan-cancer, IGF2BP2 was found to be highly expressed in HNSCC, is correlated with tumor progression and a low survival rate, and is involved in tumor immune invasion. More importantly, IGF2BP2 might be a promising candidate as a molecular biomarker for monitoring HNSCC development and provide ideal guidance for selecting therapeutic strategies. This study provides a supplement at the molecular level for epigenetic regulation in HNSCC and ideas for effective immunotherapy targets and strategies.

Data Availability Statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding authors.

Author Contributions

XCG and YHZ were the major contributors to the writing and revision of the manuscript. ZMZ and ZRZ collected the related references and participated in discussions. XCG and QL made substantial contributions to the conception or design of the work. QL, WTY, and HW approved the final version of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was financially supported by grants from the Natural Science Foundation of Hebei Province (no. H2020201050), the National Innovation and Entrepreneurship Training Plan in 2020 (no. 202014432031), Science and Technology Capacity Improvement Projects of Hebei University of Chinese Medicine in 2019 (no. KTZ2019019), Outstanding Student Scientific Research Ability Improvement Projects of Hebei University of Chinese Medicine in 2019 (no. YXZ2019002), and Graduate Innovative Ability Training Projects of Hebei Province in 2020 (no. XCXZZBS2020002) (no. hbu2020bs003).

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.