Deciphering the tumor immune microenvironment: single-cell and spatial transcriptomic insights into cervical cancer fibroblasts

Abstract

Background: Cervical cancer (CC) remains a significant global health challenge despite advancements in screening, HPV vaccination, and therapeutic strategies. Tumor heterogeneity, driven by epigenetic modifications, affects immune evasion, metastasis, and treatment response. Cancer-associated fibroblasts (CAFs) play a crucial role in CC progression and therapy resistance. Single-cell sequencing offers new insights but remains underutilized in CC research. This study integrates single-cell RNA sequencing (scRNA-seq), spatial transcriptomics, and deconvolution analysis to identify key genes and immunotherapy targets. By constructing a prognostic model and exploring the immune microenvironment, we aim to provide novel insights into CC pathogenesis and potential therapeutic strategies.

Methods: We utilized scRNA-seq, spatial transcriptomics, deconvolution analysis, and pseudotime trajectory mapping to delineate fibroblast subtypes within the tumor immune microenvironment (TIME) of CC. Functional annotations, differential gene expression profiling, cell-cell communication pathways, and transcription factor networks were systematically analyzed. A prognostic model based on bulk RNA-seq data was constructed and validated through survival analysis, with correlations to immune microenvironment characteristics. Functional experiments investigated the role of SDC1, a critical mediator of fibroblast-tumor crosstalk. Additionally, Fibroblast-tumor cell co-culture systems and functional assays were employed to investigate the paracrine role of SDC1. The CAF MYH11⁺ subpopulation was isolated via fluorescence-activated cell sorting (FACS). Multiplex immunofluorescence and immunohistochemical analyses were performed on both cultured cells and human cervical cancer tissue samples to characterize the spatial distribution and dynamic remodeling of MYH11 during stromal reorganization.

Results: Six distinct fibroblast subtypes were identified, including the C0 MYH11 + fibroblasts, which exhibited unique roles in stemness maintenance, metabolic activity, and immune regulation. Spatial and functional analyses revealed that the C0 subtype is central to tumor-fibroblast interactions, particularly through the MDK-SDC1 signaling axis. The prognostic model incorporating fibroblast-specific markers demonstrated robust predictive power for patient survival outcomes. Additionally, in vitro SDC1 knockdown significantly inhibited CC cell proliferation, migration, and invasion. Fibroblasts show spatially regulated heterogeneity, with activation markers enriched in the tumor zone and MYH11 highest in normal zones, indicating dynamic stromal remodeling. C0 MYH11 + CAF Promotes Tumor Cell Proliferation, Migration, and Inhibits Apoptosis via Soluble SDC1.

Conclusion: Our results illustrate, in some ways, the possible immunomodulatory and tumor supporting roles of CAFs in CC TIME and highlight the possibility that the MDK-SDC1 pathway is a promising therapeutic target. This study not only promotes a partially new understanding of temporal heterogeneity in CC, but also provides a possible reference base for the development of new biomarkers and immunotherapy approaches to improve clinical outcomes.

Keywords: Cancer-associated fibroblasts; Cervical cancer; Immunomodulation; SDC1; Single-cell RNA sequencing; Spatial transcriptomics; Tumor immune microenvironment.

Fig. 1

The overview of this study

Fig. 2

Characterization of fibroblast subtypes in cervical cancer. A All high-quality cells were dimensionally reduced and clustered into UMAP plots, with cell types annotated based on differentially expressed genes. The differences in G2/M score, S score, nCount RNA, and nFeature RNA across cells were also displayed. B UMAP plots showed the distribution of cells from the Neoplasm group and the Normal tissue adjacent to neoplasm group. The cell cycle phase distributions within all annotated cell types were also presented, with a bar plot illustrating their respective percentages. C The top 5 marker genes of all cells and the differentially expressed genes between the two groups were visualized using bubble plots. D GO-BP analysis was performed for all annotated cells. E GO-BP analysis was conducted for the upregulated and downregulated genes in cervical cancer fibroblasts. F GSEA analysis identified enriched pathways in fibroblasts. G An enrichment network diagram illustrated pathway analysis results based on fibroblast gene sets. H Fibroblasts were classified into six subtypes with distinct characteristics based on their marker genes. I The distribution of fibroblast subtypes between the Neoplasm group and the Normal tissue adjacent to neoplasm group was shown using UMAP plots, along with the distribution of cell cycle phases within each subtype. J Stacked bar plot depicted the percentage composition of fibroblast subtypes across different sample origins. K ST maps inferred the most likely cell types at selected points within cervical cancer tissue sections. L Differences in nFeature RNA, CNV score, and nCount RNA across fibroblast subtypes were demonstrated. M ST maps displayed nFeature RNA, CNV label, and nCount RNA for fibroblast subtypes. N A heatmap showcased differentially expressed marker genes across the six fibroblast subtypes. O GO-BP analysis was performed for the differentially expressed genes of each fibroblast subtype. P Word clouds highlighted terms strongly associated with each of the six fibroblast subtypes. Q ST maps visualized C0 MYH11 + fibroblasts. R GO-BP analysis identified upregulated genes in the Neoplasm group and the Normal tissue adjacent to neoplasm group within the C0 subtype. S GSEA analysis revealed enriched pathways for C0 MYH11 + fibroblasts

