Leveraging single-cell and multi-omics approaches to identify MTOR-centered deubiquitination signatures in esophageal cancer therapy

Gene expression RNAseq data and somatic mutation profiles in Mutation Annotation Format (MAF) were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/↗). In the TCGA dataset, gene expression profiles were quantified using Transcripts Per Million (TPM) estimates and then log2-transformed for further analysis. GEO datasets were obtained from the Gene Expression Omnibus (GEO) of the National Center for Biotechnology Information (NCBI), which is publicly accessible at https://www.ncbi.nlm.nih.gov/geo↗.

After batch effects were removed from the TCGA-ESCC and GSE53624↗ datasets using the empirical Bayesian-based R package sva, we merged the data. Deubiquitination-related genes were identified as prognostic factors using univariate Cox regression, Lasso regression, and multivariate Cox regression analyses. Based on the results of the multivariate Cox regression, a mathematical formula was developed to predict the risk score for each ESCC patient:

Subsequently, all patients were divided into high-risk and low-risk groups according to the median risk score for further analysis.

Based on gene expression profiles in ESCC, differential expression analysis was performed between high-risk and low-risk groups using the “limma” package in R to identify differentially expressed genes (DEGs). The thresholds for DEGs selection were defined as “adjusted P < 0.01 and |log(FoldChange)| > 0.5”. GO enrichment analysis was conducted using the “clusterProfiler” package in R, with a P-value less than 0.05 considered statistically significant. Data visualization was performed using the “ggplot2” package in R (21, 22). GO terms were categorized into three main ontology categories: Biological Process, Cellular Component, and Molecular Function.

Using multiple algorithms, including CIBERSORT, MCPcounter, TIMER, and xCell, we estimated the differences in the abundance of various immune cell types between high-risk and low-risk groups (23). These differences were visualized using a heatmap. Additionally, we performed a Spearman correlation analysis (r) to explore the relationship between the ESCC patient risk score and the StromalScore, ESTIMATEScore, ImmuneScore, and TumorPurity, with the results displayed in scatter plots.

The R package “Maftools” was used to analyze and visualize the gene mutation profiles and frequencies in high-risk and low-risk groups, as well as to compare the differences in ESTIMATE scores, stromal scores, immune scores, and tumor purity between these groups (24). Tumor mutation burden (TMB) was calculated as the number of mutations per million bases (mut/Mb), and an oncoplot was generated using the “oncoplot” package to create a waterfall plot. Patients were divided into four groups based on TMB (high or low) and risk score (high or low), and Kaplan-Meier survival curves were plotted to assess the impact of TMB and risk score on overall patient survival.

Single-cell sequencing data for esophageal cancer was downloaded from the GEO database (GSE160269↗) and preprocessed using the “Seurat-Req” package (25) for single-cell RNA sequencing (scRNA-seq) analysis (26–28). The “PercentageFeatureSet” function was used to assess the proportion of mitochondrial genes within the dataset. To ensure data quality and completeness, only genes expressed in at least three cells were retained. The scRNA-seq data was normalized using the “NormalizeData” function. After normalization, the data was converted into a Seurat object, and the “FindVariableFeatures” function was used to identify the top 2,000 highly variable genes.

Subsequently, these highly variable genes were scaled and subjected to principal component analysis (PCA) using the “RunPCA” tool. Dimensionality reduction and visualization in two-dimensional space were achieved through the optimization of the Shared Nearest Neighbor (SNN) module and the t-distributed Stochastic Neighbor Embedding (t-SNE) clustering algorithm. Various cell subpopulations were annotated based on marker genes. We then analyzed the expression profiles of modeled DUBGs and classified the single-cell data into high- and low-risk groups based on DUBGs, evaluating the relationship between immune cell infiltration and risk scores.

Finally, “CellChat” was employed to assess differences in signaling patterns and interaction networks between cell types in the high- and low-risk groups (29, 30).

