Identification of a venetoclax-resistance prognostic signature base on 6-senescence genes and its clinical significance for acute myeloid leukemia

Background Satisfactory responses can be obtained for acute myeloid leukemia (AML) treated by Venetoclax (VEN)-based therapy. However, there are still quite a few AML patients (AMLs) resistant to VEN, and it is critical to understand whether VEN-resistance is regulated by senescence. Methods Here, we established and validated a signature for predicting AML prognosis based on VEN resistance-related senescence genes (VRSGs). In this study, 51 senescence genes were identified with VEN-resistance in AML. Using LASSO algorithms and multiple AML cohorts, a VEN-resistance senescence prognostic model (VRSP-M) was developed and validated based on 6-senescence genes. Results According to the median score of the signature, AMLs were classified into two subtypes. A worse prognosis and more adverse features occurred in the high-risk subtype, including older patients, non-de novo AML, poor cytogenetics, adverse risk of European LeukemiaNet (ELN) 2017 recommendation, and TP53 mutation. Patients in the high-risk subtype were mainly involved in monocyte differentiation, senescence, NADPH oxidases, and PD1 signaling pathway. The model’s risk score was significantly associated with VEN-resistance, immune features, and immunotherapy response in AML. In vitro, the IC50 values of ABT-199 (VEN) rose progressively with increasing expression of G6PD and BAG3 in AML cell lines. Conclusions The 6-senescence genes prognostic model has significant meaning for the prediction of VEN-resistance, guiding personalized molecularly targeted therapies, and improving AML prognosis.


Introduction
Acute myeloid leukemia (AML) is one of the most common hematological malignant cancers, which is far more common in elderly patients (1,2).Traditional chemotherapy era, elderly AML patients (AMLs) have a much poorer prognosis, with a 5-year survival rate of only 5% after the diagnosis (3).With the recent advent of molecularly targeted therapies, such as B-cell lymphoma 2 (BCL-2) inhibitor, survival of older AMLs has been improved (4,5).
The BCL-2 protein is a key regulator of the mitochondrial apoptotic pathway and plays an important role in the survival and persistence of AML blasts (6,7).Targeting BCL-2, Venetoclax (VEN) showed an efficient strategy to promote caspase-dependent cell death in AML (4,8).In accordance with these studies, VEN has been approved for the treatment of newly-diagnosed elderly AMLs.VEN-based therapy can induce approximately 70% therapeutic responses in older AMLs.However, a significant minority of AMLs lack therapeutic response to initial induction or reinduction of VEN Monotherapy (9).The short duration of response and development of resistance have become major concerns.Previous studies have found that key contributing factors to VEN resistance include dependencies on alternative anti-apoptotic BCL-2 family proteins, selection of the activating kinase mutations, TP53 mutation, and BAX variants (10)(11)(12)(13)(14).More research is needed to explore the mechanisms of VEN resistance in AML and try to find strategies to overcome the resistance.
Senescence is the natural consequence of telomere shortening at the chromosome ends upon extensive replication, but it can also be induced by DNA damage and imbalances in cellular signaling networks (15,16).Cellular senescence response may suppress cancer progression in vivo (17-19), but could also variously stimulate tumor progression in some conditions, as well as associated with various age-related diseases (20)(21)(22).The elimination of senescent cells can delay multiple age-related symptoms, and reduce incidences of spontaneous tumorigenesis and cancer-related mortality (23).Therefore, tumor cells can undergo senescence as an evolutionary process, including both tumor-intrinsic characteristics and extrinsic immune pressure (24,25).However, a comprehensive understanding of the influences of senescence on VEN-resistance in AML is still lacking.In this study, we developed a VEN-resistance senescence prognostic model (VRSP-M) across multiple AML cohorts.The prognostic signature has significant meaning for the prediction of VENresistance, guiding personalized molecularly targeted therapies, and improving AML prognosis.

