Systematic Analysis of E3 Ligase-Related Genes Identified UBE2L3 as a Prognostic Biomarker Associated With Drug Resistance in Acute Myeloid Leukemia

Introduction

Acute myeloid leukemia is a common and highly heterogeneous malignant disease of the hematologic system. Over the past 30 years, the incidence of AML has seen a global increase, presenting a significant risk to public health due to its poor prognosis.¹ Drug resistance and relapse are important factors affecting the survival of AML patients. With the deepening of research, more and more genetic landscape features of AML have been discovered, and the molecular mechanisms of drug resistance have become clearer.² The ubiquitin-proteasome pathway regulates nearly all cellular activities through protein degradation and non-proteolytic events; it is closely associated with the occurrence and development of cancer and plays a central role in regulating protein levels and controlling crucial cellular processes such as the cell cycle, gene expression, cellular oxidative stress, proliferation, apoptosis, and DNA repair.^3,4 Thus, this pathway has a major impact on tumorigenesis and other pathology-related events, inhibition of the pathway could represent a novel strategy for the treatment and overcoming of drug resistance in chemoresistant malignancies.⁵

Over 600 E3 ligases have been identified in the human body, but only a few of targeted degron instances have been identified so far.^6,7 E3 ubiquitin ligases play an important role in the ubiquitin proteasome-pathway, and the E3 ubiquitination step determines the overall specificity of ubiquitination of target proteins, with activated ubiquitin usually transferred to lysine residues in the substrate. Depending on the structure, E3 ligases mainly include HECT, RING Finger and U-box structural domain E3s, which stimulate the transfer of ubiquitin and ubiquitin-like proteins by direct or indirect mechanisms.⁸ E2 ligases are essential partners of E3 ligases because all E3 ligases bind to E2~Ub thioesters and catalyse the transfer of ubiquitin from E2 to the substrate lysine via an amino lysis reaction, thereby determining each protein final destination.9,¹⁰ After bortezomib achieved previously unmatched efficacy in the treatment of multiple myeloma, the study of the ubiquitin-proteasome pathway in hematological tumor has received increasing attention.¹¹ Several mechanisms involving E3 ubiquitin ligase-related genes in AML have been gradually uncovered. For example, Cbl-b, by ubiquitinating and promoting the degradation of Siva1, activates the P53 pathway, potentially overcoming FLT3 inhibitor-resistant AML.¹² Triad1 antagonizes Mll-Ell induced AML progression in vivo and can act as an AML suppressor.¹³ WWP1 promotes cell cycle entry and survival of AML progenitors, maintaining leukemia cell growth; its inactivation might promote autophagy and slow disease progression.¹⁴ RNF5 regulates transcription in AML cells, promoting their proliferation, and its inhibition enhances ER stress-induced AML cell apoptosis.¹⁵ UBE2C regulates the PI3K/AKT signaling pathway to affect AML cell survival and may modulate the sensitivity of leukemia cells to ferroptosis.¹⁶ These findings suggest potential new therapeutic targets. Inhibitors of certain ubiquitin ligases have been discovered and have shown relatively good efficacy in clinical applications.^17–19

UBE2L3, as an E2 ubiquitin ligase involved in the ubiquitination of various substrate proteins. Studies have demonstrated that UBE2L3 is aberrantly expressed in cancer cell lines and tumor tissues, indicating its involvement in various cancer-related signaling pathways. Tao NN et al showed that UBE2L3 overexpression in hepatocellular carcinoma promotes apoptosis evasion by inhibiting the GSK3β/p65 pathway.²⁰ In cervical cancer, UBE2L3 and E6AP promote tumor cell proliferation by degrading p53, which blocks p53-mediated growth arrest and apoptosis.²¹ Another study indicates that the fusion of UBE2L3 with KRAS may facilitate metastatic progression in rare subgroups of prostate cancer.²² UBE2L3 drives tumor initiation, progression, and metastasis by regulating cancer-associated signaling pathways and other proteins not previously discussed.²³ However, its role in hematological malignancies remains unclear.