Fig. 3

Scoring of Stemness, Metabolism, and apCAF Features in Fibroblast Subtypes. A The stemness of fibroblast subtypes was visualized using UMAP plots and ST feature maps. A boxplot illustrated the differences in stemness scores among subtypes. B A bubble plot displayed the expression levels of differentially expressed stemness-related genes across fibroblast subtypes. C UMAP plots highlighted the expression differences of stemness-related genes, including KDM5B, CTNNB1, CD44, and CD34, across fibroblast subtypes. D The differences in stemness-related gene expression between the Neoplasm group and the Normal tissue adjacent to neoplasm group were shown. E The significant metabolic pathways for each fibroblast subtype and for the two sample groups were identified and presented. F UMAP plots demonstrated differences in oxidative phosphorylation and glutathione metabolism across subtypes. Boxplots compared these two metabolic pathways across subtypes, cell cycle phases, and sample groups. (G) The differences in pHB, pMT, and pRP scores across fibroblast subtypes were visualized using UMAP plots and boxplots. H UMAP plots showed the differences in apCAF scores across subtypes in terms of expression and density. An ST feature map displayed apCAF scores spatially, and a boxplot quantitatively compared apCAF scores among fibroblast subtypes

Fig. 4

Pseudotime heterogeneity in fibroblast subtypes. A-B CytoTRACE analysis revealed the differentiation levels of the six fibroblast subtypes. C Genes positively and negatively correlated with CytoTRACE were identified. D UMAP plots indicated the primary pseudotime regions occupied by each subtype. E Pseudotime trajectory was constructed, progressing from left to right, showing the distribution of the six fibroblast subtypes along the trajectory. F A heatmap displayed the temporal expression patterns of differentially expressed genes across subtypes along the pseudotime trajectory. G The pseudotime trajectory was divided into five states (state1–5) based on branch points, with pie charts illustrating the percentage composition of each state. H Stacked bar plot depicted the proportional representation of each subtype within states 1–5. I GeneSwitches analysis demonstrated the distribution of surface proteins, transcription factors (TFs), and all genes along the pseudotime trajectory. J GO-BP analysis was performed for upregulated and downregulated genes at different time points along states 1–4 of the pseudotime trajectory. K Genes commonly expressed along both the High C0 and Low C0 pseudotime trajectories were identified. L Genes specifically expressed along the High C0 trajectory were highlighted, including CST3 at the early stage and COL14A1, VIM, and S100A4 at the late stage. M The expression patterns of COL14A1, VIM, and S100A4 along the pseudotime trajectory were visualized. N ST feature maps depicted the spatial expression of COL14A1, VIM, and S100A4. O A heatmap illustrated the differential expression of stemness genes along the pseudotime trajectory. P CytoTRACE2 analysis revealed the expression distribution of stemness genes with varying efficacies. Q UMAP plots showed the relative density and scoring density of CytoTRACE2 across different subtypes. R Slingshot analysis constructed two trajectories, Lineage1 and Lineage2, among subtypes. The expression of the naming gene MYH11 for C0 MYH11 + Fibroblasts was displayed along both trajectories. S UMAP plots illustrated the temporal progression of Lineage1 and Lineage2 across subtypes. (T) Based on the different states, the Slingshot trajectories (Lineage1 and Lineage2) were mapped to the distributions of each subtype. U Temporal changes in the expression of stemness-related genes (KDM5B, CTNNB1, and CD44) along Lineage1 and Lineage2 were visualized

Fig. 5

