Identification and validation of a ferroptosis-related gene to predict survival outcomes and the immune microenvironment in lung adenocarcinoma
Cancer Cell International volume 22, Article number: 292 (2022)
Lung adenocarcinoma (LUAD) is a leading cause of cancer-related death worldwide. Ferroptosis, a form of cell death characterized by iron-dependent lipid peroxidation. However, the involvement of ferroptosis in the regulation of immune cell infiltration and its immunotherapeutic efficacy in LUAD remain unclear.
The Cancer Genome Atlas (TCGA) LUAD cohort was used to assess the survival prognosis of FRGs and construct a seven-gene risk signature. Correlation tests, difference tests, and a cluster analysis were performed to explore the role of FRGs in the immune microenvironment and their immunotherapeutic efficacy in LUAD. The effects of FRGs on LUAD cells were assessed by Western blot, iron assay, and lipid peroxidation assay.
The seven-gene risk signatures of patients with LUAD were established and validated. FRG clustering based on 70 differentially expressed FRGs was associated with the immune microenvironment and indicated potential immune subtypes of LUAD. The seven-gene risk signature was an independent prognostic factor for LUAD and was used to divide the LUAD cohort into a high-risk and a low-risk group. Immunocyte infiltration levels, immune checkpoints, and immunotherapy response rates were significantly different between the two groups. Patients with high risk scores had lower overall levels of immunocyte infiltration but higher immunotherapy response rates. The key gene ribonucleotide reductase subunit M2 (RRM2) was associated with LUAD prognosis, which may be related to its ability to regulate the infiltration levels of activated mast cells and activated CD4 memory T cells. In addition, RRM2 was involved in ferroptosis, and its expression was up regulated in lung cancer tissues and the LUAD cell lines. Silencing RRM2 can inhibit the proliferation and induce ferroptosis of H1975 cells suggesting that silencing RRM2 could promote ferroptosis in H1975 cells.
Our results revealed RRM2 as a promising biomarker and therapeutic target associated with tumor immune infiltration in patients with LUAD.
Pulmonary carcinoma is the most serious malignancy, representing the leading cause of cancer deaths in developed and developing countries worldwide . Lung adenocarcinoma (LUAD) is the most common histological type of non-small cell pulmonary carcinoma , and compared with other lung cancer subtypes, LUAD is more closely associated with genomic changes and is more heterogeneous . Although mortality from pulmonary carcinoma has decreased in recent years, it still causes more deaths than breast, prostate, colorectal, and intracranial tumors combined . More than half of patients with LUAD are diagnosed at advanced stages or with metastatic disease , and the patients have a poor prognosis with a low 5-year survival rate .
Ferroptosis, a recently discovered form of programmed cell death characterized by iron-dependent lipid peroxidation, is associated with a variety of diseases, particularly cancer; therefore, modulating the occurrence of ferroptosis in cancer cells could be a potential strategy for cancer therapy . According to the literature, mitochondrial-induced cysteine starvation, endoplasmic reticulum-related oxidative stress, lysosome dysfunction, and lipid peroxidation related to Golgi stress contribute to the induction of ferroptosis in LUAD . Researchers have also identified numerous genes involved in the organelle changes that induce ferroptosis. The regulatory mechanisms of these ferroptosis-related genes (FRGs) in the occurrence and development of pulmonary carcinoma have attracted increasing attention . Bufu Tang et al. found that FRGs affected pulmonary carcinoma progression and tumor immunocyte infiltration. Ribonucleotide reductase subunit M2 (RRM2), a rate-limiting protein for deoxynucleoside triphosphate (dNTP) synthesis, inhibitedM1 macrophage polarization and promoted M2 macrophage polarization. However, ferrostatin-1 treatment effectively rebalanced macrophage polarization mediated by RRM2 inhibitors . Huang et al. reported that aldo–keto reductase family 1 member C1 (AKR1C1) induced ferroptosis through multiple pathways and was associated with various cancer-infiltrating immune cells . A study by Min Wang et al. showed that long noncoding RNA (lncRNA) and competing endogenous RNA (ceRNA) networks also play important roles in tumorigenesis and ferroptosis. Endogenous microRNA 6852 (miR-6852) inhibited cell growth, and by acting as an miR-6852 sponge, LINC00336 increased the expression of cystathionine-β-synthase (CBS), thereby inhibiting ferroptosis . In addition to chromosomal genes, ferroptosis-related lncRNAs in the cytoplasm can also inhibit cancer by activating the p53 pathway. For example, studies by Chao Mao et al. showed that lncRNA P53RRA could bind to protein-binding protein 1(G3BP1) to activate G3BP1. The P53RRA–G3BP1 interaction replaced p53 in the G3BP1 complex, causing more p53 to be retained in the nucleus and leading to cell cycle arrest, apoptosis, and ferroptosis . Luo et al. summarized the regulatory mechanism of ferroptosis and the effect of ferroptosis on tumor cell metabolism and antitumor immunity .
However, the significance of FRGs in the immune microenvironment and their potential role in immunotherapy for LUAD remains unclear, and thus FRGs require further investigation. In this study, we identified seven FRGs closely related to the immune microenvironment of LUAD, of which RRM2 affected the prognosis of LUAD by regulating the infiltration levels of activated mast cells and activated CD4 memory T cells. The expression of RRM2 also affected the response rate to LUAD immunotherapy. We also demonstrated the expression RRM2 is elevated in lung cancer tissues and knockdown of RRM2 induces ferroptosis in lung cancer cells H1975.We propose potential immune subtypes of LUAD based on the FRG cluster analysis that may lead to improved treatment for LUAD.
Materials and methods
Obtaining FRGand gene expression data for LUAD
A list of FGRs was downloaded from the public database FerrDb (http://www.zhounan.org/ferrdb/); after the noncoding genes were deleted, 382 FRGs remained. Transcriptome data for LUAD in fragments per kilobase of sequence per million mapped reads (FPKM) format were downloaded from The Cancer Genome Atlas (TCGA) (https://ancergenome.nih.gov/). The data included a total of 56,530 genes from 594 samples (535 tumor samples and 59 paracancerous samples); 492 of the samples had complete clinical information, including the survival time, survival status, age, gender, stage, and tumor/node/metastasis (TNM) state. Data from the GSE310219 cohort (54,675 genes from 307 early LUAD samples, 292 of which had complete clinical information) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).The clinical data of patients with LUAD included in this study are shown in Table 1.
Differentially expressed genes (DEGs) between LUAD tissues and adjacent tissues
The mRNA data for the TCGA-LUAD and GSE310219 cohorts were standardized. The transcription sequences were then compared with the FRGs, which were analyzed using limma (an R software package). Differentially expressed genes (DEGs) with a log2 absolute value (FC) of > 1and an adjusted P-value of < 0.05 were selected.
Establishment and validation of a risk signature
Univariate Cox regression analysis was used to screen genes related to overall survival (OS), and multivariate Cox regression analysis was used to identify genes with independent prognostic effects. A P-value of < 0.05was considered statistically significant. The genes screened using multivariate Cox regression were used to construct a risk signature; subsequently, the following predictive risk-scoring model was established combining the regression coefficient (β) with the gene expression level: RiskScore = (β1 * expression level of ALOX12B) + (β2 * expression level of DDIT4) + (β3 * expression level of SLC7A5) + (β4 * expression level of TRIB3) + (β5 * expression level of IL33) + (β6 * expression level of RRM2) + (β7 * expression level of CAV1). The establishment of the risk signature was achieved using the R software packages ‘survival’ and ‘survminer’.
Prognostic potential analysis of risk signature
First, univariate and multivariate Cox regression were used to analyze whether the prognostic indicators were independent of other traditional clinical features, such as age, gender, and TNM state. The hazard ratio (HR) and 95% confidence interval (CI)of each variable were measured, and a P-value of < 0.05 was considered statistically significant. Then, the patients with LUAD in TCGA and the GEO database were divided into a high-risk group and a low-risk group with the median risk score as the boundary. The difference in the OS between the high-risk group and the low-risk group was analyzed using the Kaplan–Meier (K–M) curve. The prognostic potential and significance of the risk model were tested by a receiver operating characteristic (ROC) curve, and the area under the curve (AUC) was used to evaluate the specificity and sensitivity of the model. This analysis was conducted using the ‘survive’ package in R, and the ‘ggDCA’ package was used for the decision curve analysis (DCA) to analyze the net benefit of the risk signature. P < 0.05 represented statistical significance.
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the R packages ‘clusterProfiler’, ‘org.Hs.eg.db’, and ‘enrichplot’ to explore the differences in various molecular mechanisms between the high- and low-risk groups. GO analysis included biological processes (BPs), cell components (CCs), and molecular functions (MFs).
Construction of a predictive nomogram
Using the R packages ‘survival’ and ‘regplot’, all independent predictors of LUAD were identified and then combined with the survival information of patients with LUAD; nomograms were drawn to predict 1-,3-, and 5-year survival rates.
Analysis of tumor microenvironment
First, the StromalScore, ImmuneScore, and EstimateScore of each LUAD sample were calculated by the estimation algorithm, and the differences between the high- and low-risk groups were analyzed. Second, transcriptome data were quantitatively converted into the absolute abundance of specific cell types, including immune cells and stromal cells, using seven algorithms: CIBERSORT, TIMER, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC. Third, immunocytes with independent predictive significance for the prognosis of LUAD were identified by univariate and multivariate Cox regression analysis, and the correlation between the expression of each risk gene and the level of immunocyte infiltration was analyzed using a correlation test. Fourth, the TCGA-LUAD cohort was clustered according to the level of immunocyte infiltration, and the LUAD samples were divided into a high-infiltration group and a low-infiltration group. The differences in clinical characteristics between the two groups were compared using a difference test. Finally, the scores of 29 immune-related functions or pathways were calculated by single set gene set enrichment analysis (ssGSEA), and the differences in scores between the high- and low-risk groups were compared (statistical significance was defined as P < 0.05 for all tests).
Correlation analysis of immunotherapy
Tumor immune dysfunction and exclusion (TIDE) and microsatellite instability (MSI) scores of TCGA-LUAD participants were downloaded from the TIDE database (http://tide.dfci.harvard.edu/), and the tumor mutational burden (TMB) of patients with LUAD was downloaded using the R packages ‘TCGAbiolinks’ and ‘maftools’. Then, the Wilcoxon rank sum test was used to determine whether there were significant differences in the scores of five types of immunotherapies between the high- and low-risk groups. Moreover, the data of participants in the LUAD cohort (GSE126044) with anti-PDL1 treatment history were downloaded from the GEO database, and the differences in responses to anti-PDL1 immunotherapy between patients with high and low risk were compared. Additionally, the expression matrix of common immune checkpoints was extracted from the TCGA-LUAD cohort, and the expression differences between the high- and low-risk groups were compared. Finally, the correlation between first-line targeted therapy driving genes and risk genes was analyzed using a correlation test. P < 0.05was considered statistically significant.
The risk genes were compared with the reference data set of the CMap database (https://portals.broadinstitute.org/cmap/), and the drugs were ranked according to the enrichment of risk genes in the reference gene expression profile. The drugs that significantly inhibited the expression of high-risk genes were screened for those with a P-value of < 0.05 and an enrichment score of < 0.The chemical structures of the selected drugs were searched in the PubChem database (https://pubchem.ncbi.nlm.nih.gov/).Using the ‘pRRophetic’ algorithm, a ridge regression model was constructed to predict the IC50sof these drugs according to the Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/) cell line expression profile and TCGA gene expression profile, and antitumor drugs with significantly lower IC50s for high-risk LUAD samples were screened.
Cell culture and transfection
The normal lung cell line BEAS-2B and seven lung cancer cell lines (H1299, A549, H460, H23, H838, PC-9, and H1975) were purchased from the American Type Culture Collection (ATCC) and maintained in 1640 medium containing 10% fetal bovine serum (FBS) (Gibco, Grand Island, NY, USA) and 1% penicillin–streptomycin (Gibco, Grand Island, NY, USA) and cultured at 37 °C with 5% CO2. The siRNA of RRM2 was purchased from Sangon Biotech (Shanghai, China) was used to silence the expression of RRM2. In this study, the sequence of RRM2 siRNA was sense: 5'- GCGAUUUAGCCAAGAAGUUTT-3'. H1975 cells were seeded in 12-well plates at a density of 6 × 104 cells/well and transfected 24 h later. RRM2 siRNA or negative control (NC) siRNA was transfected to a final concentration of 20 nM using an siRNA transfection reagent (Polyplus, France). Finally, the transfected cells were detected after48h.
Proteins were extracted from cells, and cell lysates were prepared using RIPA lysate with phenylmethylsulfonyl fluoride (PMSF) (Solarbio, Beijing, China). Protein quantification was performed using a BCA protein assay kit (Sangon Biotech, Shanghai, China). Proteins were then separated by 12% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS–PAGE) and transferred to nitrocellulose membranes. The membranes were blocked with 5% bovine serum albumin (BSA) for 1 h and then incubated with primary antibody overnight at 4 °C. The next day, after washing the membrane three times with tris-buffered saline with Tween 20 (TBST), the membrane was incubated for 1 h at room temperature with horseradish peroxidase-labeled secondary antibody (1:4000), followed by three washes with TBST. Finally, the colors were developed using BeyoECL Moon (Beyotime Biotechnology, Shanghai, China).
Detection of lipid peroxidation and apoptosis
C11-BODIPY 581/591 (10 μM; aBclone, Wuhan, China) was added to the H1975 cells, and the cells were incubated at 37 °C and 5% CO2 for 1 h. Then, the cells were washed twice with phosphate-buffered saline (PBS). After trypsin digestion, the cells were resuspended with PBS containing 5% FBS and analyzed by flow cytometry. An apoptosis detection kit (BD Biosciences) was used to prepare fluorescent dyes containing propidium iodide (PI) and fluorescein isothiocyanate (FITC) according to the manufacturer’s instructions. The H1975 cells were incubated in the dark at room temperature for 15 min and then analyzed by flow cytometry.
FerroOrange (Dojindo, China) was used to detect the ferrous level of the cells. Following the manufacturer’s instructions, cells were incubated with FerroOrange for 0.5 h. The fluorescence intensity was evaluated using a rotating disk super-resolution laser confocal microscope.
From January 2020 to January 2021, 12 samples of LUAD tissues and normal tissues in the Department of Cardiothoracic Surgery, Affiliated Hospital of Guangdong Medical University, were selected. The obtained tissue was transferred to a −80 °C refrigerator in an ice box for subsequent immunohistochemical (IHC) experiments. All included research participants were assessed by experts in the pathology department of our hospital, and none of the patients studied had received radiotherapy, chemotherapy, targeted therapy, or immunotherapy before surgery. All patients signed informed consent before surgery. The tissue was fixed with 4% paraformaldehyde and embedded in paraffin to prepare slices with a thickness of 5 μM. The slices were then dewaxed with xylene and dehydrated in gradient concentrations of alcohol. Subsequently, 0.01 M sodium citrate (pH 6.0) was used for antigen repair, and endogenous peroxidase was blocked by adding 0.3% hydrogen peroxide (H2O2) and incubating in 10% goat serum albumin for 30 min. The slices were then incubated overnight with RRM2 primary antibody at 4 °C.The next day, the samples were incubated with 3,3′diaminobenzidine (DAB) after incubating with HRP-conjugated anti-rabbit secondary antibody for 1 h. The sections were stained with Mayer hematoxylin, dehydrated, removed with xylene, sealed with neutral resin, and detected by a multi-functional microporous microscope (BIOTEK).
All data were generated, processed, and analyzed in R (version 4.1.1), R-Studio, and Strawberry Perl (220.127.116.11). Student’s t-test was used to determine the difference between the two groups, and the Wilcoxon rank sum test was used to identify differences for data from non-paired groups. P < 0.05 indicated that the difference was statistically significant. *P < 0.05, **P < 0.01, and ***P < 0.001.
Obtaining DEGs of LUAD
A total of 382 FRGs were downloaded from the public database FerrDb, and noncoding genes were deleted from the data set. Difference analysis showed that a total of 70 FRGs were differentially expressed between tumors and normal tissues; the characteristics of these FRGs in the TCGA-LUAD cohort are shown in Fig. 1A. In addition, all LUAD samples were clustered by consensus clustering based on 70 differentially expressed FRGs, and the best clustering effect was achieved when K = 4 (Additional file 1: Figure S1A–D). Surprisingly, significant differences were observed in OS between the four clusters (Additional file 1: Fig. S1E).
Construction and validation of prognostic model
As shown in the forest map, 15 DEGs associated with LUAD prognosis were identified by the R software package ‘survival’ in the TCGA cohort (Fig. 1B). Subsequently, seven FRGs with independent prognostic significance for patients with LUAD were obtained by multivariate Cox regression analysis on the 15 appealed DEGs; the locations of these FRGs in the differential volcanic map are shown in Fig. 1C. The Kaplan–Meier curve revealed the survival differences between the seven FRGs at different levels (Additional file 1: Fig. S1F–L). Ultimately, we established a seven-gene risk signature using the multivariate Cox regression coefficient, modeled as follows: Risk score of each patient = (0.262297470303534 * expression level of ALOX12B) + (0.0055627086031226 * expression level of DDIT4) + (0.00646442374221167 * expression level of SLC7A5) + (− 0.0131510931059613 * expression level of TRIB3) + (− 0.0300191156801198 * expression level of IL33) + (0.015384279700975 * expression level of RRM2) + (0.0046423481143834 * expression level of CAV1). According to the risk index ranking of 504 patients in TCGA-LUAD, patients were divided into a high-risk group and a low-risk group, with the median risk index as the cutoff between high and low risk. An independent prognostic analysis of the clinical characteristics of LUAD (age, gender, stage, and risk score) was conducted (Fig. 1D, E) after screening out the clinical information. The risk score and clinical stage were independent factors affecting the prognosis of LUAD. The clinical correlation chart in Additional file 3: Figure S3A shows the proportion of patients with high and low risk scores at each stage. Further analysis showed that the OS rate of TCGA-LUAD participants in the high-risk group was significantly lower than the OS rate of those in the low-risk group (P < 0.001) (Fig. 2A). Figure 2B shows that a difference in the OS rate was also present between the low-risk and high-risk groups in the GSE30219 cohort; the OS rate of the high-risk group was significantly lower than that of the low-risk group (P < 0.001). The ROC curve (Fig. 2C) shows that our risk signature had high sensitivity and specificity in predicting the OS rate of patients with LUAD. The AUCs were 0.707, 0.685, and 0.673 for 1-, 2-, and 3-year OS, respectively. We then added clinical features, including TIDE, tumor inflammation signature (TIS), age, and stage, to compare the accuracy of OS prediction. As shown in Fig. 2D, the risk index was the best predictor of OS. The DCA (Fig. 2E) showed that the net benefit of the risk index was the largest. Finally, a nomogram based on the two independent prognostic factors was constructed to quantify the survival probability of patients with LUAD in 1, 3, and 5 years (Fig. 2F).
We conducted an enrichment analysis to explore the potential biological characteristics of 70 DEGs. The 25 pathways shown in Additional file 2: Figure S2A were enriched in the KEGG enrichment analysis. The more prominent pathways, represented by red bubbles, included fluid shear stress and atherosclerosis, HIF-1 signaling pathways, lipid and atherosclerosis, and ferroptosis. GO enrichment analysis showed that the biological processes (including responses to oxidative stress and metal ions) and the cell components (including the apical plasma membrane, cell apex, and oxidoreductase complex) were significantly enriched (Additional file 2: Figure S2B). Furthermore, GSEA enrichment analysis showed that the high- and low-risk groups were enriched incompletely different functions or pathways. Among high-risk patients, small cell pulmonary carcinoma, melanoma, cancer pathways, thep53 signaling pathway, and pathogenic Escherichia coli infection were significantly enriched (Additional file 2: Figure S2C).
Difference of immune related indexes between LUAD patients
We calculate the StromalScore, ImmuneScore, and EstimateScore for each LUAD sample usingthe estimate algorithm. The t-tests revealed significant differences in StromalScore, ImmuneScore, and EstimateScore between high- and low-risk groups (Fig. 3A–C). To explore the influence of different levels of immune-related scores on the prognosis of patients, the StromalScore, ImmuneScore, and EstimateScore of all patients with LUAD were divided into high- and low-score groups by average. The OS rates of the two groups were compared using the Kaplan–Meier method; Fig. 3D shows that a higher EstimateScore represented a better prognosis. The same was true for ImmuneScore and StromalScore (Additional file 5: Figure S5A, B). Significant differences in risk scores were also observed among the four clusters (Additional file 3: Figure S3C).
LUAD immune subtype based on FRG clustering
In addition, the t-test revealed significant differences in immune-related scores among the four clusters obtained by the consensus clustering analysis (Fig. 3E–G), suggesting that the four clusters obtained for70 differentially expressed FRGs may serve as immune subtypes. To explore the relationship between the aforementioned four clusters and the two groups clustered according to the level of immunocyte infiltration, we drew a clinical correlation heat map, shown in Fig. 3H. Significant differences were observed in immunocyte infiltration among the four clusters. For example, the samples with high immunocyte infiltration in the composition of C3 were significantly more than those with low immunocyte infiltration, whereas the samples with low immunocyte infiltration in the C2 and C4 clusters were more. The mulberry diagram in Fig. 4A shows the corresponding relationship between immune infiltration level, FRG cluster, risk level, and other clinical features. The above results show that the four clusters identified by the clustering of 70 FRGs may represent the immune subtypes of LUAD.
Further analysis of immune cell infiltration in LUAD tumor microenvironment
The CIBERSORT algorithm was used to quantify the levels of 22 immune cells in each LUAD sample, and at-test was used to compare the levels of immune cells between the high-risk and low-risk patients. As shown in Fig. 4B, significant differences were observed in the levels of B cell memory, CD4 memory T cells, activated CD4 memory T cells, regulatory T cells (Tregs), resting NK cells, monocytes, M0 macrophages, resting dendritic cells, resting mast cells, and activated mast cells between high-risk and low-risk patients. The relative proportion of 22 immune cells in all LUAD samples is shown by a barplot (Additional file 3: Figure S3D). When Kaplan–Meier curves were used to compare the effects of immunocytes with different infiltration levels on the OS rates, we found that different infiltration levels of 12 typesof cells, such as activated CD4 memory T cells and resting mast cells, were associated with different OS rates (Additional file 4: Figure S4A–L). To learn more about the differences between immune cells in high- and low-risk patients, we analyzed the data using six other algorithms (TIMER, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC); Fig. 4B shows the different infiltration levels of immunocytes in each database.
Relationship between immunocyte infiltration and clinical characteristics
The correlation heat map shown in Fig. 4C reveals the relationship between the immune cells; most modules are blue, suggesting that the immune cell infiltration in the LUAD microenvironment is negatively correlated. To understand the correlation between immune cell infiltration and clinical characteristics, we calculated the immune cell infiltration level of all LUAD samples in TCGA using the ssGSEA algorithm and then downloaded the clinical information for all LUAD samples from TCGA using the R package ‘TCGAbiolinks’. The 580 LUAD samples with clinical information were clustered into two groups according to the level of immune cell infiltration and visualized by a heat map (Fig. 5A). The left portion, dominated by blue, represents the samples with less immunocyte infiltration, whereas the right portion, dominated by red, represents the samples with higher immunocyte infiltration. The clinical information, such as survival status, stage, gender, tumor location, and related driving gene mutation state, corresponded to the heat map below. In addition, to better understand the immune differences between high-risk and low-risk groups, we quantified 29 immune-related functions using the ssGSEA algorithm. As shown in Fig. 5B, 16 immune-related functions or pathways, such as Adcs and CCR, showed significant differences. When the Kaplan–Meier curve was used to compare the effects of different levels of immune-related functions on OS rates, we found that various levels of multiple immune-related functions or pathways predicted different OS rates. Figure 5C shows that varying levels of mast cell-related immune functions or pathways corresponded to a different prognosis for patients with LUAD. After univariate Cox regression analysis, we identified activated mast cells and activated CD4 memory T cells as prognostic immune cells of LUAD; these two cell types, which had high levels of infiltration, were adverse prognostic factors for LUAD (Fig. 6A). Further analysis confirmed that the expression levels of seven risk genes significantly affected the infiltration level of immunocytes in the LUAD microenvironment, and the infiltration level of activated mast cells, which was most related to the prognosis of LUAD, was affected by the expression levels of RRM2 and IL33 (Fig. 6B). The expression of RRM2 not only regulates the infiltration of activated mast cells but is also related to the infiltration of 19 immunocytes, such as CD4 memory T cells, dendritic cells, and macrophages (Fig. 6C).
Risk signature and immunotherapy
To explore the role of the risk signature in guiding immunotherapy, we identified16 patients with pulmonary carcinoma receiving anti-PDL1 immunotherapy in the GSE126044 cohort. Additional file 5: Figure S5G shows the differences in immunocyte infiltration and response to anti-PD-L1 immunotherapy between high-risk and low-risk patients. High-risk patients had lower immunocyte infiltration levels, but their response rates to anti-PD-L1 immunotherapy were higher. In addition, we found that the expression of RRM2 was higher in patients who had no response to anti-PDL1 treatment (P = 0.052) (Additional file 5: Figure S5F). We then downloaded the TIDE and MSI scores of the LUAD cohort from the TIDE database, and we downloaded the TMB data of patients with LUAD using the R packages ‘TCGAbiolinks’ and ‘maftools’. The Wilcoxon rank sum test revealed significant differences between the five immune therapy-related scores described above in high-risk and low-risk patients (Fig. 6D–H). The patients with LUAD with high TMB had worse OS (Additional file 5: Figure S5D). Subsequent correlation analysis showed that the risk score was positively correlated with TMB (R = 0.23, P = 0.000018) (Additional file 5: Figure S5C). The correlation heat map in Additional file 3: Figure S3B shows the distribution of clinical features in the high-risk and low-risk groups. The distribution of TMB, gender, stage, T, and N was significantly different in the high-risk and low-risk groups.
To study the relationship between the seven-gene risk signature and immunotherapy, we counted the mutations of seven key genes in patients with LUAD. As shown in the waterfall diagram in Additional file 5: Figure S5E, ALOX12B had the highest the mutation probability (18%), followed by RRM2 and IL33, and missense mutation was the most common form. To further explore the relationship between the expression of the seven risk genes and immunocyte infiltration, we conducted a correlation analysis, as shown in Fig. 7A, B. The expression level of RRM2 positively regulated the infiltration abundance of activated mast cells (R = 0.12) and activated CD4 memory T cells (R = 0.41), and the infiltration level of resting mast cells was significantly negatively correlated with the expression level of ALOX12B (R = −0.13) (Additional file 3: Figure S3F). When the ssGSEA algorithm was used to analyze the scores of immune-related functions or pathways, the scores of the checkpoints between high- and low-risk LUAD samples were significantly different. The checkpoints were divided into high and low groups using the average, and the Kaplan–Meier curve showed significant differences in survival (Fig. 7C). Next, we analyzed the expression levels of 49 common immune checkpoints in high-risk versus low-risk patients, as shown in the box plot in Fig. 7D. The expression levels of ALK, ROS1, CD44, and other 27 checkpoints were different. As shown in Fig. 7E, F, the expression levels of ALK and ROS1 decreased as the risk score increased. Furthermore, the correlation analysis showed that ALK and ROS1 were significantly associated with key risk genes. The expression of ROS1 was positively correlated with that ofIL33 and negatively correlated with that ofSLC7A5, TRIB3, and RRM2 (Fig. 7G). The expression of ALK was positively correlated with the expression of IL33 and CAV1and negatively correlated with the expression of ALOX12B, SLC7A5, and TRIB3 (Fig. 7H). The expression of two additional checkpoints, PD-1 and PD-L1, were also significantly correlated with multiple prognostic FRGs (Additional file 5: Figure S5H, I).
Overexpression of RRM2 in tumor tissues
To detect the expression of RRM2 in normal lung tissues and lung cancer tissues, we performed immunohistochemical experiments for comparison. The results showed that the expression of RRM2 in tumor tissues was higher than that in normal lung tissues (Fig. 8).
Silencing RRM2 induced ferroptosis in lung cancer cells
Ferroptosis is a recently discovered form of cell death characterized by lipid peroxidation accumulation and iron dependence. Further study on the specific mechanism of ferroptosis is expected to bring new prospects for cancer treatment. To determine the clinical relevance of RRM2 expression, we detected RRM2 expression in a normal lung cell line (BEAS-2B) and seven lung cancer cell lines (H1299, A549, H460, H23, H838, PC-9, and H1975). Western blotting showed that RRM2 protein was highly expressed in H1975 cells. In addition, we detected the expression of 4HNE and ACSL4 in different lung cancer cell lines. The results showed that 4HNE and ACSL4 were highly expressed in H1975 cells (Fig. 9A), indicating that H1975 cells may be more sensitive to ferroptosis. Therefore, we used H1975 cells for subsequent experiments. After H1975 cells were transfected with RRM2siRNA for 48 h, cell death was detected by flow cytometry. The results showed that silencing RRM2 could induce cell death (Fig. 9B). To further study the effect of silencing RRM2 on ferroptosis, Western blotting was used again to detect the expression of several proteins in H1975 cells with RRM2 siRNA transfected. The results showed that the expression of RRM2 in the silencing group was significantly decreased, indicating that silencing RRM2 was successfully. In the meantime, the expression of ferroptosis-related indicators like SLC7A11 and GPX4 were decreased, while ACSL4 were increased in H1975 cells after silencing RRM2 (Fig. 9C). Then, we continued to detect the levels of ferrous ion and lipid peroxidation in lung cancer cells after silencing RRM2. The results showed that silencing RRM2 could induce an increase in ferrous ion level (Fig. 9D) and lipid peroxidation accumulation (Fig. 9E) in H1975 cells. Ultimately, we demonstrated that silencing RRM2 induced ferroptosis in H1975 cells.
Screening drugs for high-risk FRGs
To predict effective therapeutic drugs for patients with high-risk LUAD, we constructed a ridge regression model to predict the drug IC50using the GDSC cell line expression profile and TCGA gene expression profile; we screened six antitumor drugs with significantly lower IC50 in high-risk LUAD: bosutinib, dasatinib, gefitinib, tipifarnib, docetaxel, and JNK inhibitor VIII (Additional file 6: Figure S6A–F). These six antitumor drugs could kill tumor cells in patients with high-risk LUAD. In addition, we compared the risk genes with the reference data set of the CMap database and obtained the correlation of 1309 drugs according to the enrichment of risk genes in the reference gene expression profile. With P < 0.05 as the standard, 158 drugs were found to be significantly enriched in high-risk genes; of these, 86 drugs promoted the expression of high-risk genes, and 72 drugs inhibited the expression of high-risk genes. Medrysone, phenoxybenzamine, vorinostat, thioguanosine, apigenin, and chrysin were selected as the six drug candidates enriched with more prominent inhibitors of high-risk gene expression, and their chemical structures were found in the PubChem database (Additional file 6: Figure S6G-L). The detailed information on these six drugs is shown in Table 2.
Recent studies have shown that ferroptosis is closely related to many diseases, such as cancer, blood diseases, neurological diseases, kidney diseases, and local ischemia–reperfusion injury . A growing number of studies investigating ferroptosis in cancer have revealed its potential as an immunotherapy strategy . LUAD is the most common pathological type of non-small cell lung cancer, and most patients are diagnosed at advanced stages, losing their opportunities for surgical treatment. Chemotherapy, radiotherapy, and traditional Chinese medicine are common treatment choices for patients with advanced LUAD. However, the 5-year survival rate is only 15% . The benefits of targeted therapy and immunotherapy in the treatment of advanced cancer patients have begun to change the treatment of cancer, providing new treatment methods for patients with advanced tumors. In this study, the seven-gene risk signature we established was closely related to the immune microenvironment of LUAD, and the immunocyte infiltration and immune function were significantly different in the high- and low-risk groups. Moreover, the two groups of patients with LUAD with high and low infiltration (revealed by immunocyte infiltration clustering) had different ferroptosis gene clusters and risk scores, which suggests that our seven-gene signature may provide guidance for LUAD immunotherapy. The immunotherapy analysis showed that patients with higher risk scores were more likely to benefit from immunotherapy. However, previous results showed that the overall infiltration level of immunocytes in high-risk patients was lower, perhaps because low immunocyte infiltration levels are related to a good prognosis, whereas increased infiltration levels of activated mast cells and activated CD4 memory T cells lead to a worse prognosis in patients with LUAD. Ultimately, patients with LUAD with lower infiltration levels were identified as high-risk by the seven-gene signature. Immunotherapy drugs can also activate the immune system to a greater extent; increase the infiltration of immunocytes related to a good prognosis and inhibit the activation of mast cells, CD4 memory T cells, and dendritic cells; and eventually improve the response rate to immunotherapy. These results show that FRGs regulate the tumor microenvironment of LUAD and influence the efficacy of immunotherapy. Gene mutation analysis in LUAD samples showed that seven FRGs had some degree of mutations; RRM2 and ALOX12B mutations were the most significant. Further analysis also found that the expression level of RRM2 was strongly correlated with the level of immunocytes in LUAD. In addition, we found that silencing RRM2 promoted elevated ferrous ion levels and lipid peroxidation accumulation and induced ferroptosis in LUAD. Therefore, we suggest that RRM2 may play a key role in LUAD (Fig. 10).
The inhibition of RRM2 reduces the level of dNTP, which is the basis for stable senescence-related cell growth arrest induced by oncogenes .Knocking down RRM2 significantly decreased proliferation during the S phase of the cell cycle [19, 20].Interestingly, researchers found that RRM2 regulated the expression of BCL-2 (a key determinant of cell apoptosis).Rahman et al.’s study found that RRM2 consumption significantly reduced the expression of Bcl-2 protein , and Jin et al. found that by overexpressing RRM2, the activity of the Bcl-2 signaling pathway was increased and the activity of the p53 signaling pathway was decreased . These findings suggest that RRM2 plays an important role in the S phase during DNA replication, which may have potential therapeutic significance. Zhong et al. and Li et al. studied the potential of RRM2 as a therapeutic target for gastric cancer and glioblastoma [23, 24]. Another study found that human papillomavirus E7 oncoprotein increased the expression of RRM2 to promote angiogenesis in cervical cancer and that inhibiting RRM2 activity may be a new therapeutic strategy for human cervical cancer . Most importantly, treatment methods related to RRM2 have greatly improved. Shao et al. found that RRM2 inhibitor could treat breast cancer [26, 27]. Another study showed that using gambogic acid to reduce the expression of RRM2 could improve the efficacy of gemcitabine in pancreatic cancer; gambogic acid caused pancreatic cancer cells to become sensitive to gemcitabine in vitro and in vivo by inhibiting the activation of the ERK/E2F1/RRM2 signaling pathway .
To date, few studies have explored the role of RRM2 in lung cancer. RRM2, which is highly expressed in lung cancer and related to a poor prognosis, regulates the immune microenvironment of LUAD, which is consistent with our results . Cai et al. detected significant changes in RRM2 methylation in patients with non-small cell lung cancer using bioinformatics and found a significant correlation between RRM2 expression and multiple immunocyte infiltration . Knockdown of RRM2 suppressed tumor growth in xenograft tumor models, and RRM2 deficiency increased CD8 + T cells in tumor tissues and the spleen .Bufu Tang et al. found that the expression level of RRM2 was positively correlated with neutrophil and macrophage infiltration in LUAD tissue, suggesting that RRM2 promoted lung cancer progression and affected macrophage infiltration, stimulated M1 phenotype polarization, and inhibited the M2 phenotype .
Our research has several limitations. First, our risk signature was built and validated using TCGA and the GEO database; no further external validation was performed with real-world prospective LUAD RNA-seq cohorts. Second, when analyzing the efficacy of immunotherapy, the LUAD sample size of the GSE126044 cohort was relatively small. Although positive results were obtained, their accuracy remains unclear. Therefore, we are attempting to find or establish a better and larger immunotherapy cohort of patients with LUAD. Third, although we identified a key role of RRM2 in LUAD using bioinformatics technology, we have not yet verified how RRM2 affects the development of LUAD by regulating mast cells and CD4 memory T cells or other mechanisms. Fourth, we screened drugs for high-risk patients with LUAD through the CMap and GDSC databases, but the specific efficacy and mechanism need to be further explored by cellular and animal experiments as well as clinical trials. Overall, more clinical samples and prospective experiments are necessary to study the mechanism by which RRM2 influences the LUAD immune microenvironment and immunotherapy.
We established a risk prediction signature based on FRGs and found that it was closely associated with the immune microenvironment of LUAD. The key risk gene RRM2 influenced the prognosis of LUAD by regulating the infiltration of activated mast cells and activated CD4 memory T cells. Moreover, FRG-based LUAD clustering was associated with the risk signature and immune cell infiltration. In other words, the FRG clustering we provided may indicate potential LUAD immune subtypes. Furthermore, a correlation analysis supported the relationship between the signature and LUAD immunotherapy. In addition, our results provide a new perspective that FRGs, including RRM2, influence anti-PDL1 immunotherapeutic responses by regulating immune checkpoint gene expression and infiltration of immune cells in the tumor microenvironment. Finally, our results suggest that RRM2 is an important biomarker affecting the long-term survival of patients with LUAD, and silencing RRM2 could promote the occurrence of ferroptosis in lung cancer cells. Taken together, our results offer novel insights into the research and treatment of LUAD.
Availability of data and materials
The data used to support the findings of this study are included within the article.
Mattiuzzi C, Lippi G. Cancer statistics: a comparison between World Health Organization (WHO) and Global Burden of Disease (GBD). Eur J Public Health. 2020;30(5):1026–7.
Couraud S, Zalcman G, Milleron B, et al. Lung cancer in never smokers—a review. Eur J Cancer. 2012;48(9):1299–311.
Devarakonda S, Morgensztern D, Govindan R. Genomic alterations in lung adenocarcinoma. Lancet Oncol. 2015;16(7):e342–51.
Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin. 2020;70(1):7–30.
Lennon FE, Cianci GC, Cipriani NA, et al. Lung cancer-a fractal viewpoint. Nat Rev Clin Oncol. 2015;12(11):664–75.
Siegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin. 2018;68(1):7–30.
Torti SV, Torti FM. Iron and Cancer: 2020 vision. Cancer Res. 2020;80(24):5435–48.
Wu Y, Zhang S, Gong X, et al. The epigenetic regulators and metabolic changes in ferroptosis-associated cancer progression. Mol Cancer. 2020;19(1):39.
Chen P, Wu Q, Feng J, et al. Erianin, a novel dibenzyl compound in Dendrobium extract, inhibits lung cancer cell growth and migration via calcium/calmodulin-dependent ferroptosis. Signal Transduct Target Ther. 2020;5(1):51.
Tang B, Xu W, Wang Y, et al. Identification of critical ferroptosis regulators in lung adenocarcinoma that RRM2 facilitates tumor immune infiltration by inhibiting ferroptotic death. Clin Immunol. 2021;232: 108872.
Huang F, Zheng Y, Li X, et al. Ferroptosis-related gene AKR1C1 predicts the prognosis of non-small cell lung cancer. Cancer Cell Int. 2021;21(1):567.
Wang M, Mao C, Ouyang L, et al. Long noncoding RNA LINC00336 inhibits ferroptosis in lung cancer by functioning as a competing endogenous RNA. Cell Death Differ. 2019;26(11):2329–43.
Mao C, Wang X, Liu Y, et al. A G3BP1-interacting lncRNA promotes ferroptosis and apoptosis in cancer via nuclear sequestration of p53. Cancer Res. 2018;78(13):3484–96.
Luo L, Wang H, Tian W, et al. Targeting ferroptosis for cancer therapy: iron metabolism and anticancer immunity. Am J Cancer Res. 2021;11(11):5508–25.
Ding L, Getz G, Wheeler DA, et al. Somatic mutations affect key pathways in lung adenocarcinoma. Nature. 2008;455(7216):1069–75.
Mou Y, Wang J, Wu J, et al. Ferroptosis, a new form of cell death: opportunities and challenges in cancer. J Hematol Oncol. 2019;12(1):34.
Hirsch FR, Scagliotti GV, Mulshine JL, et al. Lung cancer: current therapies and new targeted treatments. Lancet. 2017;389(10066):299–311.
Aird KM, Zhang G, Li H, et al. Suppression of nucleotide metabolism underlies the establishment and maintenance of oncogene-induced senescence. Cell Rep. 2013;3(4):1252–65.
Aird KM, Li H, Xin F, et al. Identification of ribonucleotide reductase M2 as a potential target for pro-senescence therapy in epithelial ovarian cancer. Cell Cycle. 2014;13(2):199–207.
Nordlund P, Reichard P. Ribonucleotide reductases. Annu Rev Biochem. 2006;75:681–706.
Rahman MA, Amin AR, Wang D, et al. RRM2 regulates Bcl-2 in head and neck and lung cancers: a potential target for cancer therapy. Clin Cancer Res. 2013;19(13):3416–28.
Jin CY, Du L, Nuerlan AH, et al. High expression of RRM2 as an independent predictive factor of poor prognosis in patients with lung adenocarcinoma. Aging. 2020;13(3):3518–35.
Altinkilic EM, Isbir S, Gormus U, et al. RRM1, RRM2 and ERCC2 gene polymorphisms in coronary artery disease. In Vivo. 2016;30(5):611–5.
Castelblanco E, Zafon C, Maravall J, et al. APLP2, RRM2, and PRC1: new putative markers for the differential diagnosis of thyroid follicular lesions. Thyroid. 2017;27(1):59–66.
Wang N, Zhan T, Ke T, et al. Increased expression of RRM2 by human papillomavirus E7 oncoprotein promotes angiogenesis in cervical cancer. Br J Cancer. 2014;110(4):1034–44.
Shao J, Zhou B, Chu B, Yen Y. Ribonucleotide reductase inhibitors and future drug design. Curr Cancer Drug Targets. 2006;6(5):409–31.
Wadler S, Makower D, Clairmont C, et al. Phase I and pharmacokinetic study of the ribonucleotide reductase inhibitor, 3-aminopyridine-2-carboxaldehyde thiosemicarbazone, administered by 96-hour intravenous continuous infusion. J Clin Oncol. 2004;22(9):1553–63.
Xia G, Wang H, Song Z, et al. Gambogic acid sensitizes gemcitabine efficacy in pancreatic cancer by reducing the expression of ribonucleotide reductase subunit-M2 (RRM2). J Exp Clin Cancer Res. 2017;36(1):107.
Cai X, Lin L, Zhang Q, et al. Bioinformatics analysis of the circRNA-miRNA-mRNA network for non-small cell lung cancer. J Int Med Res. 2020;48(6):300060520929167.
Jiang X, Li Y, Zhang N, et al. RRM2 silencing suppresses malignant phenotype and enhances radiosensitivity via activating cGAS/STING signaling pathway in lung adenocarcinoma. Cell Biosci. 2021;11(1):74.
We thank the Public Service Platform of South China Sea for R&D Marine Biomedicine Resources for support.
This project was supported by Basic and Applied Basic Research Program of Guangdong Province (2019A1515110201); The Science and technology program of Guangdong Province (2019B090905011); Discipline Construction Project of Guangdong Medical University (4SG21004G). Innovation and Entrepreneurship Team Leads the Pilot Program of Zhanjiang (2020LHJH005).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Ferroptosis-related clustering and seven FRGs affecting LUAD prognosis. (A, B, C, D) TCGA-LUAD cohort was divided into four clusters based on 70 differentially expressed FRGs. (E)Significant differences in OS between the 4 FRGs clusters. (F-L) Different expression levels of 7 key FRGs represent different OS in LUAD patients.
Enrichment analysis. (A) KEGG enrichment, fluid shear stress and atherosclerosis, HIF-1 signaling pathways, lipid and atherosclerosis, and ferroptosis are enriched. (B) GSEA enrichment, small cell pulmonary carcinoma, melanoma, cancer pathway, P53 signaling pathway, pathogenic Escherichia coli infection were significantly enriched in high-risk patients. (C) GO enrichment.
Clinical correlation heat map and the differential infiltration level of immunocytes in LUAD patients. (A)The proportion of four clinical stages in high- and low-risk groups. (B)Distribution of several clinical features between high- and low-risk groups. (C)The proportion of risk scores of four FRGs clusters. (D)The infiltration ratio of 22 immunocytes in each TCGA-LUAD sample. (E)Infiltration of 10 immunocytes in high- and low-risk groups. (F)The change trend of infiltration level of Mast cells resting with ALOX12B expression level.
The infiltration level of immune cells reflects the OS of LUAD. (A-L) Kaplan–Meier curve. Different infiltration levels of 12 immunocytes correspond to different overall survival rates in LUAD patients.
FRGs are involved in immune response and regulate the expression of immune checkpoints. (A)The change trend of TMB with RiskScore. (B)Different levels of TMB represent significantly different OS in LUAD patients. (C, D) Different ImmunScore and StromalScore represent significantly different OS in LUAD patients. (E) Mutation analysis of 7FRGs in TCGA-LUAD cohort. (F) The difference of RRM2 expression between LUAD patients with two different immunotherapy responses. (G)The difference of immunocytes infiltration and response to anti-PDL1 immunotherapy between high- and low-risk patients. (H)Relationship between PD1 expression and 15 prognostic FRGs. (I) Relationship between PDL1 expression and 15 prognostic FRGs.
Six antitumor drugs with significantly lower IC50 and six chemical structures of the top molecules for high-risk LUAD patients selected from CMAP database. (A)Bosutinib.(B)JNK.Inhibitor.VIII.(C)Docetaxel.(D)Tipifarnib.(E)Dasatinib.(F)Gefitinib.(G)Medrysone. (H)phenoxybenzamine. (I)Vorinostat. (J)Thioguanosine. (K)apigenin. (L)chrysin.
About this article
Cite this article
Deng, B., Xiang, J., Liang, Z. et al. Identification and validation of a ferroptosis-related gene to predict survival outcomes and the immune microenvironment in lung adenocarcinoma. Cancer Cell Int 22, 292 (2022). https://doi.org/10.1186/s12935-022-02699-4