A novel necroptosis-related gene signature in acute myeloid leukemia

ABSTRACT

Objective:

Necroptosis has been reported to play an important role in different cancers, including leukemia. However, biomarkers of necroptosis-related genes (NRGs) that help predict the prognosis of AML are still lacking. Our research aims to build a novel signature of NRGs that could enhance our understanding of the molecular heterogeneity in leukemia.

Method:

Gene expression profiles as well as clinical features were downloaded from TCGA and GEO databases. Data analysis were executed using R software version 4.2.1 and GraphPad Prism version 9.0.0.

Result:

Univariate cox regression and lasso regression were applied to identify survival-specific genes. Four genes including FADD, PLA2G4A, PYCARD and ZBP1 were considered as independent risk factors that affect the prognosis of patients. Risk scores were calculated according to the coefficient of four genes. Then clinical characteristics and risk scores were enrolled to construct a nomogram. CellMiner was also used to screen potential drugs and analyze the correlations between genes and drug sensitivity.

Conclusion:

In general, we established a signature of four genes related to necroptosis that could be helpful for future risk stratification in patients with AML.

KEYWORDS:

• Acute myeloid leukemia
• mRNA
• necroptosis
• prognostic signature

Introduction

Acute myeloid leukemia (AML) is an aggressive malignant disease with uncontrolled clonal expansion of myeloid blasts. The prognosis of AML is highly heterogeneous due to the molecular differences between patients [1]. Although some novel drugs and regimens have improved therapeutic efficacy [2], there are still many patients ultimately dying from primary or secondary resistance. Thus, establishing a prognostic model that includes molecular and genetic factors may help refine risk assessment and guide therapeutic decisions for AML patients.

Necroptosis is a newly identified type of programmed cell death that can be activated in the absence of apoptosis and triggered by the activation of death receptors [3]. Emerging studies indicate that receptor-interacting protein kinases 1 and 3 (RIPK1 and RIPK3) and downstream substrate pseudokinase mixed-lineage kinase domain-like (MLKL) play a crucial role in regulating necroptosis pathway [4]. The effects of necroptosis on cancer can be complex and paradoxical. On the one hand, many pivotal molecules in necroptotic signaling are downregulated in different types of cancer, such as breast cancer, colorectal cancer [5,6]. In this way, cancer cells may evade necroptosis and maintain proliferation. This suggests that necroptosis might help eliminate antiapoptotic tumor cells [7]. On the other hand, an increasing number of studies have also shown that necroptosis can promote tumor progression and metastasis by inducing an inflammatory response and creating a tumor-promoting microenvironment with elevated ROS production [8]. These changes in the tumor microenvironment can lead to genomic instability, ultimately accelerating malignant transformation [9,10].

The role of necroptosis in acute myeloid leukemia remains undefined. RIPK3 deletion transformed FLT3-ITD-driven myeloproliferation into AML in mice by increasing the accumulation of leukemia-initiating cells [11]. A-L Nugues et al. proved RIPK3 expression is significantly decreased in CD34+ blast cells from AML patients than CD34+ hematopoietic cells from healthy donors [12]. These studies demonstrated the antitumor role of the necroptosis pathway in AML. However, Xin et al found RIPK1/RIPK3 signaling was continuously activated in leukemia cells, suggesting that inhibiting RIP1/RIP3-mediated necroptotic signaling could be an alternative treatment strategy for AML patients [13].

In this study, we constructed a 4 NRGs prognosis model using the TCGA cohort and validated it on the GEO datasets. This model could efficiently predict the prognosis of AML patients based on necroptosis-related genes.

Materials and methods

Data collection