C0 MYH11 + Fibroblasts Interacted with Tumor Cells through MDK-SDC1 Crosstalk. A Circular plots illustrated the number (top) and strength (bottom) of interactions between all fibroblast subtypes and other cell types. B The bubble plots showed incoming communication patterns (top) of target cells and outgoing communication patterns (bottom) of secreting cells, respectively. C Differences in the signal intensity of pathways between the Neoplasm and Normal tissue adjacent to neoplasm groups were observed, along with a comparison of the number and strength of cell interactions in both groups. D Circular plots depicted the interaction number and intensity of C0 MYH11 + Fibroblasts with other cell types (left) and interactions from other cells to tumor cells (right). E Interaction networks among all cell types in the MDK signaling pathway were demonstrated. F The paracrine interactions of all cells in the MDK signaling pathway and the centrality scores of the MDK signaling pathway network. G Expression of various receptor proteins in different cell types along the MDK signaling pathway was shown. H A chord diagram displayed the interactions of C0 MYH11 + Fibroblasts with all other cell types. I Interactions among all cells in the MDK-SDC1 pathway were presented

Fig. 6

TFs Regulatory Network in Fibroblast Subtypes. A A heatmap displayed the differential expression of the top 5 transcription factors (TFs) across all fibroblast subtypes. B C0 MYH11 + Fibroblasts were highlighted with dark green dots on the UMAP plot (left). The ranking of regulatory factors in C0 MYH11 + Fibroblasts was shown based on the Regulatory Score (RSS) (middle), along with the expression distribution of the top 5 TFs in C0 MYH11 + Fibroblasts. C UMAP plots illustrated the expression intensity of the top 5 TFs in C0 MYH11 + Fibroblasts across various subtypes. D Boxplots visually compared the expression levels of the top 5 TFs across the different fibroblast subtypes. E ST feature maps displayed the expression of the top 5 TFs in C0 MYH11 + Fibroblasts. F A heatmap illustrated the similarity of regulatory submodules based on SCENIC recognition modules and AUCell scores. Submodules of regulators from fibroblast subtypes in both the Neoplasm and Normal tissue adjacent to neoplasm groups were identified, resulting in the identification of two submodules based on rule similarity. G UMAP plots showed the distribution of M1 and M2 modules across fibroblast subtypes. H The value of M1 and M2 modules was displayed. I The ranking of fibroblast subtype content within M1 and M2 modules was shown. J The ranking of TF expression intensity in M1 module was presented. K Differential expression levels of five TFs—CUX1, TEAD1, HOXA10, PBX1, and FOSL2—across different fibroblast subtypes were illustrated. L A bubble plot displayed the results of GO-BP analysis based on gene sets from M1 and M2 modules

Fig. 7

Key prognostic significance of C0 MYH11 + fibroblasts in cervical cancer patients. A A forest plot showed the results of univariate Cox analysis, identifying 13 C0 MYH11 + Fibroblasts-related genes associated with cervical cancer prognosis, with P ≤ 0.05. The reference line (HR = 1) distinguished protective factors (HR < 1) from risk factors (HR > 1). B Lasso regression was applied to select 10 genes contributing to the risk score, and the results were displayed in a lambda plot (lambda.min = 0.017). C A forest plot depicted the final 10 genes associated with cervical cancer prognosis. D A bar plot showed the coefficient values of the 10 risk genes. E A curve plot illustrated the risk scores in high MFRS and low MFRS groups, with scatter plots displaying survival/death events over time for the high MFRS and low MFRS groups. F A heatmap displayed the differential expression of the 10 risk genes between the high MFRS and low MFRS groups. G ROC curves showed the AUC values for 1-year (AUC = 0.80), 3-year (AUC = 0.89), and 5-year (AUC = 0.85) survival predictions. H A scatter plot displayed the distribution of genes along PC1 and PC2 for the high MFRS and low MFRS groups. I Kaplan–Meier survival curves showed survival differences between the high and low SDC1 expression groups. J A forest plot presented the results of multivariate Cox analysis for clinical factors and risk scores in the training cohort. K The nomogram model based on the MFRS was constructed, including race, age, and grade. L Scatter and heatmaps showed the correlation between the 10 prognostic genes and OS and MFRS. M Box and scatter plots illustrated the expression differences of 8 genes from the 10 risk genes between the high MFRS and low MFRS groups

Fig. 8