Data source
The profiles of TCGA AML (n=151) were downloaded from the website of UCSC Xena (https://gdc.xenahubs.net)and exploited to build a prognostic signature for AML based on VRSGs (VEN resistance-related senescence genes).Human-related senescence genes (HRSGs, n=279, Supplementary Table S1) were acquired from the HAGR website (https://genomics.senescence.info/cells/).Ex vivo data from Beat AML cohort (26) was used to identify VRSGs and a total of 343 AMLs were enrolled to validate the relationship between prognostic signature and clinical manifestation.To test the applicability, we further verified the effects of a predictive model in non-APL (acute promyelocytic leukemia) AML (GSE106291, n=250) and normal karyotype AML (GSE71014, n=104).The inclusion criteria of Beat AML contained: expression profiles at time of the initial diagnosis, complete data of survival and ELN stratification, but excluding duplicated cases.For other three datasets, all patients with survival and expressed information were included in this study.All the data used in this study were obtained from the public program, and all processes complied with the publication guidelines.Therefore, ethical approval of local ethics committees is not required.

Identification of VRSGs
From the drug response of beat AML in Ex vivo, samples with the lowest 20% of area-under-the-curve values (AUCs) were deemed to be sensitive to VEN, while those with the highest 20% AUCs were considered as VEN-resistance.A total of 3023 differential expression genes (DEGs) were identified between VEN-resistant and -sensitive samples through the DEseq2 method (|log2FC| ≥1.0 and adjusted P value < 0.05, Supplementary Table S2).The intersections of DEGs with HRSGs were identified as VRSGs (n=51).

Prognostic model generated from VRSGs
According to the median expressed levels of VRSGs and univariate Cox analysis, 18 of 51 VRSGs were proved to be associated with AML prognosis in the modeling set (Supplementary Table S3, P<0.05).Of 18 VRSGs, the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm was performed to screen the optimal senescence genes to develop VRSP-M in TCGA-AML through "glmnet" R packet.A 10fold cross-validation method was employed to hold stability, and the minimum criteria was chosen as the optimal penalty value (l) (set.seed,2021).The risk score was calculated using the expression of model's genes as follows: To evaluate and validate VRSP-M in the modeling and validation datasets, AMLs could be divided into high-and lowrisk subtypes according to the median score.Survival analyses were applied to distinguish the difference between these two groups by "survival" package.Curves of receiver operating characteristic (ROC) and AUCs were used to assess the accuracy of the prognostic model using the "timeROC" package.Then Beat AML and GSE106291 were combined after removing batch effects by the "SVA" package.A subset containing 80% samples of the combined dataset was re-sampled 100 times and used to examine the robustness of the VRSP-M for predicting OS of AML patients.

Functional analyses
A previous study demonstrated that a monocytic clone of AML could confer VEN-resistance (27), so its markers were obtained to assess the enrichment of monocyte differentiation.Moreover, the SenMayo set (28) was also used to estimate the degree of enrichment in the senescence pathway.Therefore, four gene sets were acquired to inquire into the biological function of VRSP-M, including c2.cp.kegg.v7.5.1.symbols.gmt,c2.cp.reactome.v7.5.1.symbols.gmt,monocyte differentiation (Supplementary Table S4), and the SenMayo gene set (Supplementary Table S5).A value of false discovery rate (FDR) < 0.05, adjusted P < 0.05, and |NES| (normalized enrichment score) >1.5 were considered as significant enrichment in gene set enrichment analysis (GSEA).

Components analysis of immune cells
Using the "GSVA" and "GSEABase" R packages, the different components of immune cells between the high-and low-risk subtypes were reckoned and compared using single-sample gene set enrichment analysis (ssGSEA).The gene set was collected from the previous study (29), containing 28 types of immune cells.Then xCell (30) and ESTIMATE (31) algorithms were further utilized to impute the weights of M2 macrophages, immune and stromal score, respectively.

Predicted response of immunotherapy
The correlation with Spearman method was performed to check the link between the risk score of VRSP-M and eight immune checkpoints (SIGLEC15, TIGIT, CD274, HAVCR2, PDCD1, CTLA4, LAG3, and PDCD1LG2).We also predicted whether high-risk subtype could benefit from blockade therapy of immune checkpoints using the algorithm of Tumor Immune Dysfunction and Exclusion (TIDE) (32).

Screening marker genes of VEN-resistance
Using one dataset of CRISPR-Cas9 screens (GSE216087) (33), we checked whether model's VRSGs dysregulated in AML cells after single-VEN treatment.When sgRNAs depleted significantly on OCI-AML2 cells after VEN treatment, it meant that knockdown of these genes increased sensitivity to VEN, and high expression of these genes could decrease sensitivity to VEN and contribute to resistance.
To identify which model's VRSGs could be an effective biomarker of VEN-resistance, we first applied ABT-199 (venetoclax) to search molecular biomarkers of treatment response to VEN using a software of computational analysis of resistance (CARE) (34), the score of which indicates the correlation between gene alteration and drug efficacy.A higher positive score suggests better drug response, while a negative score demonstrates drug resistance.Furthermore, we also got data on CRISPR loss-of-function from previous VEN-resistance research (13).The negatively selected gene indicated that higher expression confers resistance to VEN.Then we verified in vitro.

Cell lines and RT-qPCR
AML cell lines, including HL-60, MOLM13, MV4-11, THP-1 OCI-AML3, and K562, were purchased from the American Type Culture Collection (Manassas, United States).All of them were cultured in RPMI 1640 with 10% FBS (Gibco, United States) as well as antibiotics (1% penicillin-streptomycin).All cells were kept at 37°C in an incubator with 5% CO2.The total RNA was extracted using RNA Isolator Total RNA Extraction Reagent (R401-01, Vazyme, China), and reverse transcribed to cDNAs using the PrimeScript ™ RT Master Kit (RR036A, Takara, Japan).Quantitative real-time PCR (qRT-PCR) was performed using 2x SYBR Green qPCR Master Mix (B21202, Bimake, United States) with 7500 real-time PCR system (Applied Biosystems, United States).The sequences of gene-specific primers are summarized in Supplementary Table S6.Gene expression levels were quantified with the 2-DCt method and GAPDH was used as endogenous control.

Western blot
Protein preparation and western blot assay were performed as described previously (35).BAG3 Ab (sc-136467, United States) was purchased from Santa Cruz Biotechnology, and GAPDH Ab (ab8245, United States) was obtained from Abcam.G6PD Ab (AF6945, China) and secondary Abs, such as HRP goat anti-mouse IgG (A0216, China) and goat anti-rabbit IgG (A0208, China), were bought from Beyotime Biotechnology.

Cell viability assay
The cytotoxic effects of ABT-199 on AML cell lines were determined by a Cell Counting Kit-8 (CCK-8; B34304, Bimake, United States) assay.ABT-199 was diluted in 100 μl of growth medium to designated doses, and leukemia cells were added to the 96-well plate (1×10 4 cells per well in 100 μl).Cultured leukemia cells were incubated in the presence of the drug for 48 hours at 37°C in a humidified 5% CO2-95% air incubator.Then, 10 mL of the CCK-8 reagent was added into each well, and optical densities at the wavelength of 450 nm were measured using the Synergy HTX Multimode Reader (BioTek, United States).The percentage of surviving cells was calculated according to the absorbance ratio of the test well and control well.IC50 values were calculated and visualized using GraphPad Prism software v9.0 (GraphPad, La Jolla, CA) based on the percent of live cells.

Statistical analysis
Categorical variables were compared by Chi-square or Fisher's exact test, and continuous variables were checked using the non-parametric test.Survival differences were quantified through the log-rank method or Cox regression analysis, and visualized using Kaplan-Meier curves.All analyses were conducted using the statistical software R (version 4.1) (http://www.R-project.org).The P-value of two-sided was set at the 0.05 significance level.

Construction of VEN-resistance senescence prognostic model
The flow chart is presented in Figure 1.A total of 3023 DEGs were found differently between VEN-sensitivity and -resistance (Figure 2A).Fifty-one genes were identified as VRSGs (Figure 2B), and 18 of which were verified with AML prognosis (Figure 2C).Through the LASSO method, 6 of 18 prognosis-related VRSGs were selected to construct VRSP-M to predict overall survival (OS) in TCGA AML (Figures 3A,  B).The risk score of each patient was determined based on the following formula (EXP indicated the expression of each gene): The distribution of risk score, survival status, and expression of 6-VRSGs are shown in Figure 3C.Associated with the elevated risk scores, mortality risk increased, while the survival time decreased.AUCs at 1, 3, and 5 years were 0.738, 0.736, and 0.864 (Figure 3D), which indicates that VRSP-M had better accuracy than random choice.Moreover, AMLs with high-risk scores presented a worse OS (P < 0.001, Figure 3E).

VRSP-M correlates with adverse features and monocyte differentiation in AMLs
Table 1 summarizes the clinical characteristics in high-and lowrisk subtypes in TCGA-AML.AMLs with high risk were older than those in the low-risk subtype, with median ages of 62 and 51 years respectively (Figure 4A, P < 0.001).More AMLs with intermediate or poor cytogenetic risk were located in the high-risk subtype (Figure 4B, P<0.001).In the distribution of the French-American-British (FAB) subtype, M5 possessed the highest risk score, while M3 had the lowest score (Figure 4C).Furthermore, GSEA analysis also confirmed the enrichment of monocyte differentiation for the high-risk subtype (Figure 4D, NES = 2.475).
The different expressions of key genes in apoptosis pathways may help to explain why the high-risk subtype was more likely to be resistant to VEN (Figure 4E).The BCL2 expression in the high-risk subtype was significantly lower than that in the low-risk subtype, but opposite results occurred in the expression of BCL2L1, BCL2A1, and MCL1.In addition, although the high-risk subtype possessed a high expression of pro-apoptotic genes (BID, BMF, BCL2L11) and apoptosis receptor (BAX, BAK1), the inhibitor of apoptotic genes (BIRC3, NAIP) was also higher in high-risk subtype than those in the low-risk group (adjusted P < 0.05).

Validating the clinical significance of VRSP-M in AML patients
Three AML cohorts were used to verify the performance of 6genes VRSP-M.In the Beat AML cohort, the clinical features are depicted in Table 2.A higher proportion of M5 subtype, non-de novo AML, and adverse risk of ELN recommendation was presented in AMLs with high-risk scores (P < 0.05).When taking Frontiers in Oncology frontiersin.orggenetic mutations into account, patients with no mutation data were excluded from the analysis.AMLs with high-risk scores were frequently accompanied by mutations of DNMT3A (P = 0.005), NRAS (P = 0.007), KARS (P = 0.066), and TP53 (P = 0.040), while the proportion of FLT3-ITD mutation in low-risk subtype was significantly higher than that in high-risk subtype (P = 0.002).
Consistent with the results in the modeling set, as risk scores increased, death risk also went up in three validation cohorts (Figures 5A-C).A worse OS was confirmed for the high-risk subtype (Figures 5D-F), and the VRSP-M also had the ability to predict AML prognosis (Figures 5G-I).Similarly, patients with high-risk scores were mainly involved in monocyte differentiation (Figures 5J-L).Leukemia heterogeneity may limit the practicability of prognostic signatures.Thus, we combined Beat AML and GSE106291 after removing batch effects by the SVA package.To check the robustness of this 6-genes prognostic signature, a total of 100 re-sampling tests were conducted randomly in 80% samples of the combined dataset.The results indicated that the P values were less than 0.05 in each Kaplan-Meier and univariate Cox analysis (Table S7), indicating a high performance for OS prediction in AML.

Functional signaling pathways
GSEA analyses were performed to better investigate the potential function of VRSP-M.The high-risk subtype was highly enriched in the senescence pathway (Figure 6A).KEGG results revealed that high-risk subtype enriched in Lysosome, Hematopoietic cell lineage, Cell adhesion molecules, and many immune-related pathways, such as cytotoxicity mediated by natural killer cell, interaction of cytokine and cytokine receptor, chemokine and T cell receptor signaling pathway (Figure 6B).Moreover, the main items of Reactome analysis were interactions between lymphoid and non-lymphoid cell, NADPH oxidases, interleukin 10 and PD1 signaling (Figure 6C).

VRSP-M correlates with immune features in AMLs
Many enriched items of immune-related pathways prompted us to explore the immune features associated with VRSP-M.We found that patients with high-risk score were often accompanied with a higher immune cell infiltration, including various subtypes of     immune activated and immunosuppressive cells (Figure 7A).The high-risk AMLs both had a higher immune and stromal score (Figure 7B).In addition, as well as Myeloid-derived suppressor cells (MDSCs) (Figure 7A), a higher proportion of M2 macrophage cells was also observed in patients with high-risk score (Figure 7C).

VRSP-M associated with immunotherapy response of AMLs
On account of the above discoveries, we speculated that PD1 plays a vital role in VRSP-M.The results showed that the risk score of the VRSP-M was positively correlated with CTLA4 (R=0.402),PDCD1 (R=0.398),HAVCR2 (R=0.328),PDCD1LG2 (R=0.283),CD274 (R=0.281),LAG3 (R=0.324), and T cell dysfunction (R=0.551),but negatively related with TIDE score (R= -0.489) and T cell exclusion (R= -0.491) (Figures 8A-D, P < 0.001).With a lower TIDE score (Figure 8E, P < 0.001), the high-risk subtype had a higher number of responders from immunotherapy (53.3% vs. 21.1%,P < 0.001) (Figure 8F).According to the expression of PDCD1 and VRSP-M's risk score, AMLs were divided into four groups.Patients both with a high level of risk score and PDCD1 presented a worse prognosis (Figure 8G, P < 0.001).These results indicated that patients with high-risk scores may benefit from the blocking therapy of immune checkpoints.

Verification of the maker genes of VRSP-M in vitro
In this signature, five senescence genes (except CDK6) were positively associated with higher VRSP-M scores.We then checked whether these five genes were dysregulated in AML cells using another dataset of CRISPR-Cas9 screens (GSE216087).The results indicated that four sgRNAs (except GRK6) were depleted significantly on OCI-AML2 cells after single-VEN treatment (Figure S1), suggesting that overexpression of these genes could confer VEN-resistance.
Whereafter, a method of CCK-8 assay was used to examine the viability of AML cell lines treated with ABT-199 for 48 hours and determined the IC50 values.The IC50 values ranged from <10 nmol/L to >1000 nmol/L (Figure 10A).The expression levels of G6PD and BAG3 were further tested in these AML cell lines, such as HL-60, MOLM13, MV4-11, and so on.The results indicated that G6PD expression was positively correlated with the IC50 of ABT-199 (Figure 10B).Except for OCI-AML3, the expression of BAG3 was also positively correlated with VEN-resistance (Figure 10C).The expression of G6PD and BAG3 were also verified in these AML cell lines using Western blot.As shown in Figures 10D, E, the IC50 levels of ABT-199 rose progressively with increasing protein levels of G6PD and BAG3.These findings indicate that G6PD and BAG3 may be effective markers for VEN-resistance in AML.

Discussion
Senescence is a complex stress response that can be grouped into different categories including genome-based failures and signaling dysfunction.However, the role of cellular senescence in cancer is controversial.In some conditions, the response of cellular senescence suppresses cancer progression (17, 18), conversely, which variously stimulates tumor progression in other ways (20)(21)(22).However, it is not quite clear whether senescence could induce VEN-resistance in AML.To better understand it, a VRSP-M was developed and validated using multiple AML cohorts, which can distinguish the prognosis of AMLs.
In this prognostic model, G6PD, BAG3, SRC, TNFSF15, and GRK6 act as risk factors, whereas CDK6 is a protective factor.Moreover, G6PD possesses the highest weight on AML prognosis, and was also proven to be an effective molecular marker of VENresistance.Previous studies have indicated that G6PD overexpression was associated with a poor prognosis in certain types of cancer, including AML, hepatocellular carcinoma, invasive breast carcinoma, and mesothelioma (36).G6PD could promote cancer progression through its effects on some metabolic pathways (37,38).Decreasing proliferation of leukemia and other cancer cells, knockdown of G6PD significantly increased apoptosis of tumor cells which are also more susceptible to oxidative stress (39,40).Acting as an effector of ATM (Ataxia telangiectasia mutated), G6PD often participates in the development of various cancers through metabolic programming and DNA repair pathways (41).As an essential enzyme in the pentose phosphate pathway (PPP), G6PD could produce more materials by this pathway to meet the high anabolic needs of tumor cells, which may make the cancers more resistant to chemotherapy.Playing an important role in AML resistance to the FLT3 inhibitor, the inactivation of G6PD increases the sensitivity of AML to FLT3 inhibitors (42).Under various stresses, multiple tumors turn metabolism to the PPP to get enough reductants to fight against reactive oxygen species (ROS) by activating G6PD rapidly.Otherwise, G6PD can also affect cancers by regulating ROS.Modulation of G6PD was proven to affect bladder cancer via ROS accumulation and the Validation of VRSP-M in three independent AML cohorts.AKT pathway in vitro (43), and knockdown of G6PD reduces ROS accumulation and enhances apoptosis of bladder cancer cells.In addition, G6PD also facilitates clear cell renal cell carcinoma invasion through the ROS−MAPK axis pathway (44).
Up-regulated ROS level induced by G6PD activation not only leads to out of control of cell growth and apoptosis of cancer, by also affects the immune microenvironment.In our study, high-risk AMLs also activate the signal of NADPH oxidases, and excessive activation of NOX (non-phagocytic cell oxidase) in cancer cells results in a large amount of ROS production.Accumulation of ROS leads to apoptosis in normal cells and mediates cellular senescence (45), which may be one of the reasons for the activation of NADPH oxidases in the highrisk subtype.Otherwise, excessive ROS also maintains the proliferation of tumor cells due to the protection of anti-oxidative stress with NADPH (39,(46)(47)(48).In fact, overexpression of G6PD could protect leukemic cells against oxidative stress by increasing NADPH production.On the other hand, not only leading to DNA damage and genomic instability, excessive of ROS could regulate signal transduction in the tumor environment, all of these are beneficial to the growth and progress of tumors (49, 50).Previous study has demonstrated that elevated ROS level induced by VEN can enhance the anti-leukemia effect of T cells (51).However, many studies also confirmed that ROS could induce the polarization of macrophages to M2 subtypes (52-56).Acting as a double-edged sword, ROS not only enhances T cells' anti-leukemia effect, but stimulates other immune cells, such as M2 macrophages and MDSCs.Hence, cellular senescence often triggers an immune response in the tumor microenvironment, facilitating tumor formation and progression (57,58), which might lead to a higher immune and stromal score in our high-risk AMLs.Moreover, a higher infiltration of MDSCs and M2 Macrophages in the high-risk subtype could prevent immune clearance of leukemia cells, and lead to poorer prognosis (59).
Shaping an unfavorable immune microenvironment, immunosenescence is also an urgent problem to be solved in the treatment of cancers.Mainly involved in PD1 signaling and strongly related to many immune checkpoints, a high-risk score of the VRSP-M was positively correlated with T cell dysfunction but negatively with T cell exclusion.T-cell dysfunction in cancer displays functional unresponsiveness, including senescence, exhaustion, anergy, and self-tolerance that is increasingly recognized as major hurdles for the success of cancer immunotherapy (60)(61)(62).So, potential approach to enhance antileukemia is to improve T and NK dysfunction, such as PD-1 inhibitor, chimeric antigen receptor T-cell therapy, and NK or gd Tbased adoptive immunotherapies (63-65).While the blocking-up of immune checkpoints has led to breakthroughs in several solid cancer therapies, research in AML remains limited (66).In the real world, AMLs have limited benefits from anti-PD-1 therapy (67, 68), which may be due to many AMLs often accompanied by T cell exclusion.The resistance to immunotherapy in AML, such as PD-1 blockade, remains one of the major challenges impeding its application in the future.So, a higher predictive response of immunotherapy in the high-risk subtype may bring a new perspective towards AML therapy, and VEN combined with blocking therapy of immune checkpoints is worthy of further exploration.
Considering VEN-resistance from the perspective of gene interaction, patients in the high-risk group presented a lower BCL2 expression, but higher levels of BCL2L1, BCL2A1, and MCL1.This means that the high-risk score often causes VENresistance by interacting with anti-apoptosis proteins, which is one of the main causes of VEN-resistance (9,69,70).Of other 5 senescence genes, as the co-variants of BCL-2, BAG3 works together with BAG1 to maintain the viability of myeloid cells, dysregulation of which could lead to physiological abnormalities (71).Notably, previous studies have indicated that BAG3 is associated with a poor AML prognosis (72), and involved in resistance to chemotherapy (73).Bound to the DR3 receptor, TNFSF15 often plays pro-inflammatory roles and regulates cytokine release (74).Furthermore, the TNFSF15/DR3 axis is involved in promoting apoptosis through Caspase pathways, but could also activate inhibitors of apoptosis proteins by regulating NF-kB pathways and inhibiting apoptosis (75).It is noteworthy that persistent inflammation could also promote the progression and resistance of tumors (76).Abnormal activation of SRC protein, one of the non-receptor tyrosine kinases, is closely related to the tumors' progression.Overexpression of which could lead to increased Src kinase activity, and play an important role in human cancer, including cell proliferation, differentiation, survival, and mortality (77).Besides, up-regulated GRK6 level is also associated with the progression and prognosis of colorectal carcinoma (78), and its role in AML is worth further exploration.Interestingly, AMLs sensitive to VEN-therapy are enriched in various gene sets of leukemic stem cells (LSCs), while the major enrichment of VEN-resistant AMLs is monocytic differentiation (27).As an essential regulatory molecule for activating LSCs (79), CDK6 is up-regulated significantly in VEN-sensitive AMLs.Consistent with these results, CDK6 levels are negatively associated with the risk score of VRSP-M, and highrisk AMLs are also enriched in the monocytic phenotype.In both mice and humans, aging is often accompanied by alteration of the monocyte function and increased production of classical monocytes expressing MHC II, which may help to explain why our high-risk subtype was mainly involved in monocyte differentiation and resistant to VEN-therapy (80).
In the context of oncogenes and clinical characteristics, patients in the high-risk subtype were frequently accompanied by mutations of NRAS, KARS, and TP53.Of note, AML with RAS mutation was associated with VEN-resistance and monocytic phenotype (27,81).As one of the most common proto-oncogenes in AML, a gain of function in KRAS/NRAS could activate the pathway of RAS/MAPK, and further lead to overexpression and increased stability of MCL-1 protein, which also plays a major role in VEN-resistance.Furthermore, RAS regulatory genes such as PTPN11 usually comutate with KRAS/NRAS mutation, which have been reported refractory to VEN monotherapy in AML (12, 82).In addition, a higher TP53 mutation may also help explain why the high-risk AML subtype exhibited a lower response to VEN and a poor prognosis (4,83,84).In general, our results demonstrate that more adverse prognostic features were presented in the high-risk group, including older patients, non-de novo AML, poor cytogenetics or adverse ELN risk, and TP53 mutation.All these results support why high-risk patients in our results have worse survival.
Our study provides a new perspective and potential therapeutic targets based on senescence aid to explore pathogenesis of VENresistance in AML; however, there are still several limitations in the current study.More independent AML cohorts are needed to validate it.Moreover, further investigation is needed to explore the underlying mechanisms.The risk score of the prognostic model is significantly associated with VEN-resistance, immune features, and immunotherapy response in AML.We also verified that G6PD and BAG3 could be effective biomarkers of VEN-resistance in vitro.
In conclusion, the 6-senescence genes prognostic model has significant meaning for the prediction of VEN-resistance, guiding personalized molecularly targeted therapies, and improving AML prognosis.

FIGURE 1 Work
FIGURE 1Work flow of the current study.

2 3
FIGURE 2 Identification of Venetoclax (VEN) resistance-related senescence genes (VRSGs).(A) Volcano map of differential expression genes (DEGs) between VEN-resistant and -sensitive samples through DEseq2 method.(B) Intersections of VEN-resistant related DEGs and Human related senescence genes (HRSGs).(C) Eighteen VRSGs were confirmed with AML prognosis in the modeling set using univariate Cox regression analysis.

FIGURE 4
FIGURE 4 Correlation of clinical features and VRSP-M in the modeling set.(A) Difference of age distribution between high-and low-risk subtypes based on VRSP-M.(B, C) Distribution of risk score between different risk of cytogenetics and French-American-British (FAB) subtype, including M0, M1, M2, M3, M4 and M5.*** indicate P<0.001.(D) Curve of gene set enrichment analysis (GSEA) between high-and low-risk subtypes in monocyte differentiation.(E) Coexpression of key genes in apoptosis pathways, the difference was compared by non-parametric test.IAP, inhibitor of apoptosis proteins.

7 6
FIGURE 7The different components of immune cells between high-and low-risk subtypes of VRSP-M.(A) The difference of 28 types of immune cells, evaluating by single-sample gene set enrichment analysis (ssGSEA).(B) The difference of immune and stromal score, calculated by ESTIMATE method.(C) Subtypes of macrophages between high-and low-risk AMLs, estimated by xCell algorithm.* indicate P<0.05, ** indicate P<0.01, *** indicate P<0.001.

8
FIGURE 8The correlation of VRSP-M's risk score with 8 immune checkpoints and predicted response of immunotherapy.(A-D) The relationship of VRSP-M's risk score with expression of 8 immune checkpoints, Exclusion, Dysfunction, and TIDE score, tested by Spearman method.(E, F) The difference of Exclusion, Dysfunction, TIDE score, and predicted response of immunotherapy between high-and low-risk subtypes.(G) Survival analysis of classification based on PDCD1 expression and VRSP-M's risk score.*** indicate P<0.001.

9 Screening
FIGURE 9 Screening Marker genes from VRSP-M.(A, B) The correlation analysis of VRSP-M's risk score with IC50 values (log10 transformation) and AUC values (log10 transformation) of VEN in Ex vivo data from Beat AML sample, tested by Spearman method.(C) The importance of G6PD and BAG3 in VENresistance, and the scores were evaluated using by CARE algorithm and CRISPR data.

10
FIGURE 10The experimental verification of the Key Markers.(A) The IC50 values of ABT-199 (Venetoclax) in different AML cell lines, tested by CCK-8 assay.(B, C) The mRNA levels of G6PD and BAG3 expression in different AML cell lines, using PCR method.(D, E) Protein levels of G6PD and BAG3 in different AML cell lines, tested by Western blot.

TABLE 1
Clinical characteristics of AML patients in TCGA cohort.

TABLE 2
Clinical characteristics of AML patients in Beat cohort.

TABLE 2 Continued
, Peripheral Blood; ELN, European LeukemiaNet.When taking genetic mutations into account, cases without genetic mutations data were excluded from analysis.The meaning of the bold value is P<0.05. PB