Identification of m6A-Regulated Ferroptosis Biomarkers for Prognosis in Laryngeal Cancer

doi:10.21203/rs.3.rs-3890785/v1

Download PDF

Research Article

Identification of m6A-Regulated Ferroptosis Biomarkers for Prognosis in Laryngeal Cancer

https://doi.org/10.21203/rs.3.rs-3890785/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

N6-methyladenosine (m6A) RNA methylation and ferroptosis are associated with laryngeal cancer (LC) development. Accordingly, further research on related molecular mechanisms and pathology of LC is necessary.

Methods

Weighted gene co-expression network analysis (WGCNA) and correlation analysis were used to identify differentially expressed m6A-related ferroptosis genes (DE-MRFG) in LC. Univariate Cox and least absolute shrinkage and selection operator (LASSO) regression were utilized for feature selection and risk model construction. Then, a nomogram was built based on the independent prognostic factor identified using univariate and multivariate Cox regression. Mutation analysis, immune-related analysis, and drug sensitivity prediction were applied to analyze the utility of the risk model as much as possible. Additionally, qRT-PCR and western blot were performed to detect the TFRC, RGS4, and FTH1 expression.

Results

We identified 83 genes as DE-MRFG in LC. Three model genes (TFRC, RGS4, and FTH1) were identified to build a risk model using the univariate Cox and LASSO regression algorithms. Receiver operating characteristic (ROC) analysis verified the accuracy of the risk model. Furthermore, calibration curves and ROC analysis indicated the good performance of the nomogram in predicting overall survival (OS). Moreover, the mutation analysis indicated that multiple genes were mutated in the high- and low-risk groups. Based on the analysis of the immune reaction in LC, immune checkpoint PD-L1 was significantly related to the risk score and was up-regulated in the high-risk group. Tumor Immune Dysfunction and Exclusion (TIDE) and Estimation of STromal and Immune cells in MAlignant Tumors using the Expression data (ESTIMATE) algorithm showed a positive relationship between risk score and TIDE or ESTIMATE score. Furthermore, drug sensitivity prediction found that 19 chemotherapy drugs were strongly correlated with a risk score. TFRC, RGS4, and FTH1 exhibited high expression levels in 30 laryngeal carcinoma tissues and cell lines (TU212, TU686, and AMC-HN-8). Notably, TFRC and FTH1 expression levels were significantly associated with patient prognosis.

Conclusion

Three prognostic genes, TFRC, RGS4, and FTH1, were identified as m6A-regulated ferroptosis biomarkers in LC, providing insights into LC treatment and prognosis.

laryngeal cancer

m6A

ferroptosis

biomarker

prognosis

Laryngeal cancer (LC), a prevalent malignant tumor of the head and neck region, is the third most common malignant tumor in this area. The primary symptoms of LC are dyspnea, cough, hoarseness, and hemoptysis (1). The progression of advanced LC significantly impacts the quality of life and survival of affected individuals (2). Surgery is a primary therapeutic approach for LC, effectively eliminating lesion tissue and enhancing the short-term prognosis of patients (3, 4). Nevertheless, it harbors numerous nerve tissues and lymph nodes due to the distinctive anatomical positioning of the larynx, rendering it susceptible to lymph node metastasis and resulting in suboptimal long-term prognostic outcomes (5). Advancements in molecular biology have led to the emergence of immunotherapy as a prominent area of focus in treating malignant tumors. This therapeutic approach can enhance tumor drug sensitivity, augment the effectiveness of anti-tumor medications, and mitigate adverse effects. Consequently, investigating novel biomarkers and their underlying mechanisms holds promise for informing clinical treatment strategies (6, 7).

N6-methyladenosine (m6A) is the most prevalent RNA modification, initially discovered in mRNA and subsequently identified in precursor mRNAs, circular RNA (circRNAs), and long non-coding RNAs (lncRNAs) (8, 9). m6A is crucial in various biological processes, including regulating mRNA stability, pre-mRNA splicing, translation, and DNA damage repair (10, 11). Our previous investigation demonstrated that RBM15-mediated m6A modification of TMBIM6 mRNA enhances TMBIM6 stability via IGF2BP3-dependent mechanism. Furthermore, RBM15 promotes the malignant progression of LC by upregulating TMBIM6 (12). Iron is a crucial trace element necessary for the oxidative respiration of mitochondria within biological cells. Research has demonstrated that elevated intracellular iron levels can trigger a distinct type of cellular demise known as ferroptosis, distinguishing it from other cell death mechanisms, such as autophagy and cell death (13, 14). The key features of this phenomenon encompass the buildup of intracellular iron, heightened production of reactive oxygen species (ROS), and lipid peroxidation (15). Ferroptosis impedes the tumor cell growth and enhances the responsiveness to chemotherapeutic agents (16, 17). Ferroptosis is significantly associated with head and neck squamous cell carcinoma. Suppressing glutaredoxin 5 gene expression induces iron death in cisplatin-resistant cells (18, 19). A substantial research foundation exists to explore the relationship between m6A, ferroptosis, and LC progression and prognosis.