Tumor boundary characterization, prognostic model differences, and immune drug treatment prediction. A ST feature maps displayed the distribution of cervical cancer tumor boundaries (Bdy), malignant regions (Mal, Mal1), and normal tissue (Normal) in selected tissue sections. B Kaplan–Meier survival curves showed survival differences between the high MFRS and low MFRS groups. C Expression differences of 10 risk genes between the high MFRS and low MFRS groups were illustrated. D A pie chart displayed the proportions of different statuses, ages, races, stages, and grades in the study. E Boxplots showed the IC50 values of four immune drugs in the high MFRS and low MFRS groups

Fig. 9

Enrichment analysis of differential genes in high and low MFRS Groups. A A heatmap displayed the expression profiles of differential genes between the high MFRS and low MFRS groups. B A volcano plot showed the differentially upregulated and downregulated genes in the high MFRS and low MFRS groups. C-E Bar plots presented the GO-BP, GO-CC, and KEGG analysis of differential genes in the high MFRS and low MFRS groups. F Detailed GSEA entries based on differentially upregulated and downregulated genes were presented. G A heatmap displayed the GSVA analysis of differential gene sets in the high MFRS and low MFRS groups. H A radar chart compared the differences in entries obtained through different enrichment methods between the high MFRS and low MFRS groups. I A radar chart showed the differences in the content of various immune cell types in different states between the high MFRS and low MFRS groups

Fig. 10

Knockdown of SDC1 Inhibits Proliferation, Migration, and Invasion of Cervical Cancer Cells. A The relative expression of SDC1 mRNA in cervical cancer cells was measured after knockdown in HeLa and SiHa cells. B The relative expression of SDC1 protein in HeLa and SiHa cells following SDC1 knockdown was evaluated. C-D CCK-8 assays showed a significant decrease in cell viability in HeLa and SiHa cells after SDC1 knockdown. E EDU staining experiments indicated that downregulation of SDC1 hindered cell proliferation in both HeLa and SiHa cells. F Colony formation assays demonstrated that cells with reduced SDC1 expression formed significantly fewer colonies compared to the NC group. G-H Transwell assays revealed that SDC1 knockdown significantly reduced the migration and invasion of cervical cancer cells. I Scratch wound healing assays showed that reduced SDC1 expression significantly decreased the wound healing rate. *P < 0.01, **P < 0.001, ***P < 0.0001, and ****P < 0.00001

Fig. 11

Characterization of fibroblast heterogeneity and spatial distribution. A Flow cytometric gating strategy identifying live, single CD45⁻CD31⁻α-SMA⁺ fibroblasts and MYH11-defined subpopulations from isolated normal fibroblasts (NFs) and cancer-associated fibroblasts (CAFs). B Comparative expression of fibroblast markers showing higher MYH11 expression in NFs, whereas CAFs exhibit increased Calponin, Vimentin, and Fibronectin levels. C Spatial distribution of fibroblast activation markers across normal zone (NZ), junctional zone (JZ), and tumor zone (TZ), revealing a progressive increase in α-SMA, Calponin, Vimentin, and Fibronectin, and a reciprocal decrease in MYH11. D Representative low- and high-magnification mIHC images demonstrating zonal co-localization patterns: α-SMA, Calponin, Vimentin, and Fibronectin are predominantly expressed in TZ, partially overlapping in JZ, and minimally detected in NZ, while MYH11 is enriched in NZ with reduced expression in JZ and absent in TZ

Fig. 12

C0 MYH11 + CAF promotes tumor cell proliferation, migration, and inhibits apoptosis via soluble SDC1. A Experimental procedure for co-culture of C0 MYH11 + CAF with high- and low-invasive cervical cancer cells. Supernatants were collected and used to treat tumor cells. B shRNA-mediated silencing of SDC1 in C0 MYH11 + CAF resulted in reduced SDC1 expression. Western blot showing SDC1 levels in control and silenced groups. C Recombinant SDC1 treatment (200 ng/mL) promoted tumor cell proliferation, whereas 500 ng/mL inhibited proliferation. D Spheroid formation assay showing enhanced stemness and malignancy in tumor cells co-cultured with C0 MYH11 + CAF or treated with its supernatant. E Migration assay indicating increased migratory ability of tumor cells exposed to C0 MYH11 + CAF or its culture supernatant. F-G qRT-PCR analysis of apoptosis-related gene expression showing upregulation of anti-apoptotic and downregulation of pro-apoptotic genes in tumor cells treated with C0 MYH11 + CAF supernatant, suggesting reduced apoptosis