We obtained a uniformly standardized pan-cancer dataset (TCGA, GTEx) and retrieved UBD gene expression data across tumor and normal tissues from the UCSC Xena Browser (http://xenabrowser.net/) (11). Mutational profiles, copy number alterations (CNAs), and gene fusion events were analyzed using the cBioPortal platform (https://www.cbioportal.org/) based on the TCGA Pan-Cancer dataset (12). Additionally, gene-level Copy Number Variation (CNV) data and Level 4 gene expression profiles for all TCGA samples, processed via GISTIC software (DOI: 10.1186/gb-2011-12-4-r41), were downloaded and integrated from the GDC portal (https://portal.gdc.cancer.gov/) (13, 14). Differential UBD expression between tumor and normal tissues was assessed using R software. Raw expression matrices underwent log2(x + 0.001) transformation to stabilize variance, followed by batch effect correction using the ComBat-seq algorithm. Statistical significance was defined as P < 0.05.

The Kaplan–Meier (KM) survival analysis technique is a commonly - used statistical method for comparing survival differences among diverse cohorts. In the present study, we utilized the survival package within R to carry out KM survival analysis on patient groups with high and low UBD expression across 33 cancer types, which covered overall survival (OS), disease - specific survival (DSS), and progression - free interval (PFI) (31). The Cox regression model was applied to calculate p - values and assess the statistical significance of survival differences. Through the use of the Survminer and ggplot2 packages, we computed hazard ratios (HRs), 95% confidence intervals (CIs), and p - values, and presented these results visually. The KM survival analysis utilized the “surv_cutpoint” function from the R package survminer to determine the optimal cutpoint for continuous variables. This function calculates maximally selected rank statistics based on the maxstat package, and identifies the optimal cutoff value by maximizing this statistic to define the most discriminative grouping.

We systematically evaluated the influence of UBD expression levels on the extent of immune cell infiltration across a diverse array of cancer types within the TCGA database using seven methods, namely XCELL, CiberSort_ABS, CiberSort, EPIC, QUANTISEQ, MCPCOUNTER, and TIMER. The R packages utilized for quantifying immune infiltration included CIBERSORT, xCell, IOBR, MCPcounter, and the quanTIseq package (15).

The study “The Immune Landscape of Cancer” conducted a large-scale immunogenomic analysis of over 10,000 tumor samples from 33 different cancer types available in the TCGA database (16). In this pan-cancer analysis, researchers identified six distinct immune subtypes based on the following criteria: macrophage or lymphocyte markers, the ratio of Th1 to Th2 cells, the range of tumor intergenetic heterogeneity, aneuploidy, neoantigen burden, the overall cell atlas, the expression of immune regulatory genes, and prognosis. The six subtypes are described as follows:

C1 (Wound Healing): Elevated expression of angiogenesis-related genes, high proliferative fraction, and Th2-skewed adaptive immune infiltration. C2 (IFN-γ Dominant): Highest M1/M2 macrophage polarization, strong CD8+ T cell signaling, and similar to C6, the highest T cell receptor (TCR) diversity. C3 (Inflammatory): Elevated Th17 and Th1 gene expression, inability to effectively restrain tumor cell proliferation, and like C5, fewer aneuploidies and overall copy number alterations compared to other subtypes. C5 (Immune Silent): The least lymphocytic infiltration, highest macrophage response, and M2 macrophage predominance. C6 (TGF-β Dominant): A smaller group with the highest TGF-β signature and high lymphocytic infiltration. It has an equal distribution of type I and type II T cells.

Data on somatic variants and DNA copy number alterations (CNA) for a pan - cancer analysis were obtained from the cBioPortal website (17). The Spearman correlation between UBD expression levels and DNA copy number alterations was calculated to evaluate the association between somatic copy number alterations (SCNA) and UBD expression. The results were presented in the form of a heatmap.

Patients from the TCGA dataset were grouped into high and low UBD expression groups according to their UBD expression levels. Gene Set Enrichment Analysis (GSEA) was used to assess the modulation of Hallmark gene sets and KEGG pathway gene sets across different expression levels in various tumors (18). Additionally, the relationship between UBD mRNA levels and protein expression measured by Reverse Phase Protein Array (RPPA) in the TCPA database was evaluated using Rank - based association analysis. The results across all tumors were visualized by heatmaps (19).

The z - score algorithm proposed by Lee et al. was used to reflect the activity of specific pathways by integrating the expression of feature genes (20). Gene sets containing genes related to Angiogenesis, Epithelial - to - Mesenchymal Transition (EMT), Cell Cycle, Apoptosis, Hypoxia, Inflammation, Invasion, Metastasis, Proliferation, Quiescence, Stemness, Differentiation, DNA Damage, and DNA Repair were subjected to the z - score algorithm implemented in the GSVA R package. These gene sets related to the aforementioned tumor pathways were derived from the Cancer Single - cell Atlas (CancerSEA) database (21).

Total RNA was extracted using TRIzol reagent (Ambion, USA) following standard protocols. cDNA synthesis was performed using PrimeScript™ RT Master Mix (Takara, Japan) according to the manufacturer’s instructions. Gene expression analysis was conducted via real-time PCR with ChamQ SYBR qPCR Master Mix (Vazyme, China), with the 2^-ΔΔCT method applied for relative quantification, using GAPDH as the endogenous control.

Specific primer sequences were designed as follows:

UBD: Forward 5′-CCGTTCCGAGGAATGGGATTT-3′, Reverse 5′-GCCATAAGATGAGAGGCTTCTCC-3′. GAPDH: Forward 5′-AACAGCCTCAAGATCATCAGC-3′, Reverse 5′-GGATGATGTTCTGGAGAGCC-3′.

Three independent experimental replicates were performed to ensure data reliability, with cycle threshold (Ct) values averaged across technical triplicates. Quantitative measurements of target gene expression were normalized against the housekeeping gene GAPDH to account for potential variations in RNA input.

Protein extraction was performed using RIPA lysis buffer. Protein concentration determination was conducted with a BCA protein assay kit (Beyotime Biotechnology) following standard protocols. Electrophoretic separation was carried out on 10% SDS-PAGE with 50 μg protein samples loaded per well. After transferring the resolved proteins onto PVDF membranes, blocking was achieved with 5% non-fat milk in TBS-T to minimize background signals.

Primary antibodies included:

Membranes were incubated with HRP-conjugated goat anti-rabbit IgG secondary antibody (Abcam;1:5000) at 4°C for 60 minutes after thorough TBS-T washing. Protein band visualization was achieved using an enhanced chemiluminescence detection system, which optimizes the sensitivity of chemiluminescent reactions through substrate optimization.

The esophageal cancer cell line TE-11 was obtained from Shanghai Fuheng Biological Technology Co., Ltd. The TE-11 cell line was maintained in Dulbecco’s Modified Eagle’s Medium (Gibco, Grand Island, USA) containing 10% fetal bovine serum (Cyagen, Suzhou, China).

For the assessment of cell proliferation and colony-forming potential, a Colony Formation Assay was conducted. Logarithmically growing cells stably expressing the gene of interest were seeded at 500 cells/mL in 6-well plates, with 1 mL of cell suspension per well, and cultured under 37°C, 5% CO2. The medium was refreshed every two days based on color change to maintain nutrient supply. After 7 days, colonies were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, and quantified using ImageJ for colony number and area analysis.

Cell proliferation and viability were assessed via Cell Counting Kit-8 (CCK-8) assay. Cells were seeded at 1*10^3 cells/well in 96-well plates and incubated at 37°C, 5% CO2. Over five days, CCK-8 solution was added daily, followed by a 1-hour incubation, and absorbance measured at 450 nm using a microplate reader.

Transwell migration assays were performed to evaluate cell migration capabilities. Cells were suspended at 2x10^5/mL in 1% FBS medium, and 100 μL of this suspension was added to Transwell chambers prehydrated with serum-free medium. Chambers were placed in 24-well plates containing 600 μL of 20% FBS medium. After 24 hours, cells were fixed with methanol, stained with 0.1% crystal violet, and non-migrated cells were gently removed. Images were captured under a microscope.

Scratch assays were utilized to investigate cell migration dynamics. Cells were plated at 6x10^5 cells/well in 6-well plates and allowed to grow until >90% confluence. Post serum-starvation for 4 hours, scratches were made with a 200 μL pipette tip, washed with PBS, and maintained in serum-free medium for 24 hours. Wound closure was monitored and analyzed using ImageJ.

This retrospective cohort study was approved by the Institutional Ethics Committee of Lu’an People’s Hospital (Ethics Approval No. 2023LLKS012). Formalin-fixed paraffin-embedded (FFPE) tissue blocks were systematically collected from patients with ESCA who underwent curative resection between January 2019 and August 2025. The inclusion criteria were as follows: histologically confirmed primary esophageal carcinoma, tumor cellularity of at least 30% in representative sections, and complete clinicopathological records. Exclusion criteria included receipt of neoadjuvant therapy and insufficient tissue for comprehensive analysis. Detailed patient characteristics are provided in Supplementary Table 1. Immunohistochemical analysis of UBD expression was performed using Rabbit anti-UBD polyclonal antibody (Thermofisher, #PA5-102790;1:2000). Automated optical density (AOD) values were calculated as the ratio of positive staining area to total tissue area using ImageJ software (National Institutes of Health, v1.53). Five randomly selected 400× fields per specimen. Patients were stratified into three distinct prognostic groups based on AOD tertile distribution.

The LIMMA differential analysis identified the top 500 most up- or down-regulated genes between UBD-high and UBD-low groups across different cancer types, which were used as a UBD-associated gene signature. An RData file containing 1,288 compound-related signatures was downloaded from the database website (https://www.pmgenomics.ca/bhklab/sites/default/files/downloads) for matching score calculation. The analytical procedure followed methods outlined in previous publications (1, 22, 23). The results across 31 cancer types were summarized and visualized using the pheatmap package in R.

The CB-Dock2 web server (https://cadd.labshare.cn/cb-dock2), which employs AutoDock Vina’s algorithm, was utilized for molecular docking studies. Default parameters were applied throughout the simulations. Drug molecular structures were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/), and the protein target was downloaded from RCSB PDB (https://www.rcsb.org/).

MD simulations were performed using GROMACS 2022. Force field parameters were obtained using the pdb2gmx tool in GROMACS and the AutoFF web server. The CHARMM36 force field was applied to the receptor protein, while the CGenFF force field was used for the ligand molecules. The system was solvated with a cubic TIP3P water box with a margin of 1 nm around the system. Using the gmx genion tool, ions were added to achieve electrostatic neutrality of the system. Long-range electrostatic interactions were treated using the Particle Mesh Ewald (PME) method, with a cutoff distance of 1 nm. Constraints on all bonds were handled using the SHAKE algorithm, and the Verlet leap-frog algorithm was employed with an integration time step of 1 fs.

Prior to the MD simulation, energy minimization was carried out. This involved 3000 steps of steepest descent minimization followed by 2000 steps of conjugate gradient minimization. The energy minimization protocol included the following steps: first, constraining the solute while minimizing the energy of the water molecules; next, constraining the counterions and performing energy minimization; finally, performing unconstrained energy minimization on the entire system.

The MD simulations were conducted under NPT ensemble conditions at a temperature of 310 K and constant pressure, with a total simulation time of 50 ns. During the simulations, the tools g-rmsd, g-rmsf, g-hbond, g-Rg, and g-sasa were used to calculate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), hydrogen bonds (HBonds), radius of gyration (Rg), and solvent-accessible surface area (SASA), respectively.

All data were processed using R version 4.2.1. Pearson’s correlation coefficient was applied to normally distributed data, while Spearman rank correlation was used for non - normally distributed data. Comparisons between two variables were evaluated using the Wilcoxon signed - rank test and the Wilcoxon rank - sum test, respectively. The Kruskal - Wallis test was used to analyze variations among multiple variables. The diagnostic capability of UBD was evaluated using ROC analysis with the ‘pROC’ R package (24). Statistical significance was defined as a p - value less than 0.05, with high significance indicated by a p - value less than 0.0001 (denoted as *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001).

We utilized TCGA cohorts containing both normal and tumor tissue samples to assess the differential expression of UBD across various cancer types (Figure 1A). Our analysis revealed that UBD was significantly upregulated in 14 cancer types and downregulated in two malignancies (KICH: Kidney Chromophobe; THCA: Thyroid Carcinoma). In the TCGA cohort, analysis of paired tumor and adjacent normal tissues revealed that UBD was significantly overexpressed in 10 tumor types, while it remained significantly downregulated in THCA (Figure 1B).

ROC analyses were performed to explore the diagnostic potential of UBD in various types of tumors. The results indicated that the gene expression levels of UBD exhibited strong diagnostic efficacy for gastrointestinal tumors, whether assessed solely in the TCGA cohort or when combined with normal tissue samples from the GTEx database (Figures 2A, B).

Tumor samples from the TCGA cohort (n = 9,104) were stratified into high and low UBD expression groups based on the median expression value of UBD. These two groups exhibited significantly different immune subtypes (Figure 3A). Notably, the proportion of patients with tumors of the C2 (IFN-γ Dominant) subtype was twice as high in the high UBD expression group compared to the low UBD expression group (38% vs. 19%). Additionally, the proportion of patients with tumors of the C1 (Wound Healing) subtype was significantly lower in the high UBD expression group compared to the low UBD expression group (20% vs. 33%). After evaluating the expression levels of UBD in tumors with varying microsatellite instability (MSI) statuses, it was found that UBD expression was significantly associated only with the MSI status in colorectal adenocarcinoma (COAD). In COAD samples, UBD expression levels were significantly higher in those with high microsatellite instability (MSI-H) compared to those with low microsatellite instability (MSI-L) or microsatellite stability (MSS, Figure 3B).

Further investigation into the relationship between UBD expression and immune subtypes in COAD revealed a strong association between high UBD expression and the C2 (IFN-γ Dominant) subtype. Specifically, the proportion of the C2 (IFN-γ Dominant) subtype was significantly higher in COAD tumors with high UBD expression compared to those with low UBD expression (32% vs. 7%, Supplementary Figure 1).

The z-score algorithm proposed by Lee et al. was employed to explore the relationships between UBD and various tumor-related pathways across different cancer types. Previous studies have indicated that the Inflammation and Apoptosis pathways are closely related to the C2 (IFN-γ Dominant) subtype. In this study, the expression levels of UBD exhibited the highest positive correlation with the Inflammation pathway (R = 0.52, Figure 4A). Additionally, UBD expression levels also showed a significant positive correlation with the Apoptosis pathway (R = 0.38). We further evaluated the correlations between UBD expression and the Inflammation and Apoptosis pathways across different tumor types (Figures 4B, C). In more than two-thirds of the tumor types, UBD expression was significantly positively correlated with the Inflammation pathway, with correlation coefficients exceeding 0.5. Among these, the highest positive correlation was observed in UVM, with a correlation coefficient of R = 0.81. In KICH, THCA, and UVM, the Apoptosis pathway was significantly positively correlated with UBD expression, with correlation coefficients exceeding 0.6.

Given the close association between the C2 (IFN-γ Dominant) subtype and UBD expression in COAD, we further investigated the relationships between UBD and various pathways in COAD (Figure 4D). As expected, UBD exhibited the highest positive correlation with the Inflammation pathway (R = 0.52). Additionally, UBD showed a significant positive correlation with the Apoptosis pathway (R = 0.34).

Consistent with the characteristics of the C2 (IFN-γ Dominant) subtype, several inflammation-related Hallmark pathways, including IL2-STAT5 Signaling, IL6-JAK-STAT3 Signaling, Inflammatory Response, Interferon Alpha Response, and Interferon Gamma Response, exhibited significant positive correlations with UBD expression levels (Supplementary Figure 2).

Additionally, the relationship between UBD mRNA levels and functional proteins expression measured by RPPA in the TCPA database was evaluated using Rank - based association analysis. In UVM, UBD exhibited strong positive correlations with JAK2, S6, CMET, and JAB1. Conversely, UBD showed significant negative correlations with c-MYC, NOTCH1, YAP, HER2 (phosphorylated at Y1248), and RICTOR (phosphorylated at T1135; Figure 5A). In testicular germ cell tumors (TGCT), UBD expression was significantly positively correlated with SYK, PI3K p85, IRF1, PKC-panbeta II (phosphorylated at S660), and STAT5α.

Additionally, we identified an intriguing protein cluster composed of CASPASE 7 (cleaved at D198), PDL1, LCK, PREX1, SYK, and PI3K p85. These proteins exhibited significant positive correlations with UBD mRNA expression across more than five cancer types. These proteins may be closely related to the biological functions of UBD.

Caspase-7 (cleaved at D198) is a key effector caspase in the apoptosis cascade, responsible for cleaving and inactivating cellular substrates to execute apoptosis. Caspase-7 (cleaved at D198) exhibited a significant positive correlation with UBD, further suggesting a potential close association between UBD and apoptosis (Figure 5B). Proteomic enrichment analysis from the TCPA database indicated that the Apoptosis pathway was in an “activated” state among patients with high UBD expression across 17 cancer types (Figure 5C).

By integrating various algorithms for calculating immune infiltration scores, we found that the infiltration levels of CD8+ T cells and M1 macrophages were consistently and significantly positively correlated with UBD gene expression across nearly all cancer types (Figure 6). The C2 (IFN-γ Dominant) subtype is characterized by a robust immune response with significant infiltration of CD8+ T cells and M1 macrophages. These findings further substantiate the strong association between UBD and the C2 immune subtype of tumors.

Our study evaluated the prognostic value of UBD expression for OS, DFI, DSS, and PFI across various cancers using a univariate Cox regression model (Figures 7A–D). As shown in Figure 7A, elevated UBD expression was significantly associated with improved OS in breast cancer (BRCA; p=0.009, HR = 0.949) and exerted a protective effect in melanoma (SKCM; p<0.001, HR = 0.891) and sarcoma (SARC; p=0.006, HR = 0.899). Elevated UBD levels are associated with poorer OS in UVM (p<0.001, HR = 1.298), KIRP (p=0.016, HR = 1.155), and PAAD (p=0.012, HR = 1.143).

For DSS, elevated UBD levels were significantly associated with shorter DSS in KIRP (p=0.002, HR = 1.257), thymoma (THYM; p=0.013, HR = 2.732), and UVM (p=0.002, HR = 1.268). Conversely, UBD expression exerted a protective effect in bladder cancer (BLCA; p=0.012, HR = 0.927) and SKCM (p<0.001, HR = 0.884; Figure 7B).

Analysis of DFI revealed that high UBD expression was associated with a higher DFI rate in BLCA (p=0.007, HR = 0.856; Figure 7C). Additionally, elevated UBD expression correlated with improved PFI in SKCM (p=0.001, HR = 1.323), but with worse PFI in THYM (p=0.001, HR = 1.323) and UVM (p=0.009, HR = 1.206; Figure 7D).

To further validate the prognostic utility of UBD, we examined additional datasets containing prognostic information and found consistent results (Supplementary Figure 3).

Using the R package maxstat to determine the optimal cutoff values for UBD, KM analysis revealed that UBD is a prognostic factor for OS in KIRP, UVM, ESCA, THYM, PAAD, and LGG, acting as a protective factor. In contrast, UBD is associated with poorer prognosis in OV, CESC, UCEC, BRCA, BLCA, SKCM, and SARC (Figure 8A). The prognostic value of UBD across pan-cancer was visualized using a heatmap (Figure 8B).

In most cancer types, the methylation levels of UBD across different DNA methylation regions are downregulated in tumors (Figure 9A). We visualized the correlation between methylation levels in different regions and UBD gene expression levels across various cancer types using a bubble chart (Figure 9B). In BLCA, the mean β values of the TSS1500 (R = 0.36), shelf (R = 0.51), and shore (R = 0.58) methylation regions were significantly positively correlated with UBD gene expression levels (Figure 9C). In BRCA, the mean β values of the shelf (R = 0.52) and shore (R = 0.58) methylation regions also exhibited significant positive correlations with UBD gene expression levels. In UCEC, the mean β values of the shelf (R = 0.48), and shore (R = 0.55) methylation regions were significantly positively correlated with UBD gene expression levels. In KIRP, the mean β values of the opensea (R = -0.47), Promoter (R = -0.47), and TSS200 (R = -0.51) methylation regions were significantly positively correlated with UBD gene expression levels.

Using XSum, CMap can significantly enrich true positive drug-indication pairs through a novel matching algorithm. The lower the XSum relative score, the more likely the drug is to exert anti-tumor effects. When samples were divided into high and low UBD expression groups based on the median value, the results showed that two drugs (imatinib and TTNPB) had therapeutic potential in the high UBD expression group across more than 20 cancer types (Figure 10). Imatinib, a tyrosine kinase inhibitor, is commonly used for the treatment of chronic myeloid leukemia (CML) and gastrointestinal stromal tumors (GIST). TTNPB, a retinoic acid receptor (RAR) agonist, has the ability to bind to nuclear RARs with high affinity and can induce G1 cell cycle arrest.

We performed molecular docking of the UBD protein with Imatinib (Figure 11A) and TTNPB (Figure 11B). We present the five lowest-energy docking poses for each complex, as determined by Vina scoring. The Vina scores for all five docking poses between the UBD protein and Imatinib were lower than –7 kcal/mol, with the C3 pose having the lowest score of –8.9 kcal/mol (Figure 11C). Therefore, compared to TTNPB, Imatinib shows better docking affinity with UBD. Accordingly, we conducted MD simulations on the molecular docking complex of Imatinib and UBD.

The root mean square deviation (RMSD) is a good indicator of the conformational stability of proteins and ligands, as well as a measure of the deviation of atomic positions from their initial positions. A smaller deviation indicates greater conformational stability. Therefore, RMSD was used to evaluate the equilibration of the simulation system. As shown in Figure 11D, the complex system reached equilibrium after 10 ns and eventually fluctuated around 3.7 Å. Thus, Imatinib exhibits high binding stability with UBD. The radius of gyration (Rg) is a measure that describes overall structural changes and can be used to characterize the compactness of a protein structure. A greater variation in Rg indicates a more expanded system. Further analysis revealed that the complex system exhibited slight fluctuations during the simulation and gradually stabilized. This suggests that conformational changes occurred in the Imatinib–UBD complex during the simulation (Figure 11E). The solvent-accessible surface area (SASA) is a metric used to evaluate the surface area of a protein. In this simulation, the SASA between Imatinib and UBD was calculated (Figure 11F). The results show that the complex system exhibits slight fluctuations and gradually stabilizes over time. This demonstrates that the binding of the small molecule affects the local microenvironment and leads to a certain degree of change in SASA. Hydrogen bonds play an important role in the binding of ligands to proteins. The number of hydrogen bonds between Imatinib and UBD during the molecular dynamics simulation is shown in Figure 11G. The number of hydrogen bonds ranged from 0 to 6, and in most cases, the complex formed approximately 5 hydrogen bonds, indicating strong and favorable hydrogen bonding interactions between the ligand and the target protein. The root mean square fluctuation (RMSF) indicates the flexibility of amino acid residues within a protein. As shown in Figure 11H, the RMSF values for this complex are relatively low (mostly below 0.9-2.7 Å), indicating that the residues have lower flexibility and higher stability. The free energy landscape (FEL) illustrates the free energy distribution calculated based on RMSD and RG during molecular dynamics simulations of protein-ligand interactions. Color gradients are employed to represent free energy levels, transitioning from red (high energy) to blue (low energy). The dynamic simulation process is depicted in Figure 11I.

In summary, the complex exhibits stable binding and has favorable hydrogen bonding interactions. Therefore, the binding interaction between Imatinib and UBD is strong and effective.

RT-qPCR analysis revealed a significant upregulation of UBD gene expression in ESCA tumor tissues, with expression levels positively correlating with tumor staging (Supplementary Figure 4A). The highest UBD expression was observed in stage IV ESCA. WB analysis demonstrated upregulated UBD protein expression levels in ESCA tumor samples, with expression levels progressively increasing in correlation with advancing tumor stages (Supplementary Figures 4B, C). IF assays further confirmed elevated UBD expression in ESCA tumor tissues (Supplementary Figure 4D).

A total of 30 paraffin-embedded tissue blocks from ESCA patients were included in this study. The clinical characteristics of the patients are summarized in Supplementary Table 1. Based on the IHC-based AOD values of UBD, the 30 enrolled ESCA patients were stratified into high, moderate, and low UBD groups using tertiles (Supplementary Figure 4E). The low UBD group exhibited the highest proportion of stage I (40%) and stage II ESCA (50%) cases and the lowest proportion of stage III cases (10%), while all stage IV patients were exclusively categorized into the high UBD group. OS was defined as the interval from the date of initial diagnosis to the occurrence of all-cause mortality or the last follow-up contact. Survival data were primarily obtained through the hospital’s electronic health records system. In instances where death documentation was unavailable, survival status verification was performed via structured telephone interviews. Kaplan-Meier curve analysis demonstrated that patients in the high UBD group exhibited significantly poorer OS outcomes (Supplementary Figure 4F).

Previous results indicated that esophageal cancer patients with high UBD expression had poorer prognoses, prompting us to further explore the biological functions of UBD in esophageal cancer. First, we established a stable TE-11 cell line overexpressing UBD. Both UBD mRNA and protein levels were upregulated by more than 5-fold (Supplementary Figures 5A, B). Plate colony formation assay (Supplementary Figure 5C) and CCK8 cell proliferation curves (Supplementary Figure 5D) both indicated that the proliferative capacity of TE-11 cells overexpressing UBD was significantly enhanced. Wound healing assay (Supplementary Figure 5E) and Transwell migration assay (Supplementary Figure 5F) both demonstrated that the migratory ability of TE-11 cells overexpressing UBD was significantly increased.

To investigate the underlying mechanisms responsible for these observed differences in proliferation and migration, we performed transcriptome sequencing on TE-11 cells transduced with empty vector control lentivirus (n=3) and TE-11 cells overexpressing UBD (n=3).

Using thresholds of absolute fold change >2 and false discovery rate (FDR) <0.05, volcano plots showed that 221 genes were upregulated and 328 genes were downregulated in UBD-overexpressing TE-11 cells (Supplementary Figures 5G, H). The heatmap displays the top 20 upregulated and downregulated genes (Supplementary Figure 5I). Pathway enrichment analysis using all differentially expressed genes was performed with Metascape, and the results indicated that the TP53 signaling pathway and TP53-related pathways (such as cell cycle regulation and cellular senescence) were significantly enriched in TE-11 cells overexpressing UBD (Supplementary Figure 5J). To further validate our findings, we used Western blotting to examine key proteins related to TP53 and the cell cycle, including c-myc, cyclin B1, CDK1, CDK4, p53, and p21. The expression levels of c-myc, cyclin B1, CDK1, and CDK4 were upregulated in TE-11 cells overexpressing UBD, while the expression levels of p53 and p21 were significantly downregulated (Supplementary Figure 5K). These results are consistent with our enrichment analysis. Therefore, the promotion of esophageal cancer proliferation and migration by UBD may be mediated through the TP53 signaling pathway.