Gene expression data and the corresponding clinical data were obtained from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). To identify the key genes involved in venetoclax resistance, GSE223598 data (from OCI-Ly1 cells) and GSE252306 data (from SU-DHL-16 cells) were downloaded to analyze the differentially expressed genes (DEGs) in venetoclax-resistant cells compared to parental, and DEG threshold was set as: |log2-fold change (FC)|>2 and an adjusted P-value<0.05. A Venn diagram was generated showing 214 genes co-regulated in the two datasets. Among the 214 genes, 175 recorded in the GSE10846 dataset were selected as candidate genes for the prognostic models. There were 412 cases with complete gene expression and prognostic information in the GSE10846(training cohort), 470 in the GSE31312, and 221 in the GES87371 (validation cohort) datasets. The characteristics of patients in the GEO datasets have been described previously (10).

Gene expression data normalized to RMA in GSE10846 dataset were downloaded. Then, a LASSO regression analysis with ten-fold cross-validation recognized 36 of the 175 Bcl-2RGs as candidate gene sets for the prognostic signature. Multifactor Cox regression was used to identify the risk model containing Bcl-2RGs which was used to develop the Bcl-2RG signature. The risk score for patients in the cohort was calculated by the formula: Risk score =∑Ni=1=(exp×coef), where N is the number of model genes, exp represents the gene expression value of each gene, and coef represents the coefficient index. The ‘time ROC’ package was applied to draw the time-dependent receiver operating characteristic (ROC) curve and determined the area under the curve (AUC) of the 2-, 3-, and 5-year overall survival (OS) rates. The risk scores of the GSE31312 and GES87371 datasets were computed and classified into high and low groups to validate the prognostic model using Kaplan–Meier survival analysis.

The ‘WGCNA’ package was used to identify Bcl-2 inhibitor resistance gene clusters based on the GSE10846 dataset. The weighted adjacency matrix was converted into a topological overlap matrix according to the optimal soft threshold (β=10), and then hierarchical clustering analysis was performed to detect the correlation between gene modules (minmodulesize=30; mergecutheight=0.25). Interaction strength was assessed using the heatmap toolkit, and gene significance and module membership were calculated to assess the relationship between the module and resistance characteristics.

SU-DHL-8, SU-DHL-4, OCI-Ly3 and OCI-Ly10 cell lines were procured from Shanghai Zhong Qiao Xin Zhou Biotechnology Co., Ltd. Cells were maintained in Roswell Park Memorial Institute 1640 (RPMI) medium (Gibco, Waltham, MA, USA) supplemented with 10% fetal bovine serum (Gibco). Cell cultures were incubated at 37 °C in a humidified atmosphere containing 5% CO₂, following the standard protocols provided by the supplier. Cell viability was detected using a Cell Counting Kit-8 (CCK-8; Protein Bio, Nanjing, China) according to the manufacturer’s protocol. Optical density was measured at 450 nm using a full-spectrum microplate reader. Cell proliferation was analyzed using a BeyoClick™ EdU Cell Proliferation Kit with Alexa Fluor 594 (C0078L, Beyotime Biotechnology) according to the manufacturer’s protocol. Apoptosis was evaluated using an Annexin V-FITC Apoptosis Detection Kit (C1062L, Beyotime Biotechnology) as previously described (11).

Cellular proteins were extracted and homogenized using RIPA lysis buffer (NCM, Suzhou, China). Protein samples were separated using SDS-PAGE (Epizyme, Shanghai, China) and transferred onto polyvinylidene fluoride membranes. Membranes were subsequently incubated with 5% skim milk to block non-specific binding sites, treated with primary and secondary antibodies, and protein bands were visualized using enhanced chemiluminescence reagents (NCM, Suzhou, China) and detected using an imaging system.

Data analysis, including basic statistical comparisons and IC50 determination was performed using GraphPad Prism v9.5 (GraphPad Software, CA, USA). Flow cytometry data were analyzed using FlowJo v10.8.1 (Tree Star, OR, USA). Bioinformatics analyses, including LASSO regression analysis, Multifactor Cox regression, Kaplan–Meier survival analysis, WGCNA analysis, C-index and DCA analysis were performed using R software v4.4.1 with the limma v3.60.4, survival v3.6-4, WGCNA v1.73, survcomp v1.54.0, ggplot2 v4.0.1, ggDCA v1.2 package. Western blot images were quantified using ImageJ software (NIH). GEO2R were used for Screening for DEGs, and jvenn for Venn diagram drawing. Statistical significance was set at P< 0.05. significant.

The differentially expressed gene data from parental and venetoclax-resistant clonal cell lines (GSE223598 and GSE252306) were downloaded and analyzed to identify 214 Bcl-2RGs in DLBCL (Figures 1a, b). To clarify the relationship between Bcl-2RGs and the prognosis of DLBCL, 175 recorded Bcl-2RGs in the GSE10486 dataset (training cohort) and the OS of patients were analyzed. LASSO regression and ten-fold cross-validation identified 36 Bcl-2RGs as candidate gene sets for the prognostic signature (Figure 1c). Multifactor Cox regression results identified 13 Bcl-2RGs that were used to construct the prognostic model; their regression coefficients are shown in Figure 1d. The resistance score was obtained based on 13 Bcl-2RGs.The resistance score was estimated as follows:

