RGS16-driven cancer-associated fibroblasts promote esophageal squamous cell carcinoma progression via the MDK-SDC1 axis-mediated intercellular crosstalk

Abstract

AbstractSection Background

See More

Cancer-associated fibroblasts (CAFs) are critical components of the tumor microenvironment (TME) and have important roles in carcinogenesis and metastasis by facilitating cancer cell proliferation, angiogenesis, extracellular matrix (ECM) remodeling, and drug resistance. The purpose of this study is to extensively describe CAFs in esophageal cancer (ESCA) and develop a CAF-based prognostic risk score to predict clinical outcomes in patients with ESCA.

AbstractSection Methods

Data for single-cell RNA sequencing (scRNA-seq) were collected from the GEO database. The GEO and TCGA databases provided the bulk RNA-seq data and microarray data for esophageal squamous cell carcinoma (ESCC), respectively. The Seurat R program was used to identify CAF clusters from scRNA-seq data. Univariate Cox regression analysis was then used to identify prognosis-related genes linked to CAFs. Using 65 machine learning algorithms, a risk signature was created using CAF-related prognosis genes. We investigated the relationships between the CAF-related gene risk score and clinical features, mutation landscape, immunotherapy response, and medication sensitivity.

AbstractSection Results

Using scRNA-seq analysis in ESCC, 4856 genes associated with CAF clusters were identified, 14 of which were selected to construct a prognostic risk signature. Functional validation revealed that overexpression of the regulator of G-protein signaling 16-nuclear factor (RGS16) in CAFs co-cultured with the ESCC cell line KYSE520 significantly increased cancer cell proliferation, invasion, and migration. The secreted midkine (MDK) coupled to its receptor syndecan1 (SDC1) on ESCC cells, further helping their malignant tendencies.

AbstractSection Conclusions

Comprehensive assessment of CAF heterogeneity in ESCC sheds light into the mechanisms of immunotherapy resistance and suggests prospective options for creating novel treatment therapies. Furthermore, we demonstrate that RGS16+ CAFs promote ESCC progression through the nuclear factor kappa B(NF-κB)-MDK-SDC1 axis, highlighting their crucial involvement in tumor-stromal interaction. These findings emphasize the therapeutic potential of targeting RGS16+ CAFs, which presents a promising method for disrupting the tumor-promoting microenvironment and improving clinical outcomes in ESCC patients.

AbstractSection Clinical trial number

Not applicable.

Introduction

Esophageal cancer is the seventh leading cause of cancer-related mortality worldwide, accounting for 445,129 deaths in 2022 which was 4.6% of all cancer-related deaths [1]. Two main histological subtypes are identified: esophageal adenocarcinoma (EAC) and ESCC, with ESCC accounting for nearly 85% of all newly diagnosed cases in 2020 [2]. ESCC prognosis remains poor due to high rates of recurrence and metastasis despite intensive surgical approaches, due to tumor cell proliferation, migration, invasion, and epithelial-mesenchymal transition (EMT) [3,4,5]. However, the processes that control ESCC progression are not completely understood.

The TME is a multicellular environment characterized by intricate interactions between the stroma and the tumor which is extensively recognized [6]. The TME is critical for driving important oncogenic processes like proliferation, angiogenesis, apoptosis inhibition, immune suppression, and immune evasion. The TME undergoes ongoing cellular and molecular remodeling due to reciprocal communication between neoplastic cells and stromal or immune components, boosting the cancer cells’ metastatic potential [7, 8]. CAFs are a critical stromal element of the TME and are found in a variety of solid tumors, making them promising therapeutic targets [9]. Activated CAFs stimulate tumor development, angiogenesis, invasion, metastasis, ECM remodeling, and chemotherapy resistance through various methods. CAFs communicate with immune cells and other components of the tumor immune microenvironment (TIME) by releasing cytokines, growth factors, chemokines, exosomes, and other effector molecules. This causes the creation of an immunosuppressive TIME, which allows cancer cells to elude immune surveillance [10]. Consequently, CAFs are common targets in anti-tumor immunotherapy [11]. However, the methods by which CAFs influence anti-tumor immune responses in solid tumors remain poorly understood.

Recent advances in scRNA-seq technology have transformed biological system exploration [12]. This method provides large-scale, parallel characterization of various cell populations at the transcriptomic level, bypassing the limits of conventional transcriptomics and allowing for a more in-depth examination of TME characteristics [10, 13, 14]. Due to this benefit, various studies have focused on identifying novel biomarkers for malignant tumors by integrating scRNA-seq with bulk RNA sequencing (RNA-seq).