Gene expression and drug sensitivity data for the same samples were downloaded from the CellMiner website (https://discover.nci.nih.gov/cellminer/↗). Drug sensitivity data were recorded after clinical laboratory validation and FDA standard certification. Pearson correlation analysis was used to determine the correlation between the risk score model and drug sensitivity.

The human ESCC cell line (KYSE30) and the mouse ESCC cell line (mEC25) were obtained from the National Biomedical Experimental Cell Bank (Beijing, China) (31). The cells were cultured in 1640 medium (Gibco, USA)supplemented with 10% fetal bovine serum (Hyclone, USA)and 1% penicillin/streptomycin (Gibco, USA)in a humidified incubator with 5% CO₂ at 37°C.

The day before transfection, plate the cells to ensure they reach 50-70% confluence at the time of transfection. On the day of transfection, prepare the transfection complex as follows: first, dilute the siRNA (sequences shown in Table 1 for both human and mouse MTOR) and the transfection reagent separately in serum-free medium. Next, combine the diluted siRNA solution with the transfection reagent, gently mix, and incubate the mixture at room temperature for 15-20 minutes to allow the complex to form (GenePharma, Shanghai, China). Next, remove the culture medium from the cells, add the complex to the cells, and supplement with fresh culture medium. Incubate the cells at 37°C in a 5% CO₂ incubator for 24-72 hours. To assess the efficiency of siRNA-mediated knockdown, total RNA was extracted from the transfected cells using a standard RNA isolation kit (e.g., TRIzol reagent) according to the manufacturer’s instructions. The RNA was then reverse-transcribed into cDNA using a reverse transcription kit. Quantitative real-time PCR (qPCR) was performed to measure the mRNA levels of the target gene (MTOR) using specific primers. The relative expression levels were normalized to GAPDH as an internal control, and knockdown efficiency was calculated using the 2^^(-ΔΔCt) method.

Transfected cells (5×10³) were seeded into the wells of a 96-well plate and incubated for 6 hours to allow adhesion (considered the starting point of the experiment, 0 h). Measurements were taken every 24 hours by adding 10 μL of CCK-8 reagent (MedChemExpress, USA) and 90 μL of 1640 medium to each well, followed by a 2-hour incubation at 37°C. The absorbance (OD) at 450 nm was then measured using a microplate reader (Thermo, USA). After a total of four measurements, a growth curve was plotted based on the OD values.

The C57/BL6 mice used in the experiments were housed at the Xuzhou Medical University Animal Experiment Center. The animal ethics for this study were reviewed and approved by the Experimental Animal Ethics Committee of Xuzhou Medical University. All mice were kept in an SPF (Specific Pathogen-Free) environment with free access to food and water. All experimental procedures were conducted in strict accordance with ethical guidelines, ensuring that the design and execution of the experiments adhered to the highest standards of animal welfare and ethical considerations.

To establish a subcutaneous tumor model using the mEC25 mouse esophageal cancer cell line in 6-8 weeks-old C57BL/6 mice (SPF (Beijing)BIOTECHNOLOGY CO., LTD). Begin by culturing the mEC25 cells until they reach the logarithmic growth phase. Harvest the cells and resuspend them in sterile PBS at a concentration of 1x10^7 cells/ml. Using an insulin syringe, inject 100 µL of the cell suspension subcutaneously into the right flanks of each mouse. Every 2-3 days, measure the tumor’s length (L) and width (W) using calipers, and calculate the tumor volume using the formula V = (L x W²)/2.

Data analysis and visualization were conducted using R version 4.0.2. Categorical variables were compared using the chi-square test, while group differences were assessed using the Student’ t-test and the Wilcoxon rank-sum test. The correlation between two parameters was evaluated using Spearman correlation analysis. All statistical tests were two-sided, and a P-value of less than 0.05 was considered statistically significant.

The complete workflow is shown in Figure 1. To explore the potential of utilizing DRGS for clinical decision support in ESCC, we combined the TCGA-ESCC (n=86) and GSE53624↗ (n=119) datasets to develop a prognostic model for ESCC. Initially, batch effects between the two datasets were removed (Figures 2A–C), followed by univariate Cox analysis to identify DUBs significantly associated with overall survival (OS) in ESCC patients (Figure 2D). To address overfitting risks and refine gene selection for OS prediction, we performed LASSO regression analysis, and stepwise multivariate Cox analysis subsequently identified USP2, ITCH, ESR1, AXIN1, MTOR, USP22, TRAF2, USP37, AKT1, AR, OTUD6B, ZC3H12A, and SMAD3 as independent prognostic factors. These genes were used to construct the DRGS for ESCC patients. The risk score was calculated by summing the expression levels of individual genes, each weighted by their respective regression coefficients (Figure 2H). Patients were then divided into high- and low-risk groups based on the median risk score (Figures 2E–G). Kaplan-Meier survival analysis revealed that ESCC patients in the high-risk group had significantly shorter survival probabilities, with 1-, 3-, and 5-year AUCs of 0.67, 0.74, and 0.75, respectively, indicating a robust predictive capability of the DRGS for ESCC patients (Figures 2I–J).

Next, we explored the functional classification of DRGS in ESCC patients and their potential roles and impacts in cancer immunity. As shown in Figure 3A, the Gene Ontology (GO) analysis indicates that DRGS may influence ESCC progression by participating in the regulation of biological processes such as Pattern Specification Process, Epidermal Cell Differentiation, Epidermis Cell Differentiation, Response to Other Organism, Response to External Biotic Stimulus, and Antimicrobial Humoral Response, thereby affecting patient prognosis. Figure 3B also suggests that DRGS may be involved in processes such as the release of cancer antigens, the presentation and activation of cancer antigens, immune cell recruitment, and tumor killing.

To explore the molecular mechanisms driving the abnormal expression of these 14 DUBGs, we examined the differences in TMB between high-risk and low-risk groups of ESCC patients. The results indicate that the high-risk group has a higher mutation frequency in several genes (such as TP53, TTN, NFE2L2) compared to the low-risk group. Additionally, the high-risk group exhibits generally lower ESTIMATE Score, Stromal Score, and Immune Score (Figure 4A). The high-risk group also has a higher standardized Tumor Mutation Burden (TMB, Figure 4B), and there is a positive correlation between the risk score and TMB, though not statistically significant (Figure 4C). The combination of TMB and risk score significantly influences patient survival, with the combination of low TMB and low risk being associated with better survival outcomes (Figure 4D).

The previous analyses suggest that DUBGs may influence immune cell infiltration in ESCC patients. Therefore, we conducted a more detailed investigation to explore the relationship between risk scores and immune cell infiltration, immune-related gene expression, and tumor purity in ESCC patients. Using different computational methods such as CIBERSORT, MCPcounter, and xCell, we observed that the expression levels of various immune cells, including Macrophages_M0, NK_cells, Dendritic_cells, B-cells, CD4+ T-cells, Memory_B-cells, and Monocytes, were lower in the high-risk group compared to the low-risk group, particularly in the xCell analysis (Figure 5A). Additionally, correlation analysis revealed that risk scores were negatively correlated with Stromal Score, ESTIMATE Score, and Immune Score, but positively correlated with Tumor Purity, indicating that higher risk scores are associated with higher tumor purity (Figures 5B–E). Moreover, key genes such as USP37, USP22, TRAF2, SMAD3, and MTOR showed a positive correlation with immune-related genes HHLA2 and CD48 (Figure 5F), and the expression levels of HHLA2 and CD48 were significantly higher in the high-risk group (p < 0.05, Figure 5G). These analysis underscores the potential role of DUBGs in modulating immune cell infiltration and highlights the relationship between risk scores, immune gene expression, and tumor microenvironment characteristics in ESCC patients.

To gain a deeper understanding of the underlying mechanisms by which DUBGs influence the function and gene regulation in esophageal cancer cells, we conducted a single-cell level analysis. We first performed dimensionality reduction, clustering, and cell identification on the single-cell sequencing data of esophageal cancer from the GSE160269↗ dataset. A total of 13 cell types were identified: B cells, CD4Tconv cells, CD8Tex cells, DC cells, Endothelial cells, Fibroblasts, Malignant cells, Mast cells, Mono/Macro cells, Pericytes, Plasma cells, T prolif cells, and Treg cells (Figures 6A–C).

Subsequently, we analyzed the expression profiles of DUBGs within the single-cell data, observing distinct expression patterns of these genes across different cell populations (Figure 6D). The extracellular matrix is a supportive tissue structure composed of complex molecules such as collagen and fibronectin. It plays a crucial role in normal cell growth and function, while also influencing the infiltration and activity of immune cells within tumors. The extracellular matrix environment surrounding tumor cells can affect the immune response to the tumor, sometimes even forming a barrier that limits the activity of immune cells (32–35). Additionally, there were notable differences in the composition of immune and malignant cells between ESCC patients with high and low-risk scores. High-risk ESCC patients exhibited a higher proportion of malignant cells and fibroblasts, with lower infiltration of CD8 T cells and CD4 T cells compared to those with low-risk scores (Figure 7A). Moreover, the number of inferred cell interactions in the high-risk group (4633) was significantly greater than that in the low-risk group (3861) (Figure 7B, Supplementary Figure 3A). Furthermore, the interaction strength and patterns among certain cell types differed markedly between the high-risk and low-risk groups, with more complex and intense cell-type interactions observed in high-risk patients, which may be associated with their prognosis (Figure 7C).

We further identified the cell communication signaling pathways in high-risk and low-risk ESCC patient groups and generated relative information flow diagrams (Figures 7D, E) and signal pattern maps (Figure 7F). The relative information flow diagrams (Figures 7D, E) reveal that the high-risk group exhibits more information flow pathways compared to the low-risk group, with significant enrichment in pathways related to MHC-II, CD34, ICAM, and LIGHT. The signal pattern map (Figure 7F, Supplementary Figure 3B) shows that the high-risk group has increased interactions among certain cell types, particularly among endothelial cells, fibroblasts, and malignant cells. Moreover, cells in the high-risk group, such as fibroblasts, malignant cells, and macrophages, demonstrate stronger interaction intensities, which may be associated with more complex signaling within the tumor microenvironment. These findings suggest that the tumor microenvironment in high-risk patients is characterized by more active cell communication and more intricate signaling mechanisms.

Through the analysis of key DUBGs, we found that genes such as CFTR, USP2, ITCH, ESR1, AXIN1, USP37, AKT1, OTUD6B, ZC3H12A, and SMAD3 are significantly upregulated in tumor tissues compared to normal tissues (Figures 8A, C). Additionally, ESR1, OTUD6B, USP2, and USP37 showed a positive correlation with risk scores, while USP22 exhibited a significant negative correlation with risk scores (r = -0.7, p < 0.0001, Figures 8B, D), suggesting a potential role for these genes in ESCC tumor progression and patient prognosis. Subsequently, a drug sensitivity analysis revealed differences in the sensitivity to certain common drugs between high-risk and low-risk groups. Drugs such as BMS-536924, BDP-00009066, AUZ-12345, and AZ6102 showed significantly higher response levels in the high-risk group compared to the low-risk group, indicating the potential therapeutic advantages of these drugs in high-risk patients (Supplementary Figure 2A).

We incorporated a total of 14 DUBGs (USP2, ITCH, ESR1, AXIN1, MTOR, USP22, TRAF2, USP37, AKT1, AR, OTUD6B, ZC3H12A, and SMAD3) into our risk model. Among these, the MTOR gene exhibited a higher Cox coefficient in our prognostic model. Mammalian target of rapamycin (mTOR), the final component of the PI3K/mTOR pathway, plays a crucial role in regulating cell proliferation, growth, metabolism, and protein synthesis (36). Previous studies have highlighted the significant role of the MTOR gene in tumorigenesis and cancer progression. Mutations in the MTOR gene can result in the persistent hyperactivation of the mTOR signaling pathway, and over 30 mTOR gene mutations have been identified across various cancer types (37). These mutations are not only associated with enhanced mTOR protein activity but have also been linked to resistance mechanisms against mTOR inhibitors (38).

To investigate the role of MTOR in esophageal cancer, we constructed MTOR knockdown cell lines, KYSE30 and mEC25, using siRNA technology. Five days post-transfection, qRT-PCR was performed to assess MTOR expression levels, confirming the effectiveness of MTOR knockdown in the KYSE30 and mEC25 cell lines (Figures 9A, D). Functional assays revealed that MTOR knockdown significantly reduced colony formation efficiency and cell proliferation in esophageal cancer cells (Figures 9B, C, E, F). Additionally, in vivo experiments using a subcutaneous tumor model in mice, established with the mEC25 cell line, demonstrated that MTOR knockdown slowed tumor growth and prolonged survival (Figure 9G–K). These findings strongly suggest that the MTOR gene is a risk factor for esophageal squamous cell carcinoma (ESCC). Further research is warranted to develop more effective diagnostic and therapeutic strategies for ESCC patients.