Comprehensive analysis based on the ubiquitination- and deubiquitylation-related genes reveals the function of NEURL3 in esophageal squamous cell carcinoma

The TCGA database TCGA-ESCC (https://portal.gdc.cancer.gov↗) was used to download the gene expression data (TPM) of 95 ESCC samples with corresponding patient’s clinical information. In addition, we obtained two datasets GSE53624↗ and GSE53622↗ from the GEO database (https://www.ncbi.nlm.nih.gov/geo↗). Among them, the GSE53624↗ dataset contained gene expression profiles of 119 ESCC samples and 119 matched normal tissues with complete survival information, and it was mainly used for preliminary analysis of this study. For external validation of prognostic features, we used the TCGA-ESCC dataset and GSE53622↗ dataset including microarray data from 60 ESCC samples.

GeneCards (http://www.genecards.org/↗) (17) is a searchable, integrative database that provides comprehensive, user-friendly information on all annotated and predicted human genes. Utilizing this platform, we filtered for protein-coding URGs and DRGs with a relevance score above 7.00 using the keyword “Ubiquitination” and “Deubiquitylation”. By combining these genes, we obtained a total of 977 URGs and DRGs. The gene symbols were listed in Supplementary Data 1. The research flowchart was shown in Supplementary Figure S1.

We analyzed the GSE53624↗ dataset using the limma R package (https://cran.r-project.org/package=limma↗) to identify differentially expressed genes (DEGs) with | log2 fold change | > 0.7 and P < 0.05. Meanwhile, the univariate cox proportional hazards regression analysis was conducted to determine the potential genes associated with prognosis of ESCC patient (P < 0.05) using survival R package (https://cran.r-project.org/package=survival↗). The gene symbols were listed in Supplementary Data 1. Subsequently, a total of 13 URGs and DRGs that were both DEGs and potential prognosis genes were identified using the Venn diagram. Then, we employed the LASSO regression analyses using glmnet R package (https://cran.r-project.org/package=glmnet↗) to further compress the number of genes and finally determined 11 URGs and DRGs to form the prognosis signature.

The risk score of each patient based on the URGs and DRGs prognosis signature was calculated by performing the survival and nomogramFormula R packsge(https://cran.r-project.org/package=nomogramFormula↗), and the ggrisk R package (https://cran.r-project.org/package=ggrisk↗) was employed to illustrate the association between risk curves and patient survival distributions. Meanwhile, the best cutoff value of the risk score was determined by cutoff R package (https://cran.r-project.org/package=cutoff↗). Then, we used the best cutoff value to classify ESCC patients into two groups (low- and high-risk groups). Subsequently, the OS time of two groups was compared by the Kaplan-Meier analysis and log-rank test through the survminer R package. The URGs and DRGs prognosis signature’s discrimination ability was evaluated using the concordance index (C-index) and Receiver Operating Characteristics (ROC) analysis by the survival, survminer (https://cran.r-project.org/package=survminer↗) and timeROC R packages (https://cran.r-project.org/package=timeROC↗).

We performed unsupervised clustering analysis using the ConsensusClusterPlus R package (https://www.bioconductor.org/packages/release/bioc/html/ConsensusClusterPlus.html↗) to develop the molecular subtypes based on the 11 URGs and DRGs. The optimal number of clusters was determined by the area under the cumulative distribution function (CDF) curve, the descending trend of CDF delta and the consistency of sample clustering, and the optimal number of clusters selected was 2. Then, we conducted the Kaplan-Meier analysis to compare the OS time of patient in two clusters using the survminer R package. Additionally, we plotted heatmap to explain the relationship between the clusters and clinical feathers of patient using the pheatmap R package. The GSVA R package (https://www.bioconductor.org/packages/release/bioc/html/GSVA.html↗) was used for Gene Set Variation Analysis (GSVA). Subsequently, the limma R package was applied to assess the difference of biological function and signaling pathway between two subtypes (adjusted P < 0.05), and the top 30 different functions and pathways were visualized using the pheatmap R package. Meanwhile, through the ggalluvial R package, we plotted Sankey plots for the relationship between patients in the low- and high-risk groups and patients with CDF subtypes.

We performed Gene Ontology (GO) functional enrichment analysis through the David database (https://davidbioinformatics.nih.gov/↗), and then revealed the potential biological process, cellular component, and molecular function of the 11 prognostic URGs and DRGs. Furthermore, the underlying pathway associated with the DEGs in the high- and low-risk groups determined by the URGs and DRGs prognosis signature was enriched through the GO, Disease Ontology (DO), and KEGG analysis using BioEnricher R package (https://github.com/Zaoqu-Liu/BioEnricher↗). In addition, several biochemical processes of the DEGs were inferred by calculating the Kyoto Encyclopedia of Genes and Genomes (KEGG) official website signaling pathway-related gene set score based on the single-sample gene set enrichment analysis (ssGSEA) algorithm. Additionally, we also calculated and compared the Epithelial-mesenchymal transition (EMT) score and angiogenesis score of high- and low-risk groups based on the ssGSEA algorithm.

This study employed the ssGSEA algorithm to accurately calculate the infiltration abundance of 28 immune cell types in each ESCC patient. Then, the correlation between risk score and immune cell’s abundance was analyzed by spearman correlation analysis. Moreover, using spearman correlation analysis, we examined the correlation between the risk score and gene markers of immune cells (18), or immune checkpoint-related markers to determine the suitability of hub genes for predicting the efficacy of immunotherapy (19). Additional, to examine the variation in immune cell infiltration between high-risk and low-risk groups, we used the TIMER, CIBERSORT-ABS, QUANTISEQ, EPIC, MCPCOUNTER, CIBERSORT and XCELL algorithms by IOBR R package (https://github.com/IOBR/IOBR↗).

The oncoPredict R packages (https://cran.r-project.org/package=oncoPredict↗) were used to assess the differences in drug sensitivity of chemotherapeutic drugs across different risk groups based on Genomics of Drug Sensitivity in Cancer (GDSC) database which provided the expression profiles and IC50 for anti-tumor drugs in more than 1000 cancer cells (20). The correlation between drug sensitivity and the risk score, as well as the expression of the 11 prognostic URGs and DRGs were calculated by Spearman correlation analysis. Additionally, using the AMDock software, the 3D structure of dacomitinib and talazoparib from PubChem database and the MAPK14 (IDENTIFIER = 1A9U) downloaded from RCSB PDB database, we performed computer molecular docking simulations to further verify the results of our drug screening.

The URGs and DRGs prognosis signature and clinical parameters associated with OS were selected by univariate and multivariate cox proportional hazards regression analysis. Next, A nomogram with endpoints of 1-, 3-, and 5-year OS was constructed based on these independent prognosis factors. The C-index and ROC analysis were used to evaluate the discrimination ability of nomogram. Time-dependent C-index curves were plotted by the pec package (https://cran.r-project.org/package=pec↗). Calibration curve was used to evaluate the average prediction probability of the nomogram by comparing observed outcomes with predicted probabilities, and decision curve analysis (DCA) was conducted to estimate the clinical utility of the nomogram. The Kaplan-Meier method and log-rank test were applied for comparing the variations in OS time of different risk groups.

KYSE30 and KYSE150 cells were maintained in RPMI-1640 medium (Thermo Fisher Scientific) supplemented with 10% fetal bovine serum (GIBCO) and 1% penicillin/streptomycin (C100C5, NCM Biotech). Cells were routinely cultured at 37°C in a humidified atmosphere containing 5% CO₂. The single siRNA oligonucleotides targeting human NEURL3 (GenePharma) and negative control siRNA were transfected by using Lipofectamine™ RNAiMAX (13778150, Life Technologies).

Cells were lysed in RIPA lysis buffer (R0010, Solarbio) supplemented with Protease and Phosphatase Inhibitor Cocktail (P002, NCM Biotech). The protein sample was separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to polyvinylidene difluoride membranes (FFP33, Beyotime Biotechnology). After being blocked with 5% Bovine Serum Albumin V solution (A8020, Solarbio), the membranes were incubated overnight at 4°C with the following primary antibodies: anti-NEURL3 (1:1000, A20479, ABclonal) and anti-GAPDH (1:5000, HRP-60004, Proteintech). Then, the membranes were incubated with goat anti-rabbit antibody (7076, 7074, Cell Signaling Technology) at room temperature for 1 hour. Protein bands were visualized using Western blotting luminol reagent (sc-2048, Santa Cruz).

Total RNA was extracted from cells using the TRIzol reagent (15596018, Invitrogen). Then, the extracted RNA was reverse transcribed into cDNA with the HiScript III RT SuperMix for qPCR (+gDNA wiper) (R323, Vazyme). qPCR was performed using the ChamQ Universal SYBR qPCR Master Mix (Q711, Vazyme). The relative expression of transcripts was analyzed using the 2-ΔΔCt method. β-actin served as the internal control. All assays were performed according to the manufacturer’s instructions. The primers used in this method were displayed as following: NEURL3, F: 5’-ATGGGACCTCAGCAACAAGGCT-3’, R: 5’-AAGACCCGCCAGGCACAGTATC-3’; β-actin, F: 5’-CAACTGGGACGACATGGAGAAA-3’, R: 5’-GATAGCAACGTACATGGCTGGG-3’.

Twenty-four hours after transfection, the cells number was countered, and 1000 cells were seeded into 96-well plates. At specific time points, 10μl of Cell Counting Kit (CCK-8; C0039, Beyotime Biotechnology) was added to each well, and the absorbance at 450 nm was measured 2 hours later using a spectrophotometer.

ESCC cells were seeded into 6-well plates and then were incubated in serum-free medium for 12 hours. The sterile 200μl tip was used to draw a straight line across the cell culture, and the medium was replaced with serum-free medium. Images were acquired at specific time points using the inverted microscope, and the wound healing area was measured using ImageJ software.

Transwell chambers equipped with 8.0μm pore inserts (353097, CORNING) or with a matrigel coating (354480, CORNING) were employed. Briefly, ESCC cells were suspended in serum-free medium and seeded into the upper chamber at a density of 5×10⁴ cells (migration assay) or 1×10⁵ cells (invasion assay) per well. The lower chamber was filled with medium supplemented with 10% FBS. After 48 or 72 hours of incubation, cells were fixed with 4% paraformaldehyde, stained with 0.5% crystal violet, and photographed under an inverted microscope.

All statistical analyses were performed using the following softwires: IBM SPSS Statistics 24, Microsoft Excel, X-tile version 3.6.1 and R version 4.2.2. The independent t-test was used for normally distributed continuous variables, while the Mann-Whitney test was employed to contrast non-normally distributed variables. Categorical variables were compared with the chi-square test or the Fisher exact test. P < 0.05 was considered statistically significant.

Using the limma package, 5031 DEGs were found between the ESCC and normal groups in the training dataset GSE53624↗ based on the criteria of |log2 fold change| > 0.7 and P < 0.05 (Figure 1A). Meanwhile, 1709 genes associated with prognosis (P < 0.05) in GSE53624↗ were identified by the univariate cox proportional hazards regression analysis. The gene symbols were listed in Supplementary Data 1. A total of 13 URGs and DRGs that were both DEGs and associated with prognosis were obtained using Venn diagram (Figure 1B, Supplementary Figure S2A). Next, we conducted the LASSO analysis to remove the overfitting genes, and the URGs and DRGs prognostic signature consisting of 11 genes was created (Figures 1C, D, Supplementary Figure S2B). As shown in Figures 1E, F, there were 5 genes (CCNF, APP, NEURL3, CCNE1, AURKA) up-regulated, and 6 genes (FBXW5, TRIM13, NBR1, ANXA11, MAPK14, and CFTR) down-regulated in the ESCC group when compared to the normal group in GSE53624↗. The expression distribution of these 11 genes was plotted in the ESCC group of GSE53624↗ (Supplementary Figure S2C).

We intended to develop a new predictive prognosis signature for ESCC based on the 11 prognostic URGs and DRGs previously identified. In training dataset GSE53624↗, the risk score of each patient based on the URGs and DRGs prognosis signature was calculated, and patients were classified into low- or high-risk groups as per the best cutoff value of the risk score (Figure 2A). The alterations in the expression of 11 prognostic URGs and DRGs between the two risk groups were illustrated in (Supplementary Figure S3A). Compared with patients in the low-risk group, patients in the high-risk group had worse prognosis (Figure 2B). Moreover, the ROC analysis was performed to evaluate the URGs and DRGs prognosis signature, which demonstrated that the AUC values for 1-, 3- and 5-year were 0.875, 0.819, and 0.805, respectively (Figure 2C, Supplementary Figure S3D). In addition, the C-index of the URGs and DRGs prognosis signature in GSE53624↗ was 0.699 (95% CI: 0.677-0.721) (Supplementary Figure S3E).

Furthermore, the predictive ability of the URGs and DRGs prognosis signature in ESCC was validated in the external validation sets (GSE53622↗ and TCGA-ESCC dataset). This ensures the validity and general applicability of the prognostic signature. In line with aforementioned analytical approach, the patients with ESCC in the validation datasets were also classified into low- or high-risk groups (Figures 2D, G). The expression of 11 prognostic URGs and DRGs changed between the high and low-risk groups in the validation datasets, as depicted in (Supplementary Figures S3B, C). The Kaplan-Meier analysis revealed that the low-risk individuals were consistently associated with more favorable prognosis in contrast to those at high-risk group (Figures 2E, H). Moreover, the URGs and DRGs prognosis signature’s AUCs of 1-, 2- and 3-year periods were 0.738, 0.756, and 0.762 in GSE53622↗, correspondingly, 0.631, 0.830, and 0.806 in TCGA-ESCC dataset (Figures 2F, I, Supplementary Figure S3D). The C-index of the URGs and DRGs prognosis signature were 0.685 (0.647-0.723) in GSE53622↗ and 0.639 (0.609-0.669) in TCGA-ESCC dataset (Supplementary Figure S3E). In summary, the above results indicated that the prognostic signature based on the 11 URGs and DRGs possessed remarkable discriminative ability and predictive accuracy for prognosis prediction for ESCC patient.

Combined with the area under the CDF curve, the descending trend of CDF delta and the consistency of sample clustering, we successfully identified two subtypes based on the 11 URGs and DRGs (k = 2, C1 = 61, C2 = 58) (Figures 3A–D). Notably, through comparison of survival time of patient in two consensus subtypes, we observed that patients in C1 exhibited significantly better OS compared to those in C2 (P = 0.05) (Figure 3E), further confirming the prognosis value of these 11 URGs and DRGs. Additionally, we also visualized the distribution of the two subtypes of samples on clinical data in a heatmap (Supplementary Figure S4A). Furthermore, we utilized GSVA to calculate the functional states and signaling pathway scores of the two subtypes in GO and KEGG catalogues. Interestingly, we found that, besides the pathways related to ubiquitination, DNA replication, damage and repair pathways and immune signaling pathways were also enriched in C2 (Supplementary Figures S4B, C). This finding implied the reason for poor prognosis of patient in C2. Finally, we revealed the patients had a worse prognosis in both the C2 subtype and high-risk group from the result of the Sankey analysis (Figure 3F).

The GO functional enrichment analysis was conducted to investigate the biological process, cellular component, and molecular function of the 11 prognostic URGs and DRGs. We noticed that these 11 genes were involved in the cell cycle, immune response, and apoptotic process in the biological process category; cytosol, cytoplasm, and early endosome in the cellular component category; enzyme binding, cyclin-dependent protein serine/threonine kinase regulator activity, and protein kinase binding in the molecular function category (). 1

The GO, KEGG and DO enrichment analyses were further conducted to investigate the latent biological roles of the DEGs between the high- and low-risk groups determined by the URGs and DRGs prognosis signature. Accordingly, the GO result revealed that the DEGs were prominently enriched in vasculature development, extracellular structure organization, and extracellular matrix organization (Figure 4A). The KEGG result suggested that the DEGs were primarily enriched in focal adhesion, cytokine-cytokine receptor interaction, and rheumatoid arthritis (Figure 4B). The DO result showed that the DEGs were connected with arthritis, bone inflammation disease, and lung disease (Figure 4C). Additionally, we calculated the KEGG official website signaling pathway-related gene set score based on the ssGSEA algorithm. The result revealed that the DEGs are significantly enriched in various signaling pathways including ECM-receptor interaction, MAPK signaling pathway, PI3K-Akt signaling pathway, TGF-beta signaling pathway, TNF signaling pathway and Wnt signaling pathway (Figure 4D). Furthermore, we also calculated the EMT score and angiogenesis score of two risk groups, and both of scores were all higher in high-risk group in contrast with that of the low-risk group (Figures 4E, F).

Given the significant impact of tumor immune microenvironment on therapeutic outcomes and patient prognosis in malignant tumors, we further investigated the potential association between risk score of patient and immune cell’s infiltration abundance. We conducted comprehensive analysis of 28 immune cell types using the ssGSEA algorithm, and the results revealed significant positive correlations between risk score and Effector memory CD4⁺ T cell, Memory B cell, Natural killer cell, Natural killer T cell, and Type 2 T helper cell (P < 0.001, R > 0.3) (Figure 5A, Supplementary Figures S6A–E). We further analyzed the link between risk score and the gene markers of immune cells or the immune checkpoint-related markers. The results showed a close association between risk score and CSF1R of monocytes, CCL2 of tumor-associated macrophages, CD163, VSIG4, and MS4A4A of M2 macrophages, ITGAM and CEACAM8 of neutrophils, KIR2DL3 and KIR3DL2 of natural killer cells, and HLA-DPB1, HLA-DRA, HLA-DPA1, and NRP1 of dendritic cells (Figure 5B), as well as closed associated with immune checkpoint markers such as NRP1, LAIR1, and TNFSF4 (Figure 5C). Additionally, the distinctions of immune cell levels between high- and low-risk groups were also explored through CIBERSORT, MCPCOUNTER, QUANTISEQ, EPIC, TIMER, CIBERSORT-ABS, and XCELL algorithms. As per the findings, the low-risk group had remarkably higher levels in most immune cells, including Plasma cells, regulatory T cells, Monocytes, Basophils, Memory B cells, Myocytes, pro B cells, B cells, and CD4⁺ T cells (Figure 5D). This may be the reason that the low-risk group has a superior prognosis.

To further understand the clinical application of the URGs and DRGs prognosis signature, we evaluated the association between drug sensitivity and the risk score of ESCC patients. A total of 87 drugs were negatively correlated with the risk score and were regarded as sensitive drugs for high-risk groups, while only 7 drugs showed significantly positive correlation that considered to be drug resistance to high-risk patient but might be appropriate treatment for low-risk patients (Figure 6A). Then, we further assessed the correlation of talazoparib and dacomitinib, of which showed the strongest negative correlation or the most prominent positive correlation, with the 11 prognostic URGs and DRGs, respectively. The results showed that the expression level of most URGs and DRGs were significantly correlated with the sensitivity to talazoparib and dacomitinib (Figures 6B, C, Supplementary Figures S7A, B). Additionally, the difference of estimated IC50 in these two drugs across different risk groups corroborated the findings of our correlation analysis (Supplementary Figure S7C). According to the results, as MAPK14 was most negatively correlated with talazoparib and positively correlated with dacomitinib, we selected it for molecular docking with dacomitinib and talazoparib, respectively. Based on the possible results analyzed by AMDock software, we found that both of these two drugs showed excellent binding ability with MAPK14 (Figures 6D, E), suggesting that they may serve as anti-ESCC drugs through targeting MAPK14.

To further complement the URGs and DRGs prognosis signature’s predictive efficiency, we conducted the univariate and multivariate Cox regression analyses in training dataset GSE53624↗. We proved that the URGs and DRGs prognosis signature, TNM stage and age independently acted as robust prognostic markers (P < 0.05) (Table 1, Supplementary Figures S8A, B), which were integrated to construct for 1-, 3-, and 5-year OS nomogram (Figure 7A). The calibration curve analysis suggested that the predicted outcomes of the nomogram were in excellent agreement with the actual OS of ESCC patients (Figure 7B). The DCA curve analysis confirmed the better clinical utility of the nomogram (Figure 7C). The C-index of nomogram was 0.744 (95% CI: 0.721-0.767), which was higher than that of the URGs and DRGs prognosis signature (0.699, 95% CI: 0.677-0.721) (Supplementary Figure S8H). The time-dependent C-index analysis also showed that the nomogram had better discrimination ability in predicting survival of patients (Figure 7D). Moreover, the ROC curve analysis manifested that the AUC values for 1-, 3- and 5-year were 0.840, 0.820, and 0.818, respectively (Figure 7E). In addition, we calculated the total point of each patient and obtained the cutoff value by cutoff package. Subsequently, the ESCC patients were subdivided into low- or high-risk groups based on the cutoff value, and Kaplan-Meier analysis revealed that patients in the low-risk group exhibited a more favorable prognosis compared with high-risk group (Figure 7F), indicating that the nomogram can effectively determine the outlook for ESCC sufferers.

Meanwhile, we also conducted a series of analyses mentioned above in the validation groups. As shown in Supplementary Figures S8C, D, the nomogram-predicted probabilities of 1-, 2-, and 3-year OS were in consistency with the actual OS, and the more outstanding clinical practical value of the nomogram were also illustrated (Supplementary Figures S8E, F). The C-index and time-dependent C-index analyses still demonstrated that the nomogram gave good discrimination (Figures 8A, B, Supplementary Figure S8H). The AUC values for 1-, 2- and 3-year of the nomogram were 0.756, 0.806, and 0.819, correspondingly, in GSE53622↗ (Figure 8C), and 0.653, 0.765, and 0.863 in TCGA-ESCC dataset (Figure 8D; Supplementary Figure S8G). Moreover, high-risk ESCC patients were consistently associated with poor outcomes in both validation groups (Figures 8E, F). These results revealed that this novel nomogram possessed tolerable accuracy and applicability for survival prediction for ESCC patients.

Among the 11 prognostic DE-URGs and DRGs, we found that NEURL3 was the upregulated expression gene with the greatest multiple of expression difference, and it was also significantly highly expressed in ESCC tissues of GSE45670↗, GSE53622↗, GSE53624↗ and GSE161533↗ (Supplementary Figure S9A). However, the role of NEURL3 remains unclear in ESCC. To explore the biological function of NEURL3 in ESCC, firstly, we detected the expression of NEURL3 in several ESCC cell lines including KYSE30, KYSE140, KYSE150 and KYSE510 by using RT-qPCR and Western blotting, and the results showed that NEURL3 was highly expressed in both the KYSE30 and KYSE150 cells (Supplementary Figures S9B, C). Then, we knocked down NEURL3 using two siRNAs in these two cell lines. Figure 9A and Supplementary Figure S9D showed that NEURL3 could be effectively silenced by two independent siRNAs. The CCK8 assay showed that NEURL3 depletion inhibited proliferation ability of KYSE30 and KYSE150 cells (Figure 9B). The transwell assay showed that NEURL3 knockdown reduced the migration and invasion ability of KYSE30, and the results were confirmed in KYSE150 cells (Figure 9C). Subsequently, in the wound-healing assay, we found that NEURL3 depletion inhibited wound closure speed in both KYSE30 and KYSE150 cells (Figure 9D). These results indicated that NEURL3 could promote ESCC cell proliferation and motility in vitro.