This study employed bioinformatics to identify m6A-related ferroptosis biomarkers in LC using publicly available databases. We developed a risk model based on ferroptosis genes regulated by m6A. The univariate Cox regression analysis results revealed that TFRC, RGS4, and FTH1 genes were significantly associated with the survival status of patients with LC. Additionally, we assessed the risk scores of the clinical characteristics of different LC samples. The risk model demonstrated a predictive relationship between the risk score and the efficacy of immune therapy and chemotherapy drugs. Besides, TFRC, RGS4, and FTH1 expression levels were examined in three laryngeal carcinoma cell lines (TU212, TU686, and AMC-HN-8) and 30 pairs of laryngeal carcinoma clinical specimens. Our results indicated that TFRC, RGS4, and FTH1 expression levels were significantly higher than normal human bronchial epithelial cells (NHBEC). Furthermore, gene expression levels were significantly higher in laryngeal carcinoma tissues than in the adjacent tissues. These findings hold substantial importance for further elucidating the underlying mechanisms of LC.

2.1 Data source

Two LC datasets, TCGA-HNSC and GSE65858, were acquired from the Cancer Genome Atlas (TCGA) (http://cancergenome.nih.gov/) and the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). The GSE65858 dataset contains gene data sequenced using the GPL10558 Illumina HumanHT-12 V4.0 expression bead chip platform (20). Forty-eight tumor samples from GSE65858 with survival information were enrolled in this study, which was used to validate the risk model. A total of 110 tumor samples and 12 control samples from the larynx were obtained from the TCGA-HNSC dataset (21), which was used as the major source for all analyses. Twenty-three m6A RNA methylation-related genes were identified based on a previous study (22). A total of 567 genes related to ferroptosis were sourced from the FerrDb database (http://www.zhounan.org/ferrdb/current/).

2.2 Function analysis of differentially expressed m6A-related ferroptosis genes (DE-MRFG)

First, the ‘DESeq2’ R package (version 1.34.0) (23) was applied to identify differentially expressed genes (DEGs) in LC patients and healthy individuals in DE-MRFG based on the data of the TCGA-HNSC dataset. The filter criteria were adjusted p-value < 0.05 and |log₂FC| > 0.5. DEGs are shown in the volcano map and heatmap. Second, LC-related genes were identified using weighted-gene co-expression network analysis (WGCNA) (24). WGCNA constructed and screened gene modules were strongly correlated to LC using the R package ‘WGCNA’ (version 1.70-3). Third, LC-related genes intersected with DEGs and ferroptosis-related genes, obtaining differentially expressed ferroptosis-related genes (DE-FG) in LC. Finally, DE-MRFG in LC was determined using Pearson’s correlation analysis between m6A-related genes and DE-FG. A p-value < 0.05 and |R| > 0.3 were the filter criteria. Gene ontology (GO) and the Kyoto encyclopedia of genes and genomes (KEGG) analysis were applied using the ‘clusterProfiler’ package (version 4.2.2) to probe the molecular mechanism of DE-MRFG (25).

2.3 Establishment and verification of a risk model

Initially, univariate Cox regression was applied using the ‘survival’ R package (version 3.2–13) to obtain biomarkers, identify prognostic genes, and build a risk model (Therneau, 2012). A survival analysis R package (2.37-2) based on the expression data of DE-MRFG in LC samples was used to obtain data from the TCGA-HNSC dataset (n = 108, two cancer samples were filtered, and patient survival was less than 30 days). A p-value < 0.05 and HR ≠ 1 were considered as the filter criteria. Then, a risk model was established using the least absolute shrinkage and selection operator (LASSO) regression and the ‘glmnet’ R package (version 4.1-2) (26). The LASSO regression formula was calculated as follows: Risk score =\(.\) Data from the LC samples from TCGA and GSE65858 were used to assess and verify the risk model. First, all LC samples were classified into two groups for use in the following analyses, according to the risk score optimal threshold calculated by the risk model. One group contained LC samples with high-risk scores, and the other contained low scores. Second, the survival state of different groups was analyzed using Kaplan-Meier (KM) survival analysis and the ‘survminer’ package (version 3.3.5) (Kassambara, 2017). Survival curves were drawn using 'ggplot2' [R package survminer version 0.2.0]). Third, the accuracy of survival prediction was assessed using receiver operating characteristic (ROC) analysis and the ‘survivalROC’ package (version 1.0.3) (Heagerty, 2013). SurvivalROC: time-dependent roc curve estimation from censored survival data). Fourth, survival state distribution maps were drawn to investigate the survival time of patients with different risk scores. Finally, a heatmap of the biomarker expression was drawn. In summary, survival, ROC, and expression analyses were used in both TCGA-HNSC and GSE65858 to assess the performance of the risk model.

2.4 Enrichment analysis of DEGs in different groups

DEG analyses were performed using the data of the TCGA-HNSC dataset and the ‘DESeq2’ package to research the potential influence of DEGs in high- and low-risk groups. DEGs from the different groups are shown in the volcano map. Genes with an adjusted p-value ≥ 0.05 and |log₂FC| ≤ 0.5 were filtered. GO and KEGG analyses were performed to investigate the function of the DEGs. The top five GO items in the three categories of GO analysis and the top 10 KEGG pathways were drawn into bubble plots.

2.5 Relationship between risk model and clinical characteristics