Gene expression profiles and corresponding clinical data for 151 AML patients were downloaded from the TCGA database (https://portal.gdc.cancer.gov). Patients with a survival time of 0 days were excluded. The GSE71014 and GSE12417 datasets were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/gds/). TCGA data were used as training set and GEO data as validation set. 159 necroptosis-related genes were collected from the Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/).

Construction and validation of the prognostic model

Univariate Cox regression analysis was performed to identify 38 candidate NRGs that were associated with OS in the TCGA LAML cohort. Genes with p-value <0.05 were further analyzed using Lasso regression (least absolute shrinkage and selector operation regression). Based on the regression coefficient (β), we calculate the risk score using the following equation: risk score = expression of gene 1* coefficient + expression of gene 2 * coefficient + expression of gene n * coefficient. Then, we classified AML patients into high-risk and low-risk groups based on the median risk scores. Kaplan–Meier survival analysis was used to compare the prognosis of these two groups. To evaluate the accuracy of this model, we calculated the area under the curve (AUC) at different time points. The calibration curve was used to assess the predictive discrimination of this risk model.

To test the predictive accuracy of risk model, we applied GSE71014 and GSE12417 as validation sets. Risk scores and AUC of each dataset were calculated in the same way as described above.

Predictive nomogram construction

Clinical parameters, including age, gender, percentage of blastcells, risk category, FAB category, leukocyte counts, hemoglobin, platelet counts and FLT3 mutation, NPM1 mutation, along with the risk score were gathered and analyzed by univariate and multivariate Cox proportional hazards regression. Due to the distinct clinical features of M3 subtype compared to other AML subtypes, patients with FAB classification were classified into M3 subtype and non-M3 subtype groups. Independent factors were selected to further construct a nomogram for predicting prognosis.

Enrichment analysis

We initially eliminated genes with low expression levels using a cut-off criterion for whose value was ‘0’ in over 50% of the samples. Subsequently, we identified differentially expressed genes (DEGs) with false discovery rate (FDR) < 0.05 and | log fold change (FC) | > 1 between the high-risk and low-risk groups. Furthermore, we conducted Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analyses with DEGs.

Quantitative reverse transcription real-time PCR

We collected bone marrow specimens from AML patients (including de novo and relapsed patients) as well as healthy donors. More details can be seen in Supplementary Table 1 and Table 2. Protocol of sample handling was performed in compliance with the Helsinki Declaration. Total RNA was extracted with TRIzol reagent and reversely transcribed into cDNA with reverse transcription reagent(Takara Primescript RT Master Mix). Then, real-time PCR was performed with SYBR Premix Ex Taq (Takara). The relative expression levels of four genes were normalized to GAPDH and then calculated by 2−ΔΔCt. The primers of four genes were shown in Supplementary Table 3.

Drug sensitivity analysis

Drug sensitivity data were downloaded from CellMiner website (https://discover.nci.nih.gov/cellminer/home.do). Target gene expression and cell sensitivity data of FDA-approved drugs were analyzed by Pearson correlation analysis.

Statistical analysis

Data analyses were performed using the R software version 4.2.1.

Cox regression and LASSO analysis were performed with the ‘survival’ and ‘glmnet’ packages. The number of folds used in cross-validation was 10 in lasso regression. The area under the ROC curve (AUC) was calculated using the packages ‘timeROC’ and ‘survivalROC’ for OS prediction.‘rms’package was utilized for nomogram and calibration diagrams. DEGs were identified using the ‘limma’ package. Function enrichment analysis containing GO and KEGG was performed by ‘clusterprofiler’ R package.

Data visualizing and statistical analysis of real-time PCR results were performed using GraphPad Prism software. T-test and Wilcoxon test were applied for data processing.

Statistically significant were considered when p-value was less than 0.05.

Results

Construction and validation of prognostic model in AML

In total, we gathered 159 necroptosis-related genes and included 129 AML samples into the training set. We first performed a univariate Cox regression analysis for preliminary screening and identified 38 NRGs (Table 1) with survival significance (p < 0.05).

Table 1. Overall survival-related NRGs in the AML patients (P < 0.05).

GENE	HR	CI95	P value
AIFM1	1.83	1.25–2.66	0.001735
BID	1.06	1.00–1.11	0.046672
CAMK2B	6.38	1.14–35.61	0.034648
CAMK2D	1.12	1.00–1.26	0.043262
CAPN1	1.04	1.02–1.06	0.000549
CHMP1A	1.12	1.03–1.22	0.009481
CHMP2A	1.05	1.02–1.08	0.002607
CHMP3	1.11	1.02–1.21	0.022011
CHMP4B	1.02	1.01–1.03	0.000734
CHMP7	1.1	1.00–1.21	0.044619
FADD	1.88	1.35–2.63	2.00E-04
FTH1	1.01	1.00–1.02	0.004134
FTL	1	1.00–1.00	0.034917
GLUD1	1.02	1.00–1.05	0.029998
GLUD2	2.31	1.32–4.04	0.00333
H2AC14	1.25	1.08–1.45	0.002369
H2AC4	1.35	1.00–1.83	0.049282
H2AJ	1.03	1.00–1.07	0.038691
H2AZ2	0.95	0.92–0.98	0.000555
IFNG	1.34	1.04–1.72	0.024342
IFNGR1	1.01	1.00–1.02	0.030796
MACROH2A1	0.98	0.97–0.99	0.003938
MACROH2A2	1.24	1.02–1.49	0.027693
PARP1	1.07	1.01–1.13	0.015206
PLA2G4A	1.08	1.04–1.12	7.00E-06
PYCARD	1.04	1.02–1.06	4.00E-06
RBCK1	1.06	1.02–1.10	0.007401
RIPK3	1.09	1.01–1.18	0.030392
SHARPIN	1.03	1.01–1.06	0.010026
SLC25A4	1.87	1.12–3.14	0.01736
SPATA2L	1.12	1.03–1.22	0.009326
STAT5B	0.96	0.94–0.99	0.004869
STAT6	1.02	1.01–1.03	0.003333
TICAM1	1.06	1.03–1.1	0.000544
TNF	0.93	0.88–0.99	0.013464
TRADD	1.2	1.08–1.32	0.00042
TYK2	1.02	1.00–1.04	0.021391
ZBP1	1.54	1.24–1.91	0.000101

The lasso regression analysis was further performed to find out the crucial survival-related NRGs (Figure 1). Four genes were finally chosen to establish an optimized prognostic signature. The risk score of this model was calculated as follows: $\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=0.4708\ast \mathrm{F}\mathrm{A}\mathrm{D}\mathrm{D}+0.088\ast \mathrm{P}\mathrm{L}\mathrm{A}2\mathrm{G}4\mathrm{A}+0.0265\ast \mathrm{P}\mathrm{Y}\mathrm{C}\mathrm{A}\mathrm{R}\mathrm{D}+0.3811\ast \mathrm{Z}\mathrm{B}\mathrm{P}1$In total, 129 patients with AML were evenly categorized into low-risk and high-risk groups according to the risk score. Our analysis revealed a positive correlation between risk score and mortality rate, as demonstrated in Figure 2(A,B). Additionally, the risk heatmap indicated differential gene expression between the two groups, as depicted in Figure 2(C). The Kaplan-Meier survival curve further confirmed that patients in the high-risk group had significantly poorer overall survival rates compared to those in the low-risk group (Figure 2(D)).

Figure 1. LASSO Cox regression performed to identify four necroptosis-related genes associated with OS.

Figure 2. (A,B) Risk scores and survival status of AML patients. (C) Heatmap depicting the different expressions of the four NRGs in different AML cohorts. (D) Kaplan–Meier survival curve for the four gene signature.

The expression profile of the four genes and their correlation with survival were examined in the TCGA LAML dataset using the GEPIA2 platform (gepia2.cancer-pku.cn) (Figure 3).

Figure 3. The survival curve of each gene from GEPIA2 platform.

ROC curve of the OS was plotted to reveal the predictive performance of the risk score. Notably, the area under the curve (AUC) value of the risk score was significantly higher than that of other clinical parameters, including age, gender, blastcell percentage, FAB classification, leukocyte counts, hemoglobin, platelet counts, risk category, FLT3 mutation and NPM1 mutation (Figure 4(A)). The AUC value of the risk score for the 1-, 2-, 3-year OS were 0.8,0.81 and 0.8 respectively (Figure 4(B)). Additionally, the calibration line for risk score showed desirable prediction accuracy at 1, 2 and 3 years (Figure 4(C)). To further substantiate the accuracy of our prognostic model, the GSE71014 and GSE12417 datasets were enrolled as validation sets. We calculate the risk score and AUC of each dataset (Figure 4(D,E)).

Figure 4. (A) ROC analysis of OS for risk score and clinical parameters. (B) Time-dependent ROC curves of training set for 1-, 2-, and 3-year overall survival. (C) Calibration curve for the four gene signature at 1, 2, 3 years. Time-dependent ROC curves of GSE71014 (D) and GSE12417 (E) for 1-, 2-, and 3-year overall survival.

In the training cohort, we gathered clinical data on patients and conducted univariate and multivariate regression analysis. The results of the univariate regression revealed that survival time was associated with risk category, risk score, age and FAB classification (Figure 5(A)). Subsequent multivariate analysis showed that age and risk score were independent risk factors (Figure 5(B)).

Figure 5. Forest plots of the univariate (A) and multivariate (B) Cox regression analysis in AML.

The Wilcoxon test was performed to explore the associations between the risk score and several clinical characteristics. Age, risk category, and NPM1 mutation were found to be related to the risk score (Figure 6).

Figure 6. The correlation between risk score and clinical characteristics (age, gender, leukocyte, risk category, FLT3 mutation and NPMc mutation).

Based on the regression results above and our clinical experience, we constructed a nomogram to help make a more accurate prognostic analysis (Figure 7(A)). Next, the time-dependent ROC of the nomogram was calculated, and the AUC values of 1-, 2- and 3-year OS were 0.81, 0.84, and 0.87, respectively (Figure 7(B)).

Figure 7. (A) Nomogram predicting 1-, 2- and 3-year OS for the AML patients. Leukocyte: <20, counts less than 20 × 10⁹/L; ≥20, counts equal to or more than 20 × 10⁹/L, risk category: 0, Favorable; 1, Intermediate; 2, Poor. (B) Time-dependent ROC curves of the nomogram points for 1-, 2- and 3-year overall survival.

DEGs, GO and KEGG

To further analyze the difference between the high-risk and low-risk subgroups, ‘limma’ packages were applied to screen the differentially expressed genes (DEGs) with false-discovery rate (FDR) < 0.05 and | log fold change (FC) | > 1. A total of 1880 DEGs were identified, with 829 upregulated and 1051 downregulated (Figure 8(A)). The GO and the KEGG pathway enrichment analysis of 1880 DEGs were summarized in Figure 8. In the GO enrichment analysis, the top enriched terms in biological processes, cell component and molecular function were pattern specification process, collagen-containing extracellular matrix and receptor-ligand activity respectively (Figure 8(B)). In the KEGG enrichment analysis, the top enriched term was cytokine-cytokine receptor interaction (Figure 8(C)).

Figure 8. Functional assays of DEGs between low and high risk groups. (A) The volcano plot of the differentially expressed genes. (B) The top 10 GO terms of BP, CC and MF. (C) The top 10 KEGG pathways.

Measuring gene expressions in AML patients by real-time PCR

Real-Time PCR was performed to verify the mRNA expression levels of 4 genes in different cohorts. FADD, PLA2G4A, PYCARD and ZBP1 were downregulated in de novo AML patients compared to healthy donors (Figure 9(A)). However, gene expression increased in the relapsed cohort (Figure 9(B)). The results gave us a hint that these genes might be associated with disease recurrence. Larger sample sizes and more experiments are needed to verify the presumption.

Figure 9. The mRNA expression levels of 4 genes (FADD, PLA2G4A, PYCARD and ZBP1) in different donors. (A) de novo patients vs healthy donors and (B) de novo patients vs relapsed patients. Error bars are indicative of ±SEM, and threshold for significance was P < 0.05. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Drug sensitivity analysis

Potential correlations between the expression of four genes and drug sensitivity in different cancer cell lines were analyzed from CellMiner database. Results with P < 0.05 and |Cor|>0.3 were considered statistically significant. The top16 most relevant correlations were shown in Figure 10. Positive correlations were found between 5 drugs (LDK-378, Alectinib, timazid, rifa and ciclosporin) and ZBP1. Similar results were identified for the relationship between PLA2G4A expression and Megestrol acetate, Isotretinoin, LOXO-101, Sunitinib. PYCARD expression was also positively associated withdrug sensitivity to Artemether, Cyclophosphamide, Hydroxyurea, ABT-199, and Parthenolide. In addition, there were negative correlations between two drugs (BLU-667, Sonidegib) and PYCARD.

Figure 10. Top16 most relevant correlations between the expression of 4 genes and drug sensitivity in different cancer cell lines from CellMiner database.

Discussion

Necroptosis is a caspase-independent form of programmed cell death mediated primarily by receptor-interacting protein kinase 1 and 3 (RIP1 and RIP3), as well as mixed lineage kinase domain-like protein (MLKL) [14]. Emerging evidence suggests that necroptosis is involved in regulating different tumors. However, the exact role of necroptosis in neoplasms, including hematologic malignancies like AML, remains unknown. Ulrike Hockendorf demonstrated that RIPK3 signaling was crucial in the mechanism of leukemia suppression and exogenous activation of necroptosis could induce LIC cells and prompt differentiation in AML [15]. On the contrary, Xin et al. showed that RIPK3 was highly expressed in the M4 and M5 subtypes of AML and illustrated that inhibition of the necroptotic pathway could be beneficial for certain patients with AML [13]. Consequently, necroptosis is likely to be a double-edged sword in AML and requires further study.

In this study, 38 NRGs were found to be significantly associated with OS in patients with AML using univariate Cox regression. Further analysis by Lasso regression led to the selection of 4 NRGs to build a prognostic model. These genes have been reported to have an essential role in the development of various cancers.

The FADD gene is located on chromosome 11q13.3 and encodes a classical apoptotic signaling adapter-FADD protein. Amplification of FADD has been reported in many malignant tumors, such as breast cancer, lung cancer, skin cancer, and esophageal cancer [16]. In some cancers, patients with FADD overexpression are prone to develop metastases, indicating a worse prognosis and poor overall survival [17].

Phospholipase A2 (PLA2s) are key enzymes that catalyze the hydrolysis of membrane phospholipids at the position sn-2 to release second messenger lipids. There are three main classes of PLA2, including secretory (sPLA2), Ca2 + independent (iPLA2), and cytosolic PLA2 (cPLA2) [18]. PLA2G4A belongs to cPLA2 which plays a fundamental role in inflammation and neoplasm. For example, researchers identified the growth dependency of AML cells on PLA2G4A by knocking out candidate genes with CRISPR/CAS9 [19].

PYCARD encodes an adaptor protein that is composed of two protein–protein interaction domains: an N-terminal PYRIN-PAAD-DAPIN domain (PYD) and a C-terminal caspase-recruitment domain (CARD). The PYCARD protein was originally identified in HL-60 by forming large aggregates called speck cells after inducing apoptosis with retinoid acid and other drugs [20]. Deswaerte et al. discovered that PYCARD was significantly upregulated in gastric cancer compared to normal gastric tissue [21].

ZBP1 was previously identified as an IFN-inducible, tumor-associated protein [22]. Baik JY et al. demonstrated that ZBP1 is a key mediator of tumor necroptosis, and that deletion of ZBP1 inhibits tumor metastasis in the MVT-1 breast cancer model [23].

Our research presents a novel predictive model for acute myeloid leukemia (AML) that may enhance risk stratification beyond the current hierarchical prognostic models based on limited number of genes as outlined by the NCCN guidelines. Our model offers the potential to address the limitations and deficiencies of existing approaches to predict patient survival in AML.

While other studies have also developed NRGs models for AML [24,25], we took a different approach to identify prognostic-related genes. Specifically, we avoided obtaining white blood cell cohort from other databases due to potential bias and tendencies that may arise from differences between data sources, given that TCGA does not provide normal blood cells data. Several limitations must be acknowledged in our study. Firstly, although the WHO AML classification is widely used currently, the data from TCGA are based on the FAB classification. Secondly, the sample sizes of the training and validation sets were far from sufficient, indicating that larger datasets are needed to validate the performance and reliability of our predictive model. Thirdly, additional biological experiments and clinical trials are necessary to verify the accuracy and effectiveness of our nomogram.

In our study, we developed a new prognostic signature utilizing bioinformatics analysis to predict the prognosis of AML patients based on necroptosis-related genes. Although our model has shown promising results, further research is required to fully elucidate the underlying mechanisms of necroptosis in AML.

Ethical statement

The study was approved by the ethics committee of the participating institutions. Written informed consent was obtained from the patients.

Authors' contribution

All authors contributed to the study conception and design. Weiyue Fang, Hongdou Lin and Junyi Chen collected and analyzed the data and wrote the manuscript. The first draft of the manuscript was written by Weiyue Fang. Wenjian Guo provided advice and reviewed the manuscript. All authors read and approved the final manuscript and agreed to the submission.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The datasets generated and/or analyzed during the current study are available in the online database TCGA (https://portal.gdc.cancer.gov/), GEO (http://www.ncbi.nlm.nih.gov/geo/) and CellMiner website (https://discover.nci.nih.gov/cellminer/home.do).

Additional information

Funding

This work is supported by the Medical Health Science and Technology Project of Zhejiang Provincial Health Commission [grant number KYYW202029] and Project of Wenzhou Science and Technology Bureau [grant number 2022ZY0005].