An increasing number of studies have shown that the ubiquitin-proteasome pathway plays an important role in tumor progression. Currently, there has been no systematic study of E3 ligases in AML, making our research highly significant. Therefore, in-depth study of the role of E3 ligase-related genes in drug resistance and progression of AML, as well as the establishment of relevant prognostic models, are of great significance for the treatment of AML. In addition, we validated the expression of UBE2L3 in K562 and HL60 cell lines and their drug-resistant strains, as well as in clinical samples; finally, we determined that UBE2L3 may be a reliable biomarker of drug resistance and poorer prognosis in AML.

Materials and Methods

Data Acquisition and Preprocessing

RNA-sequencing data were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov) for the TCGA-LAML project using the STAR workflow. Such data, were extracted in TPM format, data standardization uses log2(value+1). Beat AML 2.0 data from www.vizome.org.² Prior to analysis, data preprocessing steps were performed, including data cleaning, removal of outliers, batch effect correction, and principal component analysis.

Identification of Differentially Expressed E3 Ubiquitin Ligase-Related Genes

The E3 ubiquitin ligase itself and related genes involved in its functional regulation are regarded as E3 ubiquitin ligase-related genes. These genes are from https://www.genecards.org and https://www.uniprot.org. Transcriptome sequencing was performed on wild-type and Adriamycin-resistant HL60 cells. Significance thresholds of P < 0.05 and |log2 fold change (FC)| > 1.5 were used as the screening criteria. The results showed that 3209 genes had their expression upregulated and 2509 genes downregulated. A total of 154 E3 ubiquitin ligase-related differentially expressed genes were identified by taking the intersection with the previously searched genes.

Enrichment Analysis of Differentially Expressed Genes

The input list of differentially expressed molecules was converted into IDs using the “org.Hs.eg.db“ package, followed by enrichment analysis using the ”cluster Profiler” package. The ggplot2 package was used to visualize the results of the enrichment analysis.

Identification of Prognostic Genes

One hundred and seventy-three AML samples were curated from the TCGA database, of which 139 possessed relatively complete clinical data. To investigate the relationship between the expression levels of E3 ubiquitin ligase-related genes and overall survival (OS) of AML patients, we performed univariate Cox regression analysis using the “survival” package. The significance filter criterion was set at p < 0.05 for further analysis. Next, we used LASSO Cox regression to eliminate gene collinearity and reduce the number of genes. Finally, we performed multiple Cox regression analysis based on univariate Cox regression.

Construction and Validation of a Prognostic Model Based on Differentially Expressed E3 Ubiquitin Ligase Genes

The risk score was calculated based on AML mRNA expression data from the training set (n=139), using the following formula:.Patients with AML were divided into high- and low-risk groups based on the median risk score, and the OS of the two groups was analysed. To enhance the reliability of the model, validation was performed using the test set (n=614) and the risk score was calculated using the aforementioned formula. Patients with AML in the test set were also divided into high- and low-risk groups based on the median risk score, and the OS of the two groups was compared. Subsequently, univariate and multivariate Cox regression analyses were performed to determine whether the risk score was an independent prognostic factor for OS in patients with AML in the training set. The covariates included age, sex, white blood cell count, molecular mutation and cytogenetic risk.

Construction of a Nomogram and Calibration Curves, Decision Curve Analysis, Survival Curves (Kaplan-Meier Plots)

The “RMS“ package in R software (R Foundation for Statistical Computing, Vienna, Austria) was used to construct a nomogram to predict individual survival probabilities and to plot calibration curves for predicting 1, 2, and 3-year survival rates of patients with AML. Use the glm function to construct two-category logistic models respectively, and use the “rmda“ package to calculate the corresponding net rate of return and visualize it. Proportional hazards hypothesis testing and fitting survival regression were performed using the “survival” package, and the results were visualized using the ”survminer” and ”ggplot2” packages.

Patient and Ethics Approvals