Clinical information and samples with survival data from the TCGA-HNSC dataset were obtained to explore the influence of different clinical characteristics in the risk model. First, Wilcoxon (for two groups) and Kruskal-Wallis (for more than two groups) tests were used to compare the risk scores of different clinical features, including gender, age, pathological N stage, pathological T stage, and tumor Stage. Violin plots were used to reveal the distribution of risk scores in different clinical characteristics. Then, a heatmap was used to show the expression of biomarkers with different clinical characteristics. Finally, KM survival analyses of the high- and low-risk groups were used to investigate the different survival states according to different clinical features.

2.6 A nomogram model of LC

First, univariate Cox regression analysis was applied for risk score, age, gender, tumor stage, Pathologic_T, and Pathologic_N of LC samples from the TCGA-HNSC dataset. Forest plot demonstrated the result of these factors using ‘forestplot’ (version 2.0.1) (Max and Thomas (2022). Forestplot: Advanced Forest Plot Using 'grid' Graphics, https://rdrr.io/cran/forestplot/). The factors with p-value < 0.05 were used for multivariate Cox regression analysis to identify the independent prognosis factor. Then, a nomogram was built to predict the two-, three-, and four-year overall survival (OS) of LC patients based on the independent prognostic factors and verified using calibration curves and ROC analysis.

2.7 The relationship of risk model and immunotherapy

The single-sample gene set enrichment analysis (ssGSEA) was utilized to screen out the differences in immune cells between high- and low-risk groups. During this process, expression data from TCGA-HNSC was used. Statistical differences were calculated using Wilcoxon and corrected using Benjamini and Hochberg. Then, correlation analysis of immune infiltration and biomarkers was used. The relationship analysis results were visualized using the ‘ggpubr’ package (version 0.4.0) (Aut, 2017). ggpubr: 'ggplot2' Based Publication Ready Plots). Then, the immune infiltration situation was evaluated using the Estimation of STromal and Immune cells in MAlignant Tumors using the Expression data (ESTIMATE) algorithm and the ‘estimate’ package (version 1.0.31) in LC samples (27). The ESTIMATE results were drawn into violin plots, and the statistical significance was tested using the Wilcoxon test.

The ‘maftools’ package (version 2.10.5) (28) was used to obtain information on mutation in high- and low-risk groups from the TCGA-HNSC dataset and visualize the top 20 gene-mutation results using waterfall plots. Immune checkpoints (PD-1, IFN-γ, CTLA4, PD-L1, and IL-2) were analyzed using the TCGA-HNSC dataset to probe immune therapy response in the groups classified by risk score. Initially, the relationship between immune checkpoints and risk scores was studied using Spearman correlation analysis. Then, expression and KM survival analysis of the immune checkpoint were applied. Afterward, the tumor immune dysfunction and exclusion (TIDE) algorithm was used to estimate the tumor's immune escape capability. The TIDE, dysfunction, and exclusion score of each LC sample were calculated. The correlation between them and the risk score was analyzed using Spearman correlation analysis. A higher TIDE score means that the tumor's immune escape ability is higher, and the response to immune therapy is poor in patients with LC.

2.8 Prediction of Chemotherapy Drug Sensitivity

The TCGA-HNSC dataset was used to assess the response of chemotherapy drugs in LC and investigate the relationship between drug response and risk score. A total of 150 chemotherapeutics of LC were sourced from the Genomics of Drug Sensitivity in Cancer (GDSC) database. Then, the half maximal inhibitory concentration (IC₅₀) value of each drug in different LC samples was calculated using the ‘pRRophetic’ package (version 0.5) (29). IC₅₀ reflects the drug response; the smaller the IC₅₀, the stronger the drug's effect. Twenty out of 150 chemotherapeutics of LC were contained in ‘pRRophetic’ package, which was used to conduct difference and correlation analyses between the chemotherapeutics and risk score.

2.9 Tissue samples and Cell culture

A total of 30 pairs of laryngectomy and adjacent non-tumor tissues were collected by the Department of Otorhinolaryngology at the Shaanxi Provincial People's Hospital and the Second Affiliated Hospital of Harbin Medical University. Patients with LC diagnosed before surgery and those without any previous treatment from 2018 to 2019 were included in this study. This study was endorsed by the Harbin Medical University and Shaanxi Provincial People's Hospital Ethics Committee endorsed.

2.10 qRT-PCR

Total RNA was extracted from LC tissues and cell lines using Trizol reagent (Invitrogen, Carlsbad, CA, USA) per the manufacturer's instructions. qRT-PCR was conducted as previously described (30). The 2^−ΔΔCt method was used to calculate the data. The primer sequences are listed in Table S1.

2.11 Western blot

The cells were harvested and lysed on ice for 30 min. Lysate and cells were collected and centrifuged at 12,000 rpm at 4°C. The total protein concentration in the supernatant was determined using the BCA method. The remaining sample was added to the SDS loading buffer and heated at 100 ℃ for 5 min. The protein samples were electrophoresed and transferred to the PVDF membrane at 300 mA. The PVDF membrane was blocked with 5% milk for 1 h. The PVDF membrane was incubated with a primary antibody (Abcam, Shanghai, China) at 4 ℃ overnight. Then, it was incubated with a 1:10,000 dilution of HRP-anti-mouse secondary antibody and developed with the ECL substrate.

2.12 Statistical analysis

The R programming language was applied for all the analyses. If it was not mentioned, the p-value was less than 0.05 and considered statistically significant. Paired samples were conducted using a two-tailed paired sample t-test.