Resistance Score = (0.372545 × BRSK1) + (0.304578 × ZNF486) + (0.293816 * MARCKS) + (0.304109 * ALPK1) + (0.094496 * GABRB2) + (0.278475* ZNF43) + (-0.156742 * CYTH4) + (- 0.383886 * XBP1) + (-0.349471 * ZNF141) + (-0.231252 * GOLGA8A) + (-0.2454 * LRFN3) + (- 0.110178 * GPR82) + (-0.083037 * PDE9A).

Patients in the GSE10846 dataset were classified into high- and low-drug resistance groups according to the median drug resistance score. Kaplan–Meier survival analysis revealed that patients with a high score (fen=1) had poorer OS than those with a low score (fen=0; Figure 1e). The AUCs of the 2-, 3-, and 5-year ROC curves were 0.667, 0.679, and 0.652, respectively, suggesting that the Bcl-2RG signature was highly useful for predicting outcomes (Figure 1f). Kaplan–Meier survival analysis showed that patients with high scores had shorter OS in the validation cohorts (GSE31312 and GSE87371; Figures 1g, h).

To further clarify the biological function of Bcl-2RGs in DLBCL, the WGCNA algorithm was used to identify hub genes associated with resistance (Figures 1a, b). Thirty gene modules were identified, of which ME red and ME pink were strongly associated with Bcl-2 inhibitor resistance (Figures 2a, b). Functional enrichment analysis revealed that genes in this module are primarily involved in B cell receptor signaling pathway, apoptosis, and cell cycle (Supplementary Figure S2), highlighting a network-level shift towards B cell growth and apoptosis. Although these modules contain several highly connected Hub genes (such as KRAS and BCL2), we focus on ALPK1 because it is not only strongly associated with Venetoclax resistance, but also serves as a key prognostic factor in the Lasso-Cox model. A Venn diagram was used to identify overlapping genes (ALPK1 and ZNF486) using WGCNA and LASSO analyses (Figure 2b). The ALPK1/NFκB axis plays an important role in innate immunity; however, less was known about the ALPK1/NFκB axis in DLBCL. We conducted experiments to verify the role of ALPK1 in the development and drug resistance of DLBCL. ALPK1 expression in DLBCL cell lines and patients with DLBCL was higher than that in healthy donors (controls; Figure 2c). The ALPK1 inhibitor, ALPK1-IN2, reduced SU-DHL-8 and OCI-Ly10 cell viability with IC50s of 11.39 and 11.43 uM, respectively (Figure 2d). Cell proliferation in SU-DHL-8 and OCI-Ly10 cells was evaluated by EdU incorporation assays, revealing a significantly lower percentage of EdU⁺ cells in the ALPK1-IN2 treated cells relative to untreated controls (Figure 2e; P< 0.05). DLBCL cell apoptosis increased in the ALPK1-IN2 treatment group compared to that in the control group (Figures 2f, g; P< 0.05).

To determine whether ALPK1 is involved in venetoclax resistance in DLBCL and its mechanism of action, inhibitors of ALPK1 and Bcl-2 were used to treat DLBCL cell lines (SU-DHL-8 and OCI-Ly10). Cell viability was lower in the ALPK1-IN2-treated group than the control group The IC50s of venetoclax was reduced from 5.943 uM to 1.189 uM in SU-DHL-8 cells and from 5049 nM to 99.07 nM in OCI-Ly10 cells,(Figure 3a), indicating that ALPK1-IN2 might reverse venetoclax resistance. Venetoclax treatment significantly increased apoptosis in the ALPK1-IN2 group compared to the control group (Figure 3b). The percentage of EdU⁺ cells did not change after treatment with venetoclax in SU-DHL-8 or OCI-Ly10 cells. Combined treatment with ALPK1-IN2 reduced the percentage of EdU⁺ cells compared to control group and venetoclax treatment alone (Figure 3c, d).

ALPK1 is a protein kinase that plays an important role in activating the innate immune response to bacteria via NFκB signaling. NFκB is a hub gene in the BCR signaling pathway, which promotes the occurrence and progression of B-cell lymphoma by promoting proliferation and inhibiting apoptosis. Therefore, we hypothesize that ALPK1 might sensitize DLBCL to venetoclax by regulating NFκB signaling. As hypothesized, the expression of NFκB and its downstream signaling protein, BCL-XL, were reduced in the ALPK1-IN2 group compared to the control group, and significantly lower in the ALPK1-IN2 plus venetoclax group (P< 0.05; Figure 3e). The proliferation-associated protein PCNA was significantly downregulated in the ALPK1-IN2 plus venetoclax group compared to the other groups (P< 0.05; Figure 3e).