The criteria for defining newly diagnosed or relapsed/refractory AML in this study were based on the NCCN (National Comprehensive Cancer Network) guidelines. Induction chemotherapy for patients with AML followed the DA/IA (idarubicin, 8–12mg/m²/day, for 3 days or daunorubicin, 45–60mg/m²/day, for 3 days; and cytarabine,100 mg/m²/day, for 7 days) regimen, Consolidation chemotherapy for patients with AML followed the medium-dose cytarabine (cytarabine,2.0g/m²/q12h, for 3 days), and the assessment of treatment efficacy was also based on the NCCN guidelines. We gathered data from 49 AML patients with a follow-up duration of at least 6 months, 8 of whom relapsed. During the follow-up period, none of the patients who underwent allogeneic hematopoietic stem cell transplantation experienced relapse. This study was approved by the Ethics Committee of Affiliated Chuzhou Hospital of Anhui Medical University, First People’s Hospital of Chuzhou and all procedures were performed in accordance with the ethical guidelines of the Declaration of Helsinki. Informed consent was obtained from each participant before their participation in the study.

Leukemia Cell Line Culture

K562 and HL60 cell lines were obtained from the Hematology Laboratory of the Second Affiliated Hospital of the Anhui Medical University. The K562/Adriamycin-resistant (ADR) and HL60/Adriamycin-resistant (ADR) cell lines were generously provided by Professor Li Zhao from the Center Laboratory of the First Hospital of Lanzhou University.²⁴ The cells were cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum, 100 U/mL penicillin, 100 mg/mL streptomycin, and 2.05 mm L-glutamine. All cells were cultured in a humidified incubator at 37 °C with 5% CO₂.

Real-Time PCR and Western Blotting

Total RNA was extracted using the TRIzol method and the total RNA was reverse transcribed using the HI Script III RT Super Mix for qPCR (+gDNA) kit (Vazyme Biotech, #R323-01). Real-time PCR was performed using the Roche Light Cycler 480 II PCR amplification system with a mixture of cDNA, gene-specific primers, and 2× ChamQ Universal SYBR qPCR Master Mix (#Q711-02; Vazyme Biotech). The primer sequences used were as follows: UBE2L3: F-CGCTGGCGGGCAAAGACT; R-GCCCCTTTTCGTCGATGTTTG. β-actin: F: TGGCACCCAGCACAATGAA; R: CTAAGTCATAGTCCGCCTAGAAGCA. Cells were lysed in RIPA buffer, and proteins were extracted and denatured. Proteins were then separated by SDS-PAGE gel electrophoresis and transferred onto nitrocellulose membranes. After blocking with skim milk, the membranes were incubated overnight with specific primary antibodies (Proteintech, #14415-1-AP, #66009-1-Ig) followed by incubation with horseradish peroxidase-conjugated secondary antibodies for 1 hour. Finally, an enhanced chemiluminescence system (Tanon, China) was used for visualization.

Statistical Analysis

SPSS software version 21.0 (SPSS Inc, Chicago, IL, United States) and R software version 4.3.0 (R Foundation for Statistical Computing, Vienna, Austria) were used to analyse the data. Statistical significance was set at a two-sided P < 0.05.

Results

Identification of Differentially Expressed E3 Ubiquitin Ligase-Related Genes in AML

Our detailed workflow is shown in (Figure 1). To identify the key differentially expressed genes associated with drug resistance in AML, transcriptome sequencing to evaluate the gene expression profiles of wild-type HL60 cells (HL60-WT) and Adriamycin-resistant HL60 cells (HL60/ADR) cells were performed. To identify as many differentially expressed genes as possible, we used a significance criterion with an adjusted p-value < 0.05, and |log2 (fold change) | ≥ 1.5 to evaluate differential gene expression. We identified 3209 upregulated and 2509 downregulated genes (Figure 2A). The intersection of the differential genes from transcriptome sequencing, the E3 ubiquitin ligase-related genes searched in the GeneCards database and the UniProt database (Figure 2B).

Figure 1 The flowchart of our research process. TCGA, The Cancer Genome Atlas database; BeatAML2.0, Beat Acute Myeloid Leukemia 2.0. The UniProt and GeneCards database website can be accessed at uniprot.org. and genecards.org. Transcriptome sequencing data from Shanghai Sanshu Biotechnology Co., Ltd (Shanghai, China).

Figure 2 Construction of a risk prognosis model based on E3 ubiquitin ligase-related genes in the TCGA cohort. (A) Volcano plot of the transcriptome sequencing of HL60-WT and HL60/ADR cells. (B) Venn diagram showing differential genes screened from GeneCards data, UniProt database, and sequencing data. (C) Functional enrichment analyses of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). (D) Univariate Cox regression analyses of the selected 154 genes in the training set based on the TCGA cohort, with P < 0.05 indicating statistical significance. (E) LASSO analysis identifying 5 survival-related genes. (F) Cross-validation for parameter selection in the LASSO regression. (G) Multivariate Cox regression analyses based on the univariate Cox regression analysis results. (H) Heat map of 5 differential genes in prognostic models, with red indicating high expression, blue indicating low expression.