The aim of this study is to integrate scRNA-seq and RNA-seq data by combining datasets of patients with esophageal cancer from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO). Using 65 machine learning algorithms [15], we plan to link CAF marker genes to clinical feature modules to identify CAF-related genes. Based on these genes, we plan to create a machine learning-based prognostic risk score model for esophageal cancer. We performed a comprehensive bioinformatics analysis to evaluate mutational profiles, clinical characteristics, immunotherapy responses, and drug efficacy in high- and low-risk esophageal cancer patients. Furthermore, RGS16 was identified as a key CAF-related gene and its role in promoting ESCC proliferation, invasion, and migration was established by activating CAFs and the MDK-SDC1 intercellular communication pathway through NF-κB signaling cascade. These findings provide new insights into the treatment and improved patient prognosis of esophageal cancer.

Materials and methods

Data acquisition

The scRNA-seq data for esophageal cancer was obtained from the GSE196756 database for this study (https://www.ncbi.nlm.nih.gov/geo/), including six samples (three primary tumors and their paired samples). Transcriptome data, single nucleotide variation, and copy number variation data, and corresponding clinical information for ESCA were obtained from the TCGA database (https://portal.gdc.cancer.gov/). The cases not having either survival data nor outcome status were removed from the transcriptome data, and finally, there were 160 tumor samples and 4 paraneoplastic samples. The TCGA data served as the training cohort, while the GSE53625 dataset of ESCA samples from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) was used as the validation cohort, which included 179 tumor samples and 179 normal samples.

Identification of cancer-associated fibroblast marker genes

The Seurat R program was used to analyze and process the transcriptomic data to generate a Seurat object [16]. Cell populations were annotated using a reference gene list from 713 samples in the “Human Primary Cell Atlas Data.” Cell populations were clustered but unsupervised using the UMAP function, and then visualization analysis was performed [17]. Marker genes for distinct cell types were identified using the “FindAll-Markers” function with stringent thresholds (|log2FC| > 1 and FDR < 0.05). The identified marker genes for CAFs were used for further analysis.

Functional enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed using the clusterProfiler R program to identify the functions of the selected core genes, followed by correction for multiple testing with the Benjamini–Hochberg procedure.

Construction of risk score feature models and prognostic analysis

To identify genes linked with prognosis, we used Univariate Cox regression analysis in the survival package [18]. To narrow down on the number of genes identified, the TCGA dataset served as the training dataset, and GSE53625 was used as the independent validation dataset. Model screening was conducted using 65 machine learning methods, followed by multifactorial Cox analysis to determine risk coefficients for each gene. The concordance index (C-index) was generated for the various approaches, and patients were divided into high- and low-risk groups based on the median risk score. Kaplan-Meier (KM) survival analysis was performed to investigate prognostic differences between the high- and low-risk groups, followed by a similar analysis in the validation cohort [18]. The model’s ability to predict 1-, 3-, and 5-year survival was measured using receiver operating characteristic (ROC) analysis.

Evaluation of the model

To evaluate the ability of the model to distinguish survival outcomes, we used the timeROC package to calculate the AUC values for 1-,3-, and 5-year OS in patients. Subsequently, we used the survival package to compute the C-index values across different datasets to assess the accuracy of the model in predicting patient prognosis risks. Additionally, to compare the accuracy of the current modle with other clinical features in predicting survival outcomes, we used the compareC package to calculate the C-index values for TCGA and GSE53625, compared whether there were significant differences in these C-index values compared to our current model. In this section, these R packages were loaded by using the online tool http://sxdyc.com/diagnosticmodel. This tool provides a comprehensive platform integrating numerous machine learning algorithms, which can be used for feature selection and model construction.

Single nucleotide variation and clinical characterization of risk score models

We used the maftools software package (v2.4.05) to show mutation data across risk subgroups and investigated the relationship between risk scores and tumor grading [19]. Genes with mutation frequencies > 3 in ESCC were identified and visualized using the oncoPrint function in the ComplexHeatmap package [20]. The chi-squared test was used to determine the distribution of mutations across subtypes while correlation matrices were analyzed using Spearman’s correlation coefficient [21]. Differences between the groups were compared using the Wilcoxon test [22]. Survival differences were determined using KM curves and log-rank tests. Univariate and multivariate Cox regression models were used to determine the independent prognostic significance of risk ratings in combination with clinical data.

Immunotherapy response analysis

The Tumor Immune Dysfunction and Exclusion (TIDE) score was used to predict outcomes and immunotherapy responses in cancer patients [23]. Next, the TCGA-ESCC cohort expression profiles were uploaded to the TIDE website (http://tide.dfci.harvard.edu) to obtain TIDE scores for each sample [24]. The immunotherapy dataset GSE91061 was employed to explore immunotherapy responses in both high- and low-risk groups.

Drug sensitivity prediction

The R package “oncoPredict” was used to evaluate drug sensitivity of common chemotherapeutic and targeted therapies across risk groups [25]. Therapeutic efficacy was quantified using the Genomics of Drug Sensitivity in Cancer (GDSC) dataset, by the half-maximal inhibitory concentration (IC50); lower IC50 values indicated higher drug efficacy.

Construction of a cell-to-cell interaction network

To investigate cell-cell interactions between CAFs and epithelial cells, we used the R package ‘CellChat’ (v1.1.2) [26]. T Ligand-receptor pairs in CellChat were obtained from previous studies [27].

Cell culture and transfection

The Chinese Academy of Sciences provided the ESCC cell line KYSE520 from their Cell Bank of the Type Culture Collection. Esophageal CAFs HUM-iCELL-d042 were bought from iCell (iCell Bioscience Inc, Shanghai). Cancer cells were cultured in Roswell Park Memorial Institute 1640 medium (Gibco, Sigma-Aldrich, USA) supplemented with 10% FBS (Gibco, Sigma-Aldrich) at 37 °C with 5% CO2. CAFs were cultured in the iCell Primary Fibroblast Culture System (Cat: PriMed-iCELL-003) and seeded in six-well plates. Before transfection, the medium was replaced with Opti-MEM. Small interfering RNA (siRNA; GenePharma, China) and Lipofectamine 3000 transfection reagent (Life Technologies, USA) were used to transfect cells following the manufacturer’s instructions. The target siRNA sequences were as follows: RGS16-siRNA: 5'-CCGAGAACTGACCAAGACAAA-3'; SDC1-siRNA: 5'-CCGCAAATTGTGGCTGTAAAT-3'.

Co-culture system

To create a CAF-ESCC co-culture system, a Transwell non-contact system (0.4 µm pore size) was used to allow paracrine signaling between cells without direct physical contact. ESCC cells 8 × 10⁴ were plated in the upper chamber while 8 × 10⁴ CAFs were seeded in the lower chamber. The co-culture was maintained under standard conditions (37 °C, 5% CO₂) for 3 days. Following the co-culture period, ESCC cells were removed from the upper chamber for further analysis of malignant properties.

Colony formation assay

This was conducted to assess the proliferative potential of cells. A total of 500 cells were seeded per well of a 12-well plate (Corning, USA) and cultured at 37 °C with 5% CO₂ for 2 weeks. Following incubation, the supernatant was collected, and cells were washed twice with phosphate buffered saline (PBS). Cells were fixed with 4% paraformaldehyde for 30 min at room temperature and stained with 0.1% crystal violet (Beyotime, Shanghai, China) for 5 min. Colonies in each well were visually examined and manually counted.

Transwell assay

For the cell migration assays,24-well Transwell plates with 8 μm polycarbonate membranes (Corning, USA) were utilized. The cells were suspended in 200 μL serum-free medium and 200 μL was added into the upper chamber of the Transwell plate. The lower chamber was filled with 600 μL of 20% FBS which acted as a chemoattractant. Following the incubation for 24–48 h, non-migrated cells on the upper membrane surface were removed using a PBS-wetted cotton swab. The lower surface of the migrated cells was fixed with 4% paraformaldehyde and stained with 0.1% crystal violet for 20 min. Excess stain was rinsed with slow-running water, and the membrane was air-dried. The migrated cells were visualized by a microscope by counting 3–5 random fields per membrane. In case of invasion assays, the rest of the steps remained the same as those for the migration assay, except the upper Transwell membrane was pre-coated with diluted Matrigel matrix gel (ABW, China).

Western blotting

Total protein from exosomes was extracted using RIPA buffer (Beyotime, Shanghai, China) and quantified using a BCA Protein Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA). Equal amounts of protein (30 μg/sample) were loaded and separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis and transferred onto a polyvinylidene difluoride membrane. The membrane was blocked with 5% nonfat milk, rinsed, and incubated with primary antibodies overnight at 4 °C.The primary antibodies included RGS16 (Abcam, ab119424, 1:1000), NF-κB (Abcam, ab16502, 1:1000), phospho-NF-κB (Abcam, ab86299, 1:1000), MDK (Abcam, ab52637, 1:1000), and SDC1 (Abcam, ab128936, 1:1000). Blots were probed with HRP-conjugated anti-rabbit or anti-mouse secondary antibodies and developed using an enhanced chemiluminescence reagent (Millipore, USA).

Statistical analyses

Bioinformatics analyses were performed using R software (v4.2.1; https://www.r-project.org/). Cellular experimental data were analyzed using GraphPad Prism 8.0 (San Diego, CA, USA). The analyses results are presented as mean ± standard deviation (SD). Differences between groups were analyzed using two-tailed Student’s t-test or one-way ANOVA, as appropriate. Statistical significance was defined at *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001.

Results

Identification of CAFs marker genes

We performed Uniform Manifold Approximation and Projection (UMAP) clustering analysis on tumor microenvironment cells and separated them into 11 clusters based on marker gene expression. These clusters comprised epithelial cells, fibroblasts, macrophages, NK cells, T cells, B cells, plasma cells, endothelial cells, mast cells, smooth muscle cells, and neutrophils (Fig. 1A). Marker gene expression in each cluster was presented through a bubble plot (Fig. 1B). Cells were reclassified into 11 clusters based on marker gene expression patterns, and the proportion of each cell type across samples was calculated (Fig. 1C). Markers of significant heterogeneity in cell composition across samples, revealed considerable variability in the TME. Notable differences were observed when we compared cell type proportions between tumor and paired normal tissues (Fig. 1F). Differential expression analysis for each cell subtype was performed (Fig. 1D), followed by functional enrichment analysis of the differentially expressed genes (DEGs). Fibroblast-related genes were enriched predominantly in key pathways, including tyrosine metabolism, retinol metabolism, drug metabolism-cytochrome P450, transforming growth factor-beta signaling, and cytoskeleton regulation in muscle cells (Fig. 1D). Analysis of cellular composition proportions across different samples revealed a significant increase in fibroblast proportion (Fig. 1E). These findings suggest the role of fibroblasts in modulating these pathways in the TME, potentially influencing tumor progression and therapeutic responses.

Fig. 1

Identification of CAFs marker genes. A T-distributed stochastic neighbor embedding plot illustrating cell annotation following single-cell dimensionality reduction and clustering. B Expression patterns of marker genes across different cell populations. C Differentially expressed marker genes among distinct cell types. D KEGG functional enrichment analysis of differentially expressed marker genes in various cell types. E Representation of the proportion of cell populations across different samples. F A comparison of cell population proportions in tumor and paired normal tissue samples

Construction and validation of a prognostic signature based on CAFs marker genes

Overall, 4,856 CAF-associated DEGs were identified. Among these, univariate Cox regression analysis revealed that 24 genes were significantly associated with ESCC prognosis in the TCGA dataset (p < 0.05). These genes included HIGD1A, CST3, KCNQ1OT1, RGS16, MT1E, MT1M, EGFR, PCSK5, CCL11, MEG3, PHYHD1, VASN, DCLK1, SH3PXD2B, CDC42EP1, LMAN1, S1PR2, FAM3C, ARHGAP21, PDLIM2, ATL3, TLE1, EEF2K, and TRIOBP (Fig. 2A). Sixty-five machine learning algorithms were used to identify the optimal model by adopting the TCGA dataset as the training set and GSE53625 as the validation set (Fig. 2B). Subsequently, their risk coefficients were determined through multivariate Cox regression analysis on the selected genes (Fig. 2C). The C-index was calculated for each method, revealing that the Elastic Net (Enet) method with alpha = 0.8 yielded the best average C-index; subsequently, 14 genes could be identified for model construction. The risk score was calculated as follows:

Fig. 2

Construction and validation of a prognostic signature based on CAFs marker genes. A Univariate analysis of CAFs marker genes in TCGA dataset. B Calculation of the C-index across different datasets using 65 machine learning algorithms. C Multivariate COX analysis identified 14 genes using the enet (alpha = 0.8) method. D-E KM and ROC curve analyses for the TCGA dataset. F-G KM and ROC curve analyses for the GSE53625 dataset

RiskScore = −0.4815×ATL3–0.0589×CCL11–0.2148×CDC42EP1–0.027×EEF2K–0.3163×FAM3C + 0.2622×HIGD1A +0.3078×KCNQ1OT1–0.188×LMAN1–0.2198×MEG3–0.1451×PDLIM2–0.0716×PHYHD1 + 0.3029×RGS16–0.4179×TLE1–0.708×TRIOBP. ESCC patients were stratified into high- and low-risk groups based on the median risk score. Kaplan-Meier survival analysis in the TCGA-ESCA cohort revealed significantly worse overall survival in the high-risk group relative to the low-risk group (p < 0.0001; Fig. 2D). The area under the curve (AUC) values for 1-, 3-, and 5-year survival predictions in the TCGA cohort were 0.77, 0.74, and 0.69, respectively (Fig. 2E). the risk score model was applied to the GSE53625 dataset to validate these findings. Consistent with TCGA results, compared to the low-risk group, the high-risk group exhibited significantly worse prognosis (p < 0.005; Fig. 2F). The AUC values for 1-, 3-, and 5-year survival predictions in the GSE53625 cohort were 0.67, 0.65, and 0.67, respectively (Fig. 2G). In summary, the high-risk group had consistently worse long-term survival outcomes in both cohorts, highlighting the robustness and prognostic significance of the CAFs-based gene signature.

Mutation analysis and clinical characterization of the RiskScore model

Somatic mutation profiles in high- and low-risk groups demonstrated the top 10 most frequently mutated genes: MUC17, ROBO2, APC, CHD3, ARAP2, CNTLN, FAM135B, MYH14, SNTG1, and SORCS1 (Fig. 3A). Both univariate and multivariate Cox regression analyses were performed to evaluate the prognostic value of the risk score on TCGA clinical data, including hazard ratios (HR), 95% confidence intervals (CI), and P-values. Clinical variables included age, gender, T stage, N stage, overall stage, and risk score grouping. The risk score was found to be a statistically significant prognostic factor in both univariate and multivariate analyses (Figs. 3B−D). This was further used to validate the clinical applicability of the RiskScore model. We conducted univariate and multivariate Cox regression analyses on GSE53625 clinical dataset. Similar to the TCGA results, the risk score remained a significant prognostic indicator in both analyses (Figs. 3E−G). These findings highlight the independent reliability of the RiskScore model across datasets, indicating its robustness and clinical utility.

Fig. 3

Mutation analysis and clinical characterization of the RiskScore model. A Distribution of mutation frequencies of differentially mutated genes in the risk subgroups. B Differences in distribution of existing subtypes and clinical features in the risk subgroups within TCGA dataset. C Univariate COX analysis of clinical features and RiskScore in the TCGA dataset. D Multivariate COX analysis of clinical features and RiskScore in the TCGA dataset. E Differences in distribution of existing subtypes and clinical features in the risk subgroups within the GSE53625 dataset. F Univariate COX analysis of clinical features and RiskScore in the GSE53625 dataset. G Multivariate COX analysis of clinical features and RiskScore in the GSE53625 dataset

Prediction of immunotherapy response and drug sensitivity based on RiskScore

We utilized the TIDE online tool to evaluate the potential efficacy of immunotherapy, and to predict the immunotherapy response in the TCGA dataset. Significantly lower TIDE scores were observed in the high-risk group compared to the low-risk group, suggesting that patients under high-risk may benefit more from immunotherapy (Fig. 4A). Furthermore, analysis of the correlation between risk score and TIDE scores revealed a significant negative correlation (Fig. 4B). These findings were further validated by comparing expression profiles of high- and low-risk groups with the GSE91061 dataset, which includes PD-1 treatment response data. The analysis revealed that the correlation was significant (p < 0.05) between the TCGA high-risk group and the PD-1 treatment response group, suggesting a favorable response to PD-L1 therapy in patients in the high-risk group (Figs. 4C and D). We further assessed drug sensitivity by calculating IC50 and analyzing its correlation with the risk score. Applying a threshold of |R| > 0.25 and p < 0.05, 18 drugs with significant correlations were identified, including 17 exhibiting positive correlations and 1 exhibiting a negative correlation. Notably, drugs such as Linifanib, Midostaurin, Crizotinib, and Sunitinib showed significant positive correlations with the RiskScore, while Amuvatinib exhibited a significantly negative correlation (Fig. 4E). Further analysis of the differences in IC50 values between the high-risk and low-risk groups in the TCGA dataset revealed significant differences among several drugs, including Amuvatinib, Bortezomib, Cabozantinib, Capivasertib, Crizotinib, Foretinib, JAK1, JNJ38877605, Linifanib, Midostaurin, PI3Ka, Pictilisib, Ponatinib, Sunitinib, and TGX221 (Fig. 4F).

Fig. 4

Prediction of immunotherapy response and drug sensitivity based on risk Score. A Differences in TIDE scores between high- and low-RiskScore groups in the TCGA dataset. B Correlation analysis between RiskScore and TIDE results in the dataset. C Proportion of immune therapy response groups in the RiskScore subgroups within the TCGA dataset. D Submap immune mapping between the TCGA dataset and GSE91061 dataset. E Correlation analysis between RiskScore and drug IC50 values. F Differences in drug IC50 values of the RiskScore subgroups. Statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

Identification of the key CAFs gene RGS16 and functional analysis of RGS16+CAFs vs. RGS16-CAFs

Analysis of the expression of 14 prognostic genes in normal and tumor samples from the TCGA and GSE53625 datasets revealed significant differential expression of only three genes (PDLIM2, PHYHD1, and RGS16) in the TCGA dataset, whereas 10 genes (ATL3, CDC42EP1, FAM3C, HIGD1A, KCNQ1OT1, MEG3, PHYHD1, RGS16, TLE1, and TRIOBP) exhibited significant differences in the GSE53625 dataset (Figs 5A, B). By overlapping results from both datasets, we identified RGS16 and PHYHD1 as key genes (Fig. 5C). We extracted CAF data from the single-cell dataset GSE196756 to further assess the role of RGS16 in CAFs. CAFs were divided into two groups based on RGS16 expression levels: RGS16-negative (count = 0) and RGS16-positive (count > 0). Using the FindMarkers function, 199 DEGs between the two groups were identified, including 90 upregulated and 109 downregulated genes in the RGS16+CAFs group (Fig. 5D). Functional enrichment analysis using clusterProfiler showed that these DEGs were enriched in pathways involving lipid and atherosclerosis, TNF signaling, and NF-κB signaling (Fig. 5E).

Fig. 5

Identification of the key CAFs gene RGS16 and functional analysis of RGS16+CAFs vs. RGS16-CAFs. A Expression profiles of CAF-related genes in TCGA dataset. B Expression profiles of CAF-related genes in the GSE53625 dataset. C Overlapping CAF-related genes between the TCGA and GSE53625 datasets. D Volcano plot of DEGs between RGS16+ CAFs and RGS16- CAFs groups. E KEGG pathway enrichment analysis of DEGs between RGS16+ CAFs and RGS16- CAFs groups. Statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

The role of RGS16 in CAFs on ESCC progression

The effects of RGS16-transfected fibroblasts on ESCC cell proliferation, invasion, and migration were examined in a co-culture system to further evaluate the role of RGS16 in promoting ESCC malignancy. Co-culturing ESCC cells (KYSE520) with si-RGS16 transfected CAFs significantly reduced their proliferative, invasive, and migratory properties (Figs. 6A and B); conversely, co-culturing ESCC cells with OE-RGS16 CAFs markedly increased these properties (Figs. 6C and D). These findings demonstrate that RGS16+CAFs critically promotes ESCC malignancy.

Fig. 6

The role of RGS16 in CAFs on ESCC progression. A-B WB results showed RGS16 expression in CAFs (d042 cell line) transfected with negative control (NC) and si-RGS16. C-D Western blotting results showing RGS16 expression in CAFs (d042 cell line) transfected with NC and OE-RGS16. E-H colony formation assays and transwell experiments showed the proliferation, invasion, and migration of ESCC cells co-cultured with CAFs transfected with NC, si-RGS16 or OE-RGS16. Scale bars: 100 μm. Statistical significance: *p < 0.05, **p < 0.01

Interaction of ESCC cells and RGS16+ CAFs via the MDK–SDC1 pathway

We performed intercellular communication analysis using ligand-receptor pairs to further explore interactions between malignant epithelial cells and RGS16+ CAFs. The results showed increasing interactions along the MDK-SDC1, MDK-SDC2 (syndecan2), and MDK-SDC4 (syndecan4) axes, with the MDK-SDC1 interaction being the most prominent (Fig. 7A). We transfected ESCC cells with si-SDC1 to validate the role of MDK binding to SDC1 in ESCC progression (Figs. 7B and C). Adding recombinant MDK protein to ESCC cells with knocked-down SDC1 (si-SDC1) significantly reduced their proliferative, invasive, and migratory properties (Fig. 7D). In contrast, adding recombinant MDK protein to SDC1 overexpressing ESCC cells (OE-SDC1) markedly enhanced proliferation, invasion, and migration (Fig. 7E). These findings demonstrate that the MDK-SDC1 interaction plays a crucial role in promoting the malignant characteristics of ESCC cells.

Fig. 7

Interaction of ESCC cells and RGS16+CAFs via the MDK–SDC1 pathway. A ligand-receptor pair analysis of RGS16- CAFs and RGS16+ CAFs as the source across different cell types. B-C WB results showed SDC1 expression in ESCC cells (KYSE520 cell line) transfected with negative control (NC) and si-SDC1. D-E WB results showed SDC1 expression ESCC cells (KYSE520 cell line) transfected with NC and OE-SDC1. F-I colony formation assays and transwell experiments were performed to evaluate the effects of MDK on the proliferative, invasive, and migratory properties of ESCC cells transfected with NC, si-SDC1, or OE-SDC1. Scale bars: 100 μm. Statistical significance: *p < 0.05, **p < 0.01

The RGS16-NF-κB-MDK signaling pathway mediates the progression of ESCC cells in co-culture with CAFs

To investigate how RGS16+CAFs regulates ESCC malignancy, we examined MDK expression in RGS16-transfected CAFs through western blot analysis. The results demonstrated that si-RGS16 in CAFs reduced phosphorylated NF-κB (pNF-κB) levels, subsequently reducing MDK expression (Fig. 8A). In contrast, OE-RGS16 increased pNF-κB levels, resulting in elevated MDK expression (Fig. 8B). To investigate the role of MDK in ESCC malignancy, we treated OE-RGS16 CAFs with an MDK inhibitor (iMDK) and co-cultured them with ESCC cells. The results demonstrated that iMDK treatment significantly attenuated the enhanced proliferation, invasion, and migration of ESCC cells induced by co-culture with OE-RGS16 CAFs (Fig. 8C-F). These findings suggest that RGS16 activates the NF-κB pathway, upregulates MDK expression, and promotes ESCC cell proliferation, invasion, and migration.

Fig. 8

RGS16-NF-κB-MDK signaling pathway mediates the progression of ESCC cells in co-culture with CAFs. A-B NFκB, p-NFκB, and MDK expression in CAFs transfected with NC and si-RGS16 shown by western blotting (WB). C-D NFκB, p-NFκB and MDK expression in CAFs transfected with NC and OE-RGS16 shown by WB. E-F colony formation assays and transwell experiments to assess the effects of iMDK on the proliferation, invasion, and migration of ESCC cells co-cultured with CAFs transfected with OE-RGS16.Scale bars: 100 μm. Statistical significance: *p < 0.05, **p < 0.01

Discussion

ESCC has recently garnered attention due to its high malignancy, recurrence rates, poor prognosis, and unremarkable early symptoms [28]. ESCC is often associated with visceral metastasis and results in a high mortality rate. Conventional treatments have shown limited improvement in ESCC prognosis, spurring the search for novel treatment options [29]. While advances in genetics and molecular biology have enhanced our understanding of ESCC pathogenesis, further exploration of ESCC progression mechanisms is crucial for developing effective treatments. The TME consists of non-cancerous cells, tumor components, and the molecules they produce and release [30]. The interplay between tumor cells and the TME critically determines tumor initiation, progression, metastasis, and the response to therapy [31]. Several studies have provided compelling evidence supporting the dynamic crosstalk between tumor cells and stromal cells, which plays a pivotal role in tumor progression [32]. By understanding the mechanisms of this interaction, there is an potential to develop enhanced therapies that simultaneously target multiple components of the TME, ultimately increasing the likelihood of favorable patient outcomes [33]. Considering the CAFs are known to play important roles in tumor initiation and progression [34], we performed a comprehensive characterization and classification of CAFs in ESCC using scRNA-seq data. It is worth noting that there is growing evidence confirming that CAF-related gene signatures hold significant prognostic value in ESCC [35, 36]. A novel CAF-based risk signature was developed using the genes that were found within the CAF clusters which correlated to ESCC prognosis.

We, herein, propose and validate a prognostic model that enabled us to classify ESCC patients into low- and high-risk cohorts using 14 CAF related genes. The high-risk cohort was shown to have worse prognosis, higher tumor mutational burden, and higher level of immunotherapy sensitivity, which demonstrated that CAF related genes are important indicators of ESCC progression and treatment outcome. Currently, chemotherapy remains a cornerstone treatment for advanced and metastatic ESCC [37]. However, the TME of each patient exhibits distinct characteristics due to tumor heterogeneity, leading to variations in chemotherapy sensitivity among ESCC patients [38].

Analysis of drug sensitivity showed that specific drugs, such as Linifanib and Crizotinib, which had the strongest correlation to the risk score may be better targeted therapies for high-risk patients. Among the 14 CAF-related genes, RGS16 emerged as a core gene mediating CAF-cancer cell interactions. Like other members of the RGS family, RGS16 is implicated in the regulation of signaling networks within the cell. Recent reports suggest that RGS proteins are misregulated in a number of prevalent diseases such as cancer, cardiovascular disorders, and neurodegenerative conditions [39, 40]. RGS16 protein is known to influence key processes such as immunity, cellular inflammation, and tumor proliferation and migration [41]. RGS16 was recently found to be associated with the regulation of G-protein-coupled receptor signaling and other pathways in tumors, and it is closely associated with multiple cancers, including breast cancer, pancreatic cancer, and colorectal cancer [42,43,44]. RGS16 function in gliomas, showing that its expression is highly elevated in glioma tissues and that in turn, RGS16 promotes tumor cell proliferation, migration, and the EMT process [45]. Wang et al. [46] further confirmed markedly elevated expression of RGS16 in glioma tissues and is strongly correlated with poor survival rates. Loss of RGS16 expression was shown to significantly reduce the proliferation and migration of glioma cells. However, the role of RGS16 in ESCC and CAFs has not been previously explored. Our functional analyses revealed that knockdown of RGS16 in CAFs significantly inhibited the proliferation, invasion, and migration of ESCC cells, and an opposite effect was noted when RGS16 was overexpressed. These findings suggest a pro-tumorigenic role of RGS16+ CAFs in ESCC by enhancing the malignant behavior of cancer cells.

We found that RGS16 activates the NF-κB signaling pathway in CAFs, leading to increased secretion of MDK, a heparin-binding growth factor, now recognized as a multifaceted factor that contributes to both paraneoplastic tissue homeostasis as well as disease progression. It is prominent in various human malignancies and serves as a mediator for the acquisition of critical cancer hallmarks, including cell growth, survival, metastasis, migration, and angiogenesis [47]. In this study, we used an MDK inhibitor (iMDK) and noted that proliferation, invasion, and migration capabilities of ESCC cells was greatly diminished. This observation indicates that RGS16+ CAFs promote the malignant phenotype of ESCC cells through MDK secretion. Ligand-receptor pair analysis of intercellular interactions revealed that RGS16+ CAFs exert their effects via the MDK-SDC-1 pathway in malignant epithelial cells. SDC1, one of the four members of the syndecan family, is a cell surface heparan sulfate proteoglycan consisting of 288 amino acids and ranking as the second largest in the syndecan family [48]. SDC1 has a role in modulating cell-cell and cell-matrix interactions [49] and is also involved in regulating cell proliferation and invasive growth [50]. The expression of SDC1 is highly specific for multiple myeloma [51], and its overexpression has also been reported in various cancers, including breast, pancreatic, gallbladder, endometrial, ovarian, prostate, and bladder cancers. In this study, we demonstrated that even when MDK was added to SDC1-knocked-down ESCC cells, the proliferation, invasion, and migration abilities of ESCC cells were still reduced. This suggests that MDK promotes the malignant phenotype of ESCC cells by binding to the membrane surface receptor SDC1, highlighting that targeting the MDK-SDC1 axis in ESCC holds therapeutic potential (Fig. 9). Despite the insights gleaned regarding the function of RGS16+ CAFs in ESCCs progression, some limitations still remain. To start, the functional studies were largely performed in vitro, and the role of the NF-κB-MDK-SDC1 axis in animal models is necessary to confirm its role in vivo. Additionally, the way of RGS16’s control of NF-κB activation in the CAFs is unknown. Further work should be directed towards understanding the upstream modulators of RGS16 as well as other interactions within the TME signaling network.

Fig. 9

Graphical illustration of the mechanism. RGS16+CAFs activate the NF-κB signaling pathway, resulting in the secretion of MDK. MDK then binds to the membrane surface receptor SDC1 on ESCC cells, thereby promoting the malignant phenotype of ESCC. (created in BioRender. Luo, D. (2025) https://BioRender.com/gdhxa1u)

Conclusions

In a single sentence, we developed a novel risk score model based on integrative scRNA-seq and transcriptomic data, that sifted through a CAF-related gene profile, which serves as a fruitful resource for prognosis and therapeutic targeting in ESCC. RGS16 was found as a prominent CAF-associated gene and we demonstrated RGS16+ CAFs mediating ESCC advancement through the NF-κB-MDK-SDC1 signaling axis. Manipulation of the RGS16-MDK-SDC1 axis is an approachable therapeutic option for patients with ESCC, especially those with high risk, while offering new avenues for harnessing treatment for esophageal cancer.