3.1 DE-MRFG in LC was associated with oxidation-reduction and cancer

We screened 14,742 DEGs (Table S2), including 431 up- and 1,065 down-regulated genes, between LC and control samples in the TCGA-HNSC dataset (Fig. 1A). The top 10 DEGs for up- and down-regulated gene expressions were drawn into the heatmap (Fig. 1B). Then, we obtained 5,832 LC-related genes using WGCNA. Combined with 5,832 LC-related genes, 14,742 DEGs, and 567 ferroptosis-related genes, an intersection gene, 97 DE-FG in LC, was obtained (Fig. 1C, Table S3). Among them, 83 genes are related to m6A through correlation analysis between DE-FG and m6A-related genes, named DE-MRFG. GO and KEGG enrichment analyses demonstrated the molecular mechanism of DE-MRFG. Figures 1D–F show that DE-MRFG is enriched in iron ion transport, NADPH oxidase complex, and superoxide-generating NAD(P)H oxidase activity. Combined with the KEGG analysis results (Fig. 1G), the function and signaling pathways of DE-MRFG were enriched in oxidation-reduction and cancer, such as unsaturated fatty acid biosynthesis and thyroid cancer.

3.2 A risk model based on m6A-regulated ferroptosis genes

Univariate Cox regression results indicated that three genes, including TFRC, RGS4, and FTH1, are related to the survival state of patients with LC (Fig. 2A). Based on these three genes, LASSO regression established a risk model using the LASSO regression formula: Risk score = 0.097 × TFRC + 0.269 × RGS4 + 0.303 × FTH1. This formula assigned a risk score to each LC sample in the TCGA-HNSC dataset (108 samples with survival information). We divided all LC samples into two groups with an optimal risk score threshold of 3.44. In GSE65858, the risk score of LC samples was calculated and classified into two groups according to the optimal risk score threshold of 5.88 to verify the risk model. Additionally, KM survival analysis suggested that the low-risk group usually combined with a better survival state in the TCGA and GSE65858 datasets (Figs. 2B-C). ROC analysis demonstrated that the risk model accurately predicted OS (Figs. 2D-E). The survival state distribution maps (Figs. 2F-G) suggest that shorter survival times and more deaths are usually combined with higher risk scores. The three gene expressions are shown in the heatmap. The risk score in this heatmap was sorted ascending from the left to the right (Figs. 2H–I). These three genes showed higher expression in the high-risk groups.

3.3 Risk model was related to skin and protein interaction

To further investigate the potential influence of DEGs in the high- and low-risk groups, 5,754 DEGs, including 2,225 up- and 3,529 down-regulated genes, were found (Fig. 3A, Table S4). GO and KEGG enrichment analysis revealed that the most enriched GO items and pathways were ‘epidermis development,’ ‘collagen-containing extracellular matrix,’ ‘receptor ligand activity,’ and ‘neuroactive ligand-receptor interaction’ (Figs. 3B-C).

3.4 The performance of the risk model in different clinical characteristics

The distribution of LC samples with different risk scores in clinical characteristics indicated that the risk scores of patients with LC have no significant differences in most features, such as gender, tumor stage, pathological N stage, and pathological T stage, except for age (Fig. 4A). Figure 4B displays the TFRC, RGS4, and FTH1 expressions in different clinical characteristics. Furthermore, the KM survival curve of different groups with different clinical characteristics demonstrated that people in the low-risk group had better OS in males, older age, tumor stage III–IV, pathological N2_N3, pathological T2, and pathological T3 (Figs. 4C–H). Figure S1 depicts the results without significant differences.

3.5 Establishment of a nomogram

Univariate Cox regression analysis suggested that risk score, gender, and Pathologic_N were associated with the OS of patients with LC (Fig. 5A). According to multivariate Cox regression analysis (Fig. 5B), risk score, gender, and Pathologic_N were considered independent prognostic factors to build a nomogram (Fig. 5C) that can predict the OS of patients with LC. Calibration curves and ROC analyses proclaimed the good comportment of nomograms for predicting OS (Figs. 5D–G).

3.6 The risk model predicted the effect of immune therapy

SsGSEA was used to determine immune infiltration in the LC. Figure 6A reveals that the ssGSEA scores of six immune cells, including neutrophils, mast cells, and Memory B cells, differed significantly between the high- and low-risk groups. Additionally, Figs. 6B–D present the relationship between the three prognostic genes (TFRC, RGS4, and FTH1) and the six immune cells. Prognostic genes were positively related to the other five immune cells, except for neutrophils. Finally, stromal, immune, and ESTIMATE scores were calculated (Figs. 6E-F). The stromal and ESTIMATE scores were dramatically higher in the high-risk group than in the low-risk group.

To probe the mutation in different groups, the ‘maftools’ package was used. Figures 6G–H depict the top 20 mutant genes. The High-risk group had 3,905 mutated genes, which were present in all samples (Fig. 6G). The low-risk group contained 9,030 genes, with mutations occurring in 98.85% of the samples (Fig. 6H). TP53 had the highest mutation rate in both groups.

Based on the association analyses of immune checkpoints and risk score, only PD-L1 had a significantly positive relationship with risk score (Fig. 6I), with dramatically up-regulated expression in a high-risk group (Fig. 6J). KM survival analysis indicated better OS of patients with high-expressed PD-L1, although the result did not have statistical significance (Fig. 6K). Then, the TIDE, dysfunction, and exclusion scores of all LC samples were calculated. Moreover, the relationship between immunotherapy response and risk score was analyzed. Figures 6L-M indicate a significantly positive relationship between the risk score, TIDE, and exclusion score.