Functional Enrichment Analysis

To better understand the functions of differentially expressed E3 ubiquitin ligase-related genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed. The analysis revealed that the differentially expressed genes associated with E3 ubiquitin ligases were primarily involved in the ubiquitin-proteasome pathway, which was consistent with the screening results (Figure 2C).

Construction of a Prognostic Model Based on E3 Ubiquitin Ligase-Related Genes in the Train Set

As shown in (Figure 2D), we screened out 25 E3 ubiquitin ligase-related genes with P < 0.05 using univariate Cox regression analysis. On the basis of the univariate Cox regression, we then performed LASSO regression analysis (Figure 2E and F). Next, we constructed the prognostic E3 ubiquitin ligase-related model using HBP1, RNF130, RMND5B, TRIM32, and UBE2L3 by LASSO regression. Finally, we performed multivariate Cox regression analysis and identified three E3 ubiquitin ligase-related genes, two of which were potential risk genes and one of which was potential protective genes (Figure 2G). Gene expression heat maps of the model based on transcriptome sequencing, as shown in (Figure 2H).

By establishing the model, a prognostic index for the samples was constructed using the following formula: HBP1 expression level × (−0.0355832045) + RNF130 expression level × (−0.0319662324) + RMND5B expression level × 0.0546280923 + TRIM32 expression level ×0.0008787761+ UBE2L3 expression level ×0.1736809863. Based on the median risk score, divide patients in the training set and test set into high-risk and low-risk groups. Compared with the low-risk group, the high-risk group had a higher mortality rate and shorter survival time. A high prognostic index score correlated with poor disease prognosis (Figure 3A–D). In the high-risk group, high expression of RMND5B, TRIM32 and UBE2L3, and low expression of RNF130 and HBP1 were observed (Figure 3E and F). The KM curve demonstrated a relatively poor prognosis in the high-risk group (P < 0.05, Figure 3G and H). Time-dependent ROC curve analysis showed the accuracy of OS prediction at 1, 3, and 5 years as follows: Training set:0.75514 (95% CI 0.6697–0.8406), 0.75983 (95% CI 0.6663–0.8533), and 0.74496 (95% CI 0.6321–0.8578, Figure 3I), Test set:0.65676 (95% CI 0.6106–0.7029), 0.5996 (95% CI 0.548–0.6512), and 0.59862 (95% CI 0.5436–0.6537, Figure 3J) respectively. These results suggested that the E3 ubiquitin ligase-associated features in our model may contribute to predicting the prognosis of AML with reliable accuracy.

Figure 3 Validation of prognostic models and independent prognosis analysis. (A) Patients in the training set are divided into two groups based on the median risk score. (B) Patients in the test set are divided into two groups based on the median risk score. (C) Survival status of high- and low-risk groups of patients with AML in training set. (D) Survival status of high- and low-risk groups of patients with AML in test set. (E) Heatmap showing the expression of five E3 ubiquitin ligase-related genes in training set. (F) Heatmap showing the expression of five E3 ubiquitin ligase-related genes in test set. (G) Kaplan-Meier curves showing the overall survival of high- and low-risk group patients in training set. (H) Kaplan-Meier curves showing the overall survival of high- and low-risk group patients in test set. (I) Validation of the predictive efficiency of the risk score through ROC curve analysis in training set. (J) Validation of the predictive efficiency of the risk score through ROC curve analysis in test set. (K) Univariate independent prognosis Cox regression analysis of risk score and indicated clinical characteristics. (L) Multivariate independent prognosis Cox regression analysis of risk score and indicated clinical characteristics.

Independent Prognostic Analysis of Risk Score and Clinical Features