3.7 Chemotherapy drugs were related to the risk model

The IC₅₀ values of 20 chemotherapeutics were calculated for each LC sample to estimate the sensitivity of chemotherapy drugs in different groups. Nine drugs, such as Cisplatin, had a significantly higher IC₅₀ in the high-risk group than in the low-risk group, while 10 drugs, such as Erlotinib, had a lower one (Fig. 7). Moreover, the correlation analysis between risk score and IC₅₀ in 19 drugs indicated a strong relationship (Fig. 8). For example, Bosutinib was positively related to the risk score, while Docetaxel was negatively related to the risk score.

3.8 Expression of TFRC, RGS4 and FTH1 in laryngeal carcinoma

In this study, we analyzed the TFRC, RGS4, and FTH1 expression levels in LC cell lines and NHBEC. Our findings revealed a significant increase in TFRC, RGS4, and FTH1 expressions in these cell lines than in NHBEC (Figs. 9A-B). Additionally, we examined 30 pairs of laryngeal carcinomas and paracancerous tissues and utilized qRT-PCR to assess TFRC, RGS4, and FTH1 expression. The results demonstrated significantly higher TFRC, RGS4, and FTH1 expression in laryngeal carcinoma tissues than in paracancerous tissues (Fig. 9C). Furthermore, KM analysis demonstrated a significant correlation between TFRC and FTH1 expression and patient prognosis, whereas RGS4 expression was not significantly correlated with patient prognosis (Fig. 9D).

m6A modification has recently received significant attention in cell biology (31). Extensive research has demonstrated that m6A modification is pivotal for regulating RNA stability, transcription, and translation (32). Moreover, recent investigations have revealed a strong association between m6A modification and ferroptosis, offering a novel framework for comprehending cellular iron metabolism and regulation (33). Furthermore, m6A modification has been identified as intimately linked to tumorigenesis and disease progression in LC research (12, 34). Subsequent research has revealed that the m6A modification may influence the incidence and therapeutic efficacy of LC by regulating pivotal gene expression (35, 36). In summary, m6A modification is important in cellular biology and cancer investigation, providing valuable insights into the etiology and treatment approaches for associated diseases (37). Our study identified three model genes (TFRC, RGS4, and FTH1) from the DE-MRFG dataset to construct a risk model using the univariate Cox and LASSO regression algorithms. The analysis of immune reaction in lung cancer revealed a significant association between the immune checkpoint PD-L1 and the risk score, with PD-L1 up-regulated in the high-risk group. Additionally, the TIDE and ESTIMATE algorithms demonstrated a positive correlation between risk and the TIDE or ESTIMATE scores.

Previous studies have found that the LINC00888/miR-378g/TFRC pathway may lead to the development of LC cells, making it a potential therapeutic target for LC (38). TFRC knockdown strengthened the inhibitory effect of si-LINC00888 on LC cells' malignant properties. FTH1 plays a vital role in ferroptosis resistance in HNSCCs and provides clues to target HNSCCs resistant to ferroptosis induction and/or EGFR inhibition (39). In this study, we unveiled the differential expression of ferroptosis genes (TFRC, RGS4, and FTH1) associated with m6A in laryngeal carcinoma. These genes are potential prognostic indicators of patient survival. Using multivariate Cox regression analysis, we verified that the risk score, gender, and Pathologic_N are independent prognostic factors that can be utilized to construct a nomogram to predict the OS of patients with laryngeal carcinoma. The effectiveness of this nomogram for OS prediction was demonstrated using calibration curves and ROC analyses.

Recent investigations have suggested that these genes may influence immune cell response to immunotherapy and infiltration. For instance, heightened TFRC expression may be associated with unfavorable outcomes in immunotherapy (40). Similarly, aberrant RGS4 and FTH1 expression may disrupt the infiltration of immune cells, consequently affecting the immune response against tumors. Further research is warranted to elucidate the underlying mechanisms and their implications for treating laryngeal carcinoma.

External immunohistochemical analysis showed that TFRC and SLC7A11 were more highly expressed in tumor tissues than in para-cancer normal tissues, and their expression levels were significantly associated with a more advanced stage and worse cancer-specific survival. A prognostic signature including SERPINE1, EDIL3, RGS4, and MATN3 (SERM signature) was constructed to predict OS, DFS, and drug sensitivity in gastric cancer. A total of 14,742 DEGs were identified in the TCGA-HNSC dataset, and 5,832 LC-related genes were identified using WGCNA. A set of intersection genes was identified and designated as 97 DE-FG in LC by intersecting these 5,832 LC-related genes with the 14,742 DEGs and 567 ferroptosis-related genes. Among these DE-FG genes, 83 were correlated with m6A-related genes by correlation analysis. GO and KEGG enrichment analyses were performed on DE-MRFGs. Univariate Cox regression analysis revealed a significant association between the survival status of patients with LC and three specific genes, namely TFRC, RGS4, and FTH1. Subsequently, a risk model was constructed using LASSO regression based on three genes. KM survival analysis demonstrated that individuals in the low-risk group exhibited more favorable survival outcomes. Furthermore, the distribution maps of the survival status indicated that higher risk scores were typically associated with shorter survival times and a higher number of deceased patients. The ssGSEA scores of the six immune cell types differed significantly between the high- and low-risk groups. The relationship between the three prognostic genes (TFRC, RGS4, and FTH1) and six immune cells was investigated. Except for neutrophils, the prognostic genes were positively correlated with the other five immune cells. A total of 3,905 genes were mutated in the high-risk group, and these mutations were present in all samples. TP53 exhibited the highest mutation rate in both groups. Additionally, PD-L1 displayed a significantly positive association with the risk score, as evidenced by its up-regulated expression in the high-risk group. A correlation analysis of the risk score and IC₅₀ of the 19 chemotherapy drugs indicated a strong relationship between them. Additionally, qPCR analysis of 30 pairs of LC and paracancerous tissue samples revealed consistently high expression of TFRC, RGS4, and FTH1 in LC tissues and cell lines. These gene expression levels were correlated with the prognosis of patients with LC. Despite its dynamic nature, current cancer research has some limitations. However, we have addressed these shortcomings by connecting with clinical research, validating our findings in 30 clinical LC cases and three LC cell lines, and analyzing the correlation between target gene expression and prognosis. In future investigations, we will prioritize the examination of TFRC, RGS4, and FTH1 to delve deeper into their underlying mechanisms.

LC Laryngeal Cancer

WGCNA Weighted Gene Co-expression Network Analysis

DE-MRFG m6A-related ferroptosis genes

LASSO Least Absolute Shrinkage and Selection Operator

ROC Receiver Operating Characteristic

OS Overall Survival

TIDE Tumour Immune Dysfunction and Exclusion

ESTIMATE Estimation of STromal and Immune cells in MAlignant Tumors using Expression

ROS Reactive Oxygen Species

NHBEC Normal Human Bronchial Epithelial Cells

TCGA The Cancer Genome Atlas

GEO Gene Expression Omnibus

DEGs Differentially Expressed Genes

DE-FG Differentially Expressed Ferroptosis-related Genes

GO Gene Ontology

KEGG Kyoto Encyclopedia of Genes and Genomes

KM Kaplan-Meier

ssGSEA single-sample Gene Set Enrichment Analysis

GDSC Genomics of Drug Sensitivity in Cancer

Acknowledgements

Not applicable.

Authors’ contributions

XW: Conceptualization, Methodology, Funding acquisition, Writing–original draft, Writing–review & editing. HL: Conceptualization, Data curation, Project administration, Writing–review & editing. WZ: Investigation, Writing–review & editing. KL: Data curation, Writing–review & editing. YW: Formal Analysis, Methodology, Writing–review & editing. JZ: Methodology, Writing–review & editing. XQ: Formal Analysis, Writing–review & editing. JW: Data curation, Supervision, Writing–review & editing. AL: Supervision, Validation, Writing–review & editing. GS: Formal Analysis, Writing–review & editing. YY: Data curation, Writing–review & editing. YW: Formal Analysis, Writing–review & editing. YS: Investigation, Writing–review & editing.

Funding

This work was supported by the National Natural Science Foundation of China (No. 82203114), and Technology Incubation Fund and Talent Program Projects of Shaanxi Provincial People’s Hospital (2023JY-02).

Data availability

The data used in this study came from the TCGA and GSE65858 databases, and this study was in full compliance with the published guidelines of public databases. Accordingly, no additional informed consent was provided.

Competing interests

All authors declare that they have no confict of interest.

Ethics approval and consent to participate

According to the protocol approved by the Ethical Review Committee of the Second Affiliated Hospital of Harbin Medical University and Shaanxi Provincial People's Hospital, patient tissue samples were obtained with informed consent. All patients provided signed informed consent (2022–K268).

Consent for publication

Not applicable.

Author details

¹Department of Otorhinolaryngology, Head and Neck Surgery, Shaanxi Provincial People's Hospital, 256 Youyi Road, Xi’an 710000, China

² Department of Otorhinolaryngology, Head and Neck Surgery, The Second Affiliated Hospital, Harbin Medical University, 246 Xuefu Road, harbin 150000, China

Supplementary material

The Supplementary Material for this article can be found online.

Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin (2021) 71(1):7-33. doi: 10.3322/caac.21654
Han J, Sumer BD. The Changing Demographics and Treatment of Larynx Cancer. Ann Surg Oncol (2021) 28(12):6927-8. doi: 10.1245/s10434-021-10657-z
Deng K, Yao J, Zeng S, Wen M, Huang J, Zhu T, et al. The effect of surgery plus intensity-modulated radiotherapy on treatment in laryngeal cancer: A clinical retrospective study. J Cancer Res Clin Oncol (2022) 148(2):517-25. doi: 10.1007/s00432-021-03637-z
Huang Y, Gu M, Tang Y, Sun Z, Luo J, Li Z. Systematic review and meta-analysis of prognostic microRNA biomarkers for survival outcome in laryngeal squamous cell cancer. Cancer Cell Int (2021) 21(1):316. doi: 10.1186/s12935-021-02021-8
Hoesseini A, Dronkers EAC, Sewnaik A, Hardillo JAU, Baatenburg de Jong RJ, Offerman MPJ. Head and neck cancer patients' preferences for individualized prognostic information: a focus group study. BMC Cancer (2020) 20(1):399. doi: 10.1186/s12885-020-6554-8
Lechner M, Liu J, Masterson L, Fenton TR. HPV-associated oropharyngeal cancer: epidemiology, molecular biology and clinical management. Nat Rev Clin Oncol 2022;19(5):306-27. doi: 10.1038/s41571-022-00603-7
He X, Xu C. Immune checkpoint signaling and cancer immunotherapy. Cell Res (2020) 30(8):660-9. doi: 10.1038/s41422-020-0343-4
Sendinc E, Shi Y. RNA m6A methylation across the transcriptome. Mol Cell (2023) 83(3):428-41. doi: 10.1016/j.molcel.2023.01.006
Ma S, Chen C, Ji X, Liu J, Zhou Q, Wang G, et al. The interplay between m6A RNA methylation and noncoding RNA in cancer. J Hematol Oncol (2019) 12(1):121. doi: 10.1186/s13045-019-0805-7
Shi H, Wei J, He C. Where, When, and How: Context-Dependent Functions of RNA Methylation Writers, Readers, and Erasers. Mol Cell (2019) 74(4):640-50. doi: 10.1016/j.molcel.2019.04.025
Huang H, Weng H, Chen J. m(6)A Modification in Coding and Non-coding RNAs: Roles and Therapeutic Implications in Cancer. Cancer Cell (2020) 37(3):270-88. doi: 10.1016/j.ccell.2020.02.004
Wang X, Tian L, Li Y, Wang J, Yan B, Yang L, et al. RBM15 facilitates laryngeal squamous cell carcinoma progression by regulating TMBIM6 stability through IGF2BP3 dependent. J Exp Clin Cancer Res (2021) 40(1):80. doi: 10.1186/s13046-021-01871-4
Tang D, Chen X, Kang R, Kroemer G. Ferroptosis: molecular mechanisms and health implications. Cell Res (2021) 31(2):107-25. doi: 10.1038/s41422-020-00441-1
Jiang X, Stockwell BR, Conrad M. Ferroptosis: mechanisms, biology and role in disease. Nat Rev Mol Cell Biol (2021) 22(4):266-82. doi: 10.1038/s41580-020-00324-8
Chen X, Kang R, Kroemer G, Tang D. Broadening horizons: the role of ferroptosis in cancer. Nat Rev Clin Oncol (2021) 18(5):280-96. doi: 10.1038/s41571-020-00462-0
Wang Z, Dai Z, Zheng L, Xu B, Zhang H, Fan F, et al. Ferroptosis Activation Scoring Model Assists in Chemotherapeutic Agents' Selection and Mediates Cross-Talk With Immunocytes in Malignant Glioblastoma. Front Immunol (2021) 12:747408. doi: 10.3389/fimmu.2021.747408
Niu B, Liao K, Zhou Y, Wen T, Quan G, Pan X, Wu C. Application of glutathione depletion in cancer therapy: Enhanced ROS-based therapy, ferroptosis, and chemotherapy. Biomaterials (2021) 277:121110. doi: 10.1016/j.biomaterials.2021.121110
Bersuker K, Hendricks JM, Li Z, Magtanong L, Ford B, Tang PH, et al. The CoQ oxidoreductase FSP1 acts parallel to GPX4 to inhibit ferroptosis. Nature (2019) 575(7784):688-92. doi: 10.1038/s41586-019-1705-2
Wang Y, Wu X, Ren Z, Li Y, Zou W, Chen J, Wang H. Overcoming cancer chemotherapy resistance by the induction of ferroptosis. Drug Resist Updat (2023) 66:100916. doi: 10.1016/j.drup.2022.100916
Wichmann G, Rosolowski M, Krohn K, Kreuz M, Boehm A, Reiche A, et al. The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer. Int J Cancer (2015) 137(12):2846-57. doi: 10.1002/ijc.29649
Huang C, He J, Dong Y, Huang L, Chen Y, Peng A, Huang H. Identification of Novel Prognostic Markers Associated With Laryngeal Squamous Cell Carcinoma Using Comprehensive Analysis. Front Oncol (2021) 11:779153. doi: 10.3389/fonc.2021.779153
Zhang X, Zhang S, Yan X, Shan Y, Liu L, Zhou J, et al. m6A regulator-mediated RNA methylation modification patterns are involved in immune microenvironment regulation of periodontitis. J Cell Mol Med (2021) 25(7):3634-45. doi: 10.1111/jcmm.16469
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics (2008) 9:559. doi: 10.1186/1471-2105-9-559
Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141
Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw (2010) 33(1):1-22.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747-56. doi: 10.1101/gr.239244.118
Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform (2021) 22(6). doi: 10.1093/bib/bbab260
Wang X, Yu B, Jin Q, Zhang J, Yan B, Yang L, et al. Regulation of laryngeal squamous cell cancer progression by the lncRNA RP11-159K7.2/miR-206/DNMT3A axis. J Cell Mol Med (2020) 24(12):6781-95. doi: 10.1111/jcmm.15331
Zhuang H, Yu B, Tao D, Xu X, Xu Y, Wang J, et al. The role of m6A methylation in therapy resistance in cancer. Mol Cancer (2023) 22(1):91. doi: 10.1186/s12943-023-01782-2
Pan J, Huang T, Deng Z, Zou C. Roles and therapeutic implications of m6A modification in cancer immunotherapy. Front Immunol (2023) 14:1132601. doi: 10.3389/fimmu.2023.1132601
Xia H, Wu Y, Zhao J, Cheng C, Lin J, Yang Y, et al. N6-Methyladenosine-modified circSAV1 triggers ferroptosis in COPD through recruiting YTHDF1 to facilitate the translation of IREB2. Cell Death Differ (2023) 30(5):1293-304. doi: 10.1038/s41418-023-01138-9
Yang L, Yan B, Qu L, Ren J, Li Q, Wang J, et al. IGF2BP3 Regulates TMA7-mediated Autophagy and Cisplatin Resistance in Laryngeal Cancer via m6A RNA Methylation. Int J Biol Sci (2023) 19(5):1382-400. doi: 10.7150/ijbs.80921
Guo H, Han Q, Guan X, Li Z, Wang Y, He L, et al. M6A reader YTHDF1 promotes malignant progression of laryngeal squamous carcinoma through activating the EMT pathway by EIF4A3. Cell Signal (2023) 114:111002. doi: 10.1016/j.cellsig.2023.111002
Li J, Cao H, Yang J, Wang B. CircCDK1 blocking IGF2BP2-mediated m6A modification of CPPED1 promotes laryngeal squamous cell carcinoma metastasis via the PI3K/AKT signal pathway. Gene (2023) 884:147686. doi: 10.1016/j.gene.2023.147686
Wang L, Dou X, Chen S, Yu X, Huang X, Zhang L, et al. YTHDF2 inhibition potentiates radiotherapy antitumor efficacy. Cancer Cell (2023) 41(7):1294-308 e8. doi: 10.1016/j.ccell.2023.04.019
Cui P, Dai X, Liu R, Cao H. LncRNA LINC00888 upregulation predicts a worse survival of laryngeal cancer patients and accelerates the growth and mobility of laryngeal cancer cells through regulation of miR-378g/TFRC. J Biochem Mol Toxicol (2021) 35(10):e22878. doi: 10.1002/jbt.22878
Liu S, Yan S, Zhu J, Lu R, Kang C, Tang K, et al. Combination RSL3 Treatment Sensitizes Ferroptosis- and EGFR-Inhibition-Resistant HNSCCs to Cetuximab. Int J Mol Sci (2022) 23(16). doi: 10.3390/ijms23169014
Yang C, Li J, Guo Y, Gan D, Zhang C, Wang R, et al. Role of TFRC as a Novel Prognostic Biomarker and in Immunotherapy for Pancreatic Carcinoma. Front Mol Biosci (2022) 9:756895. doi: 10.3389/fmolb.2022.756895
Zhang J, Chen S, Wei S, Cheng S, Shi R, Zhao R, et al. CircRAPGEF5 interacts with RBFOX2 to confer ferroptosis resistance by modulating alternative splicing of TFRC in endometrial cancer. Redox Biol (2022) 57:102493. doi: 10.1016/j.redox.2022.102493
Hu Y, Zheng M, Wang S, Gao L, Gou R, Liu O, et al. Identification of a five-gene signature of the RGS gene family with prognostic value in ovarian cancer. Genomics (2021) 113(4):2134-44. doi: 10.1016/j.ygeno.2021.04.012
Hu ZW, Wen YH, Ma RQ, Chen L, Zeng XL, Wen WP, Sun W. Ferroptosis Driver SOCS1 and Suppressor FTH1 Independently Correlate With M1 and M2 Macrophage Infiltration in Head and Neck Squamous Cell Carcinoma. Front Cell Dev Biol (2021) 9:727762. doi: 10.3389/fcell.2021.727762
Jia X, Chen B, Li Z, Huang S, Chen S, Zhou R, et al. Identification of a Four-Gene-Based SERM Signature for Prognostic and Drug Sensitivity Prediction in Gastric Cancer. Front Oncol (2021) 11:799223. doi: 10.1016/j.ygeno.2021.04.012

No competing interests reported.

FigureS1.tif
Figure S1. Results without significant differences.
TableS1.docx
TableS2.xls
TableS3.xls
TableS4.xls

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of m6A-Regulated Ferroptosis Biomarkers for Prognosis in Laryngeal Cancer

Status:

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

Introduction

Materials and methods

2.1 Data source

2.2 Function analysis of differentially expressed m6A-related ferroptosis genes (DE-MRFG)

2.3 Establishment and verification of a risk model

2.4 Enrichment analysis of DEGs in different groups

2.5 Relationship between risk model and clinical characteristics

2.6 A nomogram model of LC

2.7 The relationship of risk model and immunotherapy

2.8 Prediction of Chemotherapy Drug Sensitivity

2.9 Tissue samples and Cell culture

2.10 qRT-PCR

2.11 Western blot

2.12 Statistical analysis

Results

3.1 DE-MRFG in LC was associated with oxidation-reduction and cancer

3.2 A risk model based on m6A-regulated ferroptosis genes

3.3 Risk model was related to skin and protein interaction

3.4 The performance of the risk model in different clinical characteristics

3.5 Establishment of a nomogram

3.6 The risk model predicted the effect of immune therapy

3.7 Chemotherapy drugs were related to the risk model

3.8 Expression of TFRC, RGS4 and FTH1 in laryngeal carcinoma

Discussion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1