To validate whether the risk score and clinical features could serve as independent prognostic factors, univariate and multivariate independent prognostic analyses were conducted. Univariate analyses revealed significant correlations between age, cytogenetic risk stratification, risk score, and OS in patients with AML. Multivariate independent prognostic analysis indicated that the risk score may be an independent predictive factor (Figure 3K and L).

Construction of Nomogram and Calibration Curves

A nomogram was created incorporating age, sex, white blood cell count, cytogenetic risk stratification, molecular mutation and risk scores (Figure 4A). The nomogram demonstrated the risk score as an important factor among various clinical parameters. Additionally, we constructed calibration curves that showed a close match between the nomogram and the actual survival rates of patients with AML (Figure 4B–D). The decision curve analysis (DCA) and ROC curves showed a significant advantage of the risk model in predicting survival time in the training cohort (Figure 4E–H).

Figure 4 Nomogram to predict the survival probability of patients with AML. (A) Nomogram combining the risk score with AML prognostic factors such as cytogenetic risk stratification. (B-D) calibration plots for predicting the 1-, 2-, and 3-year overall survival of patients. (E) ROC curve for survival prediction using the risk score and other variables (age, gender, white blood cell count, cytogenetic risk stratification). (F-H) Decision curve analysis of the nomogram.

ROC Curve and Kaplan-Meier Plots of Prognostic Genes

Drawing ROC curves based on patient survival outcomes, we can find that UBE2L3 has better specificity and sensitivity for distinguishing the prognosis of AML patients than the other 4 genes in the model (Figure 5A–E). Kaplan-Meier plots were generated to assess the association between the expression of E3 ubiquitin ligase-related genes in the prognostic model and AML prognosis. The results revealed that patients with AML with high UBE2L3, RMND5B and TRIM32 expression exhibited a poorer prognosis (P < 0.05, Figure 5F–H), whereas those with high RNF130 and HBP1 expression showed a better prognosis (P < 0.05, Figure 5I and J).

Figure 5 Prediction of disease outcome and survival analysis by model genes. (A-E) The ROC curves of HBP1, RNF130, RMND5B, TRIM22, and UBE2L3 for predicting survival outcomes. (F-J) The survival curves of UBE2L3, RMND5B, TRIM32, RNF130, HBP1.

UBE2L3 Is Expressed in Cell Lines and in Patients With AML

In the K562 and HL60 cell lines, UBE2L3 expression in the mRNA and protein levels of the drug-resistant cell lines was higher than that of the wild-type cell lines (P < 0.05, Figure 6A–C). Compared with patients with newly diagnosed AML, the expression of UBE2L3 mRNA in the bone marrow of relapsed and refractory patients was increased (P < 0.05, Figure 6D). Through the Kaplan-Meier plots, we also can find that high expression of UBE2L3 has a poor prognosis (p<0.05, Figure 6E). The characteristics of the patients is presented in Supplemental Table 1.

Figure 6 Validation of cell lines and clinical samples. (A) Relative expression levels of UBE2L3 mRNA in wild-type leukemia cell lines (K562-WT and HL60-WT) and Adriamycin-resistant cell lines (K562/ADR and HL60/ADR). (B and C) Relative expression levels of UBE2L3 protein in wild-type leukemia cell lines and Adriamycin-resistant cell lines. (D) Changes in UBE2L3 expression newly diagnosed and after relapse in 8 patients with AML. (E) The survival curves of UBE2L3 in 49 patients with AML. **P<0.01, ***P<0.001.

Discussion

AML is a highly heterogeneous hematologic malignancy with significant variations in pathology, molecular basis, clinical treatment, and response to therapy.²⁵ The role of E3 ubiquitin ligases in the ubiquitin-proteasome pathway has gathered increasing attention in research regarding tumor onset and progression.^26–32

In this study, an AML prognostic risk model was first constructed through a systematic analysis of E3 ubiquitin ligase-related genes. A risk model was then constructed and validated using univariate Cox regression and LASSO analyses. This model included five genes (HBP1, RNF130, RMND5B, TRIM32 and UBE2L3). Ultimately, HBP1 and RNF130 were identified as protective factors for the disease, whereas RMND5B, TRIM32 and UBE2L3 was determined to be a risk factor for AML drug resistance and prognosis. HMG box protein 1 (HBP1) is a transcription factor and a potent cell cycle inhibitor in normal and cancer cells.³³ HBP1 is ubiquitinated by the CTLH E3 ubiquitin ligase complex, leading to its degradation and thus preventing cell cycle exit in G1 phase.³⁴ As an E3 ubiquitin ligase, RNF130 may be involved in programmed cell death, immune regulation, tumor progression, etc.^34,35 As a core component of the CTLH E3 ubiquitin-protein ligase complex, RMND5B selectively accepts ubiquitin from UBE2H and mediates ubiquitination and subsequent proteasomal degradation of the transcription factor HBP1.³⁶ TRIM32 is a member of the tripartite motif family, which mediates and regulates many physiological and pathophysiological processes such as cell growth, immunity, carcinogenesis, etc. through its E3 ubiquitin ligase activity.³⁷ UBE2L3 is one of the many E2 ubiquitin conjugating enzymes, it plays a biological role in conjunction mainly with HECT or RBR (RING-in- Between-RING) E3.²³ In recent years, more and more studies have shown that UBE2L3 is significantly related to the occurrence and development of tumor, involving multiple signaling pathways, such as NF-κB signaling pathway, GSK3β/p65 signaling pathway, p53 signaling pathway, DSB repair pathways.^20,38–40

Drug resistance affects the chemotherapy effect of AML patients, thereby affecting the prognosis. The drug resistance mechanism of cells mainly includes the following aspects: first, the impact of cell membrane transport proteins on the uptake of water-soluble drugs; second, changes in the cell cycle, increased repair of DNA damage, reduced cell apoptosis, and changes in drug metabolism, which affects the killing of cells by cytotoxic drugs; finally, energy-dependent efflux of hydrophobic drugs is increased.⁴¹ Sang Ah Yi et al found that in the DNA damage response induced by doxorubicin, nuclear translocation of HP1γ and binding on the UBE2L3 promoter mediate UBE2L3 silencing and increase p53 stability, which improves chemosensitivity in cisplatin-resistant cervical cancer.^21,39 Depletion of UBE2L3 stabilizes 53BP1, sensitizing cells to DNA damage and further enhancing the effects of radiotherapy or chemotherapy.⁴⁰ However, when UBE2L3 is overexpressed, the multifunctional protein p27 is stabilized and promotes cell migration.⁴² The study demonstrated that overexpression of UBE2L3 significantly activated the nuclear factor kappa B (NF-κB) signaling pathway through promoting NF-κB p65 nuclear translocation and the ubiquitination and degradation of IκBα protein that can head and neck squamous cell carcinoma progression.³⁸ UBE2L3 interacts with various E3 ubiquitin ligases to regulate the NF-κB signaling pathway. In response to IL-1R stimulation, TAX binds to UBE2L3, thereby regulating the degradation and activation of NF-κB.⁴³ In the TNFR pathway, UBE2L3 binds to LUBAC to initiate the synthesis of linear ubiquitin chains, which ultimately activate NF-κB.⁴⁴ In addition, E6AP can interact with UBE2L3 to promote the ubiquitination and degradation of p53.⁴⁵ Increasing evidence shows that UEB2L3 plays an important role in the occurrence, development, and invasion of tumor.^21,27,38,46 However, most of these studies have been conducted in solid tumors, and research on UBE2L3 in AML remains absent.

In addition, we further confirmed through leukemia cell lines and clinical samples that the relative expression of UBE2L3 in drug-resistant cell lines was significantly increased compared with wild type. In clinical samples, the expression of UBE2L3 in relapsed and refractory patients was significantly higher than that in newly diagnosed patients. However, our study has some limitations, including the lack of large-scale clinical cohort studies to validate the predictive efficacy and broad applicability of our model. In subsequent studies, we plan to increase the sample size and include more case samples to further explore the role of UBE2L3 in drug resistance, relapse, and refractory AML patients.

Conclusion

In this study, we developed a gene prognostic model based on HBP1, RNF130, RMND5B, TRIM32, and UBE2L3 to effectively predict the prognosis of AML patients. In addition, in vitro cell experiment results and clinical sample studies have shown that high expression of UBE2L3 is a reliable biomarker for drug resistance and poor prognosis of acute myeloid leukemia.