Effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer

Ovarian cancer (OC) ranks first in mortality among gynecologic malignancies¹. High-grade serous ovarian cancer (HGSOC) is the most common and deadly subtype of OC². Over 75% of HGSOC patients are diagnosed at an advanced stage with extensive malignant ascites and omental metastasis³. Complete staging surgery followed by platinum-based chemotherapy is the standard treatment for HGSOC⁴. However, over 25% of patients develop chemotherapy resistance within 6 months of initial treatment, and 70% experience recurrence within 2–3 years, ultimately succumbing to acquired drug resistance⁵. Therefore, platinum drug resistance is the primary cause of poor prognosis, recurrence, and death in HGSOC patients. Identifying clinical indicators that are closely associated with platinum chemotherapy resistance and HGSOC treatment prognosis is of great importance.

Inter- and intratumoral heterogeneity is a crucial factor influencing patient prognosis and treatment outcomes⁶. Intertumoral heterogeneity between patients is associated with diverse genetic mutations, epigenetic modifications, and transcriptional alterations. Intratumoral heterogeneity refers to the presence of different cell types within a tumor, resulting from interactions between tumor cells and the tumor microenvironment (TME)⁶.

Previous research has focused on intertumoral heterogeneity by analyzing genomic characteristics to obtain molecular subtypes of tumors, identify different patient subgroups, and select targeted therapies. For instance, the composition types of cells within the tumor microenvironment can classify tumors as either cold tumors or hot tumors³. In HGSOC, it could be classified into four molecular subtypes: mesenchymal, immunoreactive, differentiated, and proliferative⁷. Among these, the immunoreactive subtype demonstrates a better prognosis and treatment response⁷. Identifying different mutational characteristics is crucial for treatment selection in patients. BRCA1/2, TP53, and PI3K/AKT/mTOR pathway genes, which are frequently altered in ovarian cancer, serve as prime targets for clinical diagnosis and therapy⁸. However, intratumoral heterogeneity causes drug resistance, leading to therapeutic failure⁶. Thus, we need to integrate molecular features of both inter- and intratumoral heterogeneity with functional heterogeneity to improve patient subclassification and the response to therapy.

Tumor heterogeneity is thought to adhere to the fundamental principles of Darwinian evolution, where individual cells harboring heritable mutations that enhance adaptability gain a survival advantage⁹. Natural selection drives clonal expansion, leading to the emergence of subclones with varying proliferative, migratory, and invasive capabilities¹⁰. The evolution of adaptive clones occurs within the dynamic tissue environment known as the TME, forming a complex local tumor ecosystem. Changes within the TME further influence genetic diversification and phenotypic outcomes, driving the compositional tumor heterogeneity^11,12,13. Therefore, identifying the distinct tumor cell clone subtypes within HGSOC and the microenvironmental cell types they interact with, exploring the interactions between the microenvironment and tumor cells, and elucidating the patterns of tumor evolution are crucial for the treatment and prognosis prediction of HGSOC patients.

The TME, which includes stromal cells, immune cells, extracellular matrix, and a variety of soluble factors, can influence tumor progression and the response to therapy¹³. Single-cell and spatial transcriptomics technologies have significantly advanced the investigation of heterogeneity in HGSOC^7,14. For example, the presence of tumor-associated macrophages (TAMs) has been associated with poor prognosis in HGSOC, as these cells can promote angiogenesis, immunosuppression, and chemoresistance¹⁵.TGF-β(Transforming growth factor beta)driven cancer-associated fibroblasts (CAFs), mesothelial cells, and lymphatic endothelial cells are associated with a poor outcome, while plasma cells correlate with more favorable outcomes in HGSOC⁷.

Understanding the evolutionary trajectories of tumor and TME cells is crucial to develope more effective and personalized treatment strategies¹⁶. The multifaceted intratumoral heterogeneity and the TME may shape the evolutionary dynamics of OC¹⁷. Understanding the tumoral heterogeneity and evolutionary biodiversity could provide conceptual knowledge about the evolving tumor cell community to improve patient outcomes.

Preprocessing of single-cell RNA-seq data

We collected single-cell RNAseq matrix from five HGSOC datasets (GSE140819, GSE154600, GSE165897, GSE184880 and data in http://blueprint.lambrechtslab.org, Supplementary Table S1). Raw gene expression matrices were processed using the Seurat R package (version 3.2.3)¹⁸. Cells of low quality were filtered out based on two criteria: (1) cells with fewer than 1000 unique molecular identifiers (UMIs) or with fewer than 100 genes detected; (2) cells with more than 30% of their UMIs originating from mitochondrial genes. For each matrix, gene expression matrices were normalized using the LogNormalize method via the NormalizeData function, and highly variable genes were identified using the scran package¹⁹. Dimensionality reduction was performed with principal component analysis (PCA), where the number of principal components (PCs) was selected based on the JackStraw function(Seurat). Cell clusters were identified using the FindClusters function(Seurat)and visualized through uniform manifold approximation and projection (UMAP). Differential genes between clusters were calculated by the Wilcoxon rank-sum test with a family-wise error rate set at 5%.

Identification of tumor cells

Copy number variations (CNVs) were estimated in each cell by analyzing the averaged expression profiles across chromosomal intervals. The initial estimation of CNVs for each chromosomal region was conducted using the infercnv R package (version 1.0.4; inferCNV of the Trinity CTAT Project. https://github.com/broadinstitute/inferCNV). In epithelial cells, CNVs were determined by comparing the expression levels derived from the single-cell RNA-sequencing (scRNA-seq) dataset, with a cutoff value of 1 and a noise filter set at 0.2. For each sample, gene expression data were re-standardized, and the values were constrained within the range of − 1 to 1. The CNV for each cell was calculated as the sum of the squares of the CNV values across the respective regions.

Calculation of epithelial cell scores

We used the expression levels of epithelial scores based on the mean expression of five well-established markers: KRT19, KRT7, KRT18, KRT8 and EPCAM, as reported in previous studies²⁰.

Tumor evolution analysis

Selection of representative tumor cells

We input the clustering information from Seurat and simultaneously identified the tumor cell subgroups using the CNV method. We used the normalized matrix to calculate the PC scores for each cell. Subsequently, we used the Wilcoxon rank-sum test to obtain PCs that showed significant differences between classes (FDR < 0.05, false discovery rate). Finally, we selected the top 10 genes with the highest absolute coefficients of the differentially expressed PCs as the characteristic genes for the corresponding PCs.

We used the characteristic genes obtained from the tumor subgroups and applied the term frequency-inverse document frequency (TF-IDF) method. This technique is typically used in information retrieval and text mining to evaluate how important a word is to a document in a collection or corpus. In the context of gene expression analysis, TF-IDF (Term Frequency-Inverse Document Frequency) can be used to identify the most relevant genes in the tumor subgroups. TF-IDF measures the importance of a gene based on how frequently it appears in a particular subgroup compared to its occurrence across all subgroups in the dataset, thus highlighting genes that are both frequent and unique within specific tumor subgroups.

Where m and n represent the number of genes and cells in the single-cell expression matrix, respectively. \({N}_{i,j}\) is the \(m\) by \(n\) normalized gene expression matrix defined above. Within the TF-IDF method, we set the lower quartile of a gene’s expression across all cells as the threshold. If a gene’s expression was below this threshold, we set its expression to 0. This approach ensures that only genes with a certain level of expression are considered significant, effectively filtering out genes that are expressed at very low levels across the cells. After applying the TF-IDF method to obtain the gene weight matrix \({S}_{i,j}\), we calculated the 95th percentile (A) of each gene’s weight across all cells and set the lower threshold for gene weights at 0.25 × A. Any gene weight below this threshold in \({S}_{i,j}\) was set to 0, resulting in a corrected \({S}_{i,j}\) matrix. Subsequently, for each cell group, we summed the weights of the characteristic genes to derive a total weight and set the cutoff value \(\text{T}\) to establish the cell weight threshold utilizing the formula:

where Q1, Q2, and Q3 represent lower, median, and upper quartile values, respectively. Cells with a score higher than T were identified as representative tumor cells.

Clustering and refinement

Using the high-purity tumor cells obtained in the first step for each tumor subgroup, we employed consensus clustering method (SC3²¹) for consensus clustering to obtain subclassification labels for each subgroup. Subsequently, we merged the subcategories of each subgroup to identify those with statistically significant differences. Specifically, we first extracted the top 1000 highly variable genes from the PCA of high-purity malignant cells. For each subgroup’s subclassification obtained through the consensus clustering method SC3²¹, we trained a support vector machine model using the first 10 PCs and validated the classification accuracy through 10-fold cross-validation. Next, we randomly shuffled the labels 100 times to perform a permutation test to calculate the classification significance P-values between each pair of clusters. These P-values were then sorted in descending order, and an iterative merging process was conducted until the largest P-value was less than 0.05, at which point the iteration ceased. Through this method, we identified the subclusters within each tumor subgroup that were of the highest purity and had statistically meaningful classifications.

Hierarchical evolutionary analysis

We intersected the features derived from TF-IDF with those obtained through PCA and obtained tumor-associated feature genes. Using these genes, high-purity tumor cells, and the purified subpopulations of each subgroup, we computed the mean expression values of the corresponding subgroup’s feature genes to generate a pseudo-bulk dataset. We then applied the removeBatchEffect function in the Limma package (version 3.58.1) to remove batch effects across various datasets. Subsequently, we performed hierarchical evolutionary analysis between subgroups using the pvclust function from the pvclust package (version 2.2, nboot = 1000).

Following the aforementioned step, we obtained highly purified subpopulations within tumors, facilitating the analysis of intratumoral heterogeneity.

Pseudotime analysis

For each tumor, we generated a single-cell trajectory using the reverse graph embedding method offered by the Monocle R package (version 2.14.0). This method relies on the CellDataSet object to house single-cell gene expression data and the results of various analyses. Consequently, we established a CellDataSet for the single cells, setting the expression model parameter to “Family” with the value “negbinomial”. The construction of the single-cell trajectory entailed two primary phases. In the initial phase, our objective was to identify genes that play a pivotal role in defining the trajectory. For this purpose, we performed a differential gene expression analysis across different cell types. Genes exhibiting q-values of less than 0.001 were deemed significant for the pseudo-time analysis. The subsequent phase involved reducing dimensionality and constructing the trajectory with the identified genes. We employed the reverse graph embedding technique, which projects cells into a lower-dimensional space while learning a smooth, tree-like manifold. This step also includes the assignment of cells to this manifold, effectively mapping their progression or trajectory.

Cell interaction

We applied R package Cellchat (v 1.1.1) to evaluate the interaction weights among clusters, implemented with default parameters.

Cell-type enrichment analysis

Each cell-subtype analysis encompassed a range of tumor stages, for which we computed enrichment scores (EScores) to quantify the cell-subtype prevalence at different stages. These EScores reflect the ratio of cell subtype numbers at specific stages to their overall distribution and highlight when a subtype is predominantly enriched, with values greater than 1 signifying enrichment at that stage²⁰.

GSEA enrichment analysis

In this study, we utilized the `clusterProfiler` package (version 4.10) to perform gene set enrichment analysis (GSEA). First, genes were ranked based on their expression levels. We then employed the `GSEA` function to analyze the enrichment of predefined gene sets from the MSigDB database. The enrichment score (ES) and normalized enrichment score (NES) were calculated to assess the enrichment levels of gene sets, and the False Discovery Rate (FDR) was used to correct for multiple comparisons. To visualize pathway enrichment, we used the `enrichmentNetwork` function to create an enrichment network plot, which illustrates the relationships between the enriched pathways and the associated genes.

Weighted correlation network analysis in cell subclusters and acquisition of fibroblast-related genes

The normalized expression matrix was used to construct a weighted gene co-expression network via the weighted correlation network analysis (WGCNA) R package (version 1.69). To mitigate the effect of noise and outliers, analysis was performed on “pseudo cells,” which represent the average gene expression of 10 randomly selected cells within each distinct cell type²². Network construction was achieved using the “blockwiseModules” function, applying the default settings. For each identified module, a PCA was performed using the module eigengenes. The correlation between module eigengenes and cell-type metadata was calculated to evaluate the relevance of each module using Pearson’s correlation test. Subsequently, hub genes within significant modules were identified based on their modular connectivity, which refers to the absolute value of Pearson’s correlation between genes (module membership) and their relationship with clinical traits, defined as the absolute value of Pearson’s correlation between individual gene expression and cell type. We performed WGCNA in the fibroblast subpopulations of datasets GSE154600 and GSE165897 from GEO database²³ and took the intersection of characteristic genes of cell subpopulations significantly infiltrating the BR3 branch as the final fibroblast-related characteristic genes, which were used for the construction of the risk prognosis and drug resistance analysis model below.

Similarity analysis among cell subpopulations

We used the single-cell subpopulations and the characteristic genes obtained through WGCNA to calculate the mean expression levels of these genes in each subpopulation and sample phenotype (the evolutionary branch to which the sample belongs). We used the R package ggcor (https://github.com/hannet91/ggcor) to compute the correlations between cell subpopulations and between cell subpopulations and phenotypes and visualize these correlations.

Clustering of bulk RNA data samples

Bulk RNA data were retrieved from the Cancer Genome Atlas (TCGA; https://www.cancer.gov/tcga) and GEO datasets. Subsequently, non-negative matrix factorization (NMF) clustering methods were performed on the normalized expression data using the NMF R package (version 0.23).

Mapping of single-cell and bulk RNA subpopulations

To establish the correspondence between the intertumoral clusters obtained from single-cell evolutionary analysis and the clusters identified through clustering methods in bulk, we first extracted the characteristic gene sets from the single-cell BR1, BR2, and BR3 subpopulations. Subsequently, we used a hypergeometric test to calculate the enrichment scores of the bulk RNA-seq samples within these three gene sets. We then selected the connections with the highest enrichment scores that were significantly enriched as the mapping relationship between single-cell and bulk RNA-seq samples.

Friends analysis

The Friends analysis approach assesses the functional correlation among various genes within a pathway, suggesting that the interaction of a gene with others in the same pathway enhances its likelihood of expression. Using the R package GOSemSim²⁴, we calculated the functional correlations among genes linked to the prognosis of HGSOC and drug resistance.

Spatial transcriptome data analysis

We collected two HGSCO spatial transcriptome datasets including chemotherapy treatment response (GSE189843 and GSE211956, Supplementary Table S1) from GEO database. Both datasets are spatial transcriptomics obtained using the 10x Genomics platform, and they include matrices of spatial gene expression information where rows represent genes and columns represent spatial spots. Raw gene-spot matrices were analyzed with the Seurat package (version 3.2.3) in R. Spatial transcriptome data were qualitatively controlled using parameters including total spots, media UMIs/spot, median genes/spot, and median mitochondrial genes/spot. Spots used in the subsequent analysis were filtered for a minimum detected gene count of 200 genes while genes expressed in fewer than three spots were removed. Normalization across spots was performed with the SCTransform function. Dimensionality reduction and clustering were performed with PCA at a resolution of 1 with the first 30 PCs. We conducted cluster analysis using FindClusters and then used the standardized expression matrix to calculate the average expression levels of immune-related genes (PTPRC, CD2, CD3D, CD3E, CD3G, CD5, CD7, CD79A, MS4A1, and CD19)²⁵. The subpopulation with the highest expression levels was selected as the normal control. We used the InferCNV method to identify subpopulations of tumor cells. Subsequently, using the RegionNeighbours function from the R package STutility (version 1.1.1; https://ludvigla.github.io/STUtility_web_site/), we determined the cells at the tumor edge. We then used the FindMarkers function to calculate the differentially expressed genes in the region adjacent to the tumor edge.

Immune infiltration analysis

We assessed the immune score of various immune cells in HGSOC patients by employing xCell (R package, version 1.1) on RNA-seq datasets and microarray datasets²⁶. The microarray datasets underwent quantile normalization, while the RNA-seq dataset was quantified in terms of fragments per kilobase million(FPKM).

Cell transfection

To further validate the impact of the CXCL12-CXCR4 signaling pathway on cell viability, we performed siRNA knockdown of CXCR4 in HGSOC cell line SKOV3. The SKOV3 purchased from the ATCC cell bank, were digested with trypsin and resuspended. After cell counting, cells were plated at a density of 4 × 104 per well in a 24-well plate and incubated overnight in a cell culture incubator. siRNA transfection was performed according to the Lipofectamine 3000 (L3000015, Thermo Fisher Scientific (Waltham, MA, USA)) protocol and using 15 pmol siRNA and 1.5 µL lipofectamine per well. After 48 h of incubation, cells were collected for qRT-PCR and Western blot. The sequences of siRNA for the knockdown and control groups were as follows:

• siCXCR4-1 GGCAAUGGAUUGGUCAUCCUGGUCA;

• siCXCR4-2 UGGUUGGCCUUAUCCUGCCUGGUAU;

• siCXCR4-3 UGUUUCCACUGAGUCUGAGUCUUCA;

• siNC UUCUCCGAACGUGUCACGUTT.

qRT-PCR

qRT-PCR was used to detect the transcriptional level of the CXCR4 gene and the knockdown effect of different siRNA sequences. RNAiso Plus (9109, TaKaRa, Tokyo, Japan) was added to the cell samples for lysis, and the samples were processed according to the Chloroform-Trizol RNA Extraction protocol After air-drying at room temperature, RNA pellets were dissolved in 20 µL Rnase-free H₂O. The concentration and A260/280 of RNA samples were measured using a Microplate UV-Vis Spectrophotometer. Reverse transcription was performed according to the Goldenstar^® RT6 cDNA Synthesis Kit Ver.2 (TSK302S, Tsingke, Beijing, China) instructions, and 2 µg of total RNA was used as the template. The reaction mixture was incubated at 50 °C for 5 min and 85 °C for 2 min. The qPCR system was prepared according to the ArtiCanATM SYBR qPCR Mix (TSE501, instructions, Tsingke, Beijing, China) and the real-time PCR program was set as follows: holding stage, 95 °C 30 s; cycling stage (40 cycles), step 1 95 °C 15 s, step 2 60 °C 30 s; and melt curve stage: 95 °C 15 s, 60 °C 60 s, 95 °C 15 s. The primer sequences used were as follows:

• CXCR4-h-F ACTACACCGAGGAAATGGGCT;

• CXCR4-h-R CCCACAATGCCAGTTAAGAAGA;

• GAPDH-h-F TGACAACTTTGGTATCGTGGAAGG;

• GAPDH-h-R AGGCAGGGATGATGTTCTGGAGAG.

Western blot

Western blotting was used to detect the protein level of CXCR4. Cell samples were treated with RIPA (Radio-Immunoprecipitation Assay) lysis buffer (R0278, Sigma-Aldrich, St.Louis, MO, USA) containing Phenyl-methane-sulfonyl- fluoride (PMSF) and protease inhibitors; sonicated on ice; and centrifuged at 4 °C, 12,000 g for 10 min. Then, the supernatant was taken for protein quantification using the BCA (Bicinchoninic Acid Assay) method. After mixing with loading buffer and denaturing at 100 °C for 5 min, proteins were separated using 10% SDS-PAGE gel and then transferred to a PVDF (Polyvinylidene fluoride) membrane. After blocking with 5% skimmed milk for 1 h at room temperature, the membrane was incubated with CXCR4 (E3Q4B) Rabbit mAb (Cell Signaling Technology, Danvers, MA, USA) diluted at a ratio of 1:1000 in TBST (Tris-Buffered Saline with Tween) containing 5% BSA (Bovine Serum Albumin). After washing with TBST three times for 5 min each, incubate the membrane with Anti-rabbit IgG HRP-linked Antibody (7074 S, Cell Signaling Technology, Danvers, MA, USA) at a ratio of 1:2000 in TBST .Visualize the protein bands using ECL reagent (E-IR-R307, Elabscience Biotechnology, Wuhan, China) according to the manufacturer’s instructions.

CCK8 assay

The CCK8 assay was used to detect the effect of CXCL2/CXCR4 on cell viability. When cells reached 80–90% confluence, cells were digested with trypsin and resuspended. After cell counting, cells were plated at a density of 3000 cells/well in 96-well plates. After incubation in a 37 °C incubator for 24 h, the medium was replaced with complete medium containing 0, 50, 100, 200, or 300 ng/mL CXCL12 and cultured for 48 h. CCK8 reagent (10 µL) was added per well, and the absorbance at 450 nm was measured using a microplate reader after a 1-hour incubation at 37 °C. In separate experiments, siCXCR4 and siNC transfections were performed 24 h after plating. Six hours later, CXCL12 (200 ng/mL) was added, and the CCK8 assay was conducted using the same method.

Construction of a multi-omics atlas of HGSOC including chemotherapy response characteristics using single-cell and spatial transcriptomics.

To enhance the understanding of the complex interplay between intra- and intertumoral heterogeneity in HGSOC, and to develop a comprehensive tumor evolutionary model that delineates the distinct tumor and microenvironmental profiles across diverse patient cohorts, we assembled a robust dataset (Fig. 1A). This dataset comprises four single-cell datasets and two spatial transcriptomics datasets from GEO datasets (Supplementary Table S1), totaling 50 samples. Notably, 28 (57%) of these samples were obtained from patients prior to chemotherapy and who exhibited measurable responses to therapy (Fig. 1B). We also curated two bulk RNA-sequencing datasets that include chemotherapy treatment responses from GEO datasets, encompassing 24 samples (Supplementary Table S1). Our analytical approach began with the application of advanced text-mining and machine-learning algorithms to three of the single-cell datasets. This enabled the construction of a sophisticated hierarchical evolutionary model specific to HGSOC. This model was then validated and further scrutinized for insights into the relationships between intratumoral microenvironmental heterogeneity and intertumoral population heterogeneity using the remaining two datasets. Subsequently, we focused on the spatial transcriptomics datasets to dissect the spatial distribution characteristics of microenvironmental heterogeneity In the merged test dataset, which consisted of 14 samples, we processed and analyzed a total of 87,288 cells after quality control (Fig. 1C). These cells were classified into several key subtypes: B cells (5,700, 6.53%), CAFs (17,664, 20.24%), endothelial cells (ECs; 1,174, 1.35%), epithelial cells (Epi; 20,273, 23.22%), monocytes (15,217, 17.43%), and T and NK cells (T_NK; 27,259, 31.23%, Fig. 1C,D). We used the CNV method and calculated the epithelial cell score to identify 18,273 individual tumor cells (Fig. 1E). Through this integrated approach, we aimed to not only elucidate the evolutionary dynamics of HGSOC but also provide a foundation for the development of personalized therapeutic strategies that account for the unique heterogeneity observed in each patient’s tumor.

High-grade serous ovarian cancer transcriptome atlas. (A) Schematic depicting the study design. (B) Number of samples in the tumor evolution analysis of high-grade serous ovarian cancer (HGSOC). Pie chart showing the proportions of clinical treatments in the tumor evolution analysis. Number of cells and spots in the scRNA-seq datasets. (C) The t-distributed stochastic neighbor embedding (t-SNE) plots showing the major cell types in HGSOC. Clusters are distinguished by color. (D) Heatmap showing cell-type marker gene expression level in the first single-cell dataset. (E) Expression profile of epithelial cells and tumor scores in the first single-cell dataset. The colors from gray to red represent the expression level from low to high.

Full size image

Inter- and intratumoral heterogeneity in high-grade serous ovarian cancer

Patient prognosis and treatment response are significantly affected by inter- and intratumoral heterogeneity. We used text mining and consensus clustering to dissect the clonal evolution within individual patients, followed by hierarchical clustering and bootstrap analysis to establish intertumoral clonal evolutionary relationships. We identified three evolutionary branches, designated as BR1, BR2, and BR3, across 69 clusters in 14 samples from GSE140819, GSE154600, and GSE184880 (Fig. 2A). Notably, individual tumor samples tended to cluster under a specific evolutionary branch. To validate the robustness of our method, we replicated the analysis in a subset of the data and confirmed consistent evolutionary relationships among cell populations in dataset GSE184880 (Supplementary Figure S1).

Tumor branching evolution reveals intratumor heterogeneity in high-grade serous ovarian cancer. (A) Tumor phylogenetic tree constructed by hierarchical clustering of all the clusters from 14 tumors, in which BR1, BR2, and BR3 were defined according to the hierarchical relationship. (B) Bar plot showing enrichment analysis using the tumor branch evolution features via clusterProfiler. FOR: Formation, PRO: proliferation, ORG: organization, POS REG: positive regulation. (C) Bar plot showing the sample origins of three subtypes of branching evolution. (D) Distribution characteristics of intratumoral cell types obtained through tumor evolutionary analysis. Profile and uniform manifold approximation and projection (UMAP) plots showing the cell-type subgroups in the epithelial, CAF, macrophage, and CD8 cells (top). Velocity and single-cell trajectory results (rows 2 and 3) and differentially expressed genes (rows) along the pseudo-time (columns) were clustered hierarchically into five groups in the scRNA-seq dataset. Pathway enrichment scores were calculated using clusterProfiler.

Full size image

Having established the clonal evolutionary relationships across different samples, we identified key genes within these branches and conducted functional analyses (Supplementary Table S2). Key genes of BR1 were enriched in pathways related to myeloid leukocyte migration (GO:0097529), neutrophil degranulation (GO:0043312), and digestion (GO:0007586). BR2 genes were predominantly associated with cellular responses to chemokines (GO:1990869), T-cell migration (GO:0072678), regulation of response to type II interferon (GO:0060330), and autophagy (GO:0006914). BR3 genes were enriched in pathways such as the generation of amyloid fibrils (GO:1990000), cellular response to hypoxia (GO:0071456), positive regulation of growth (GO:0045927), and collagen fibril organization (GO:0030199) (Fig. 2B). These findings suggest that BR3 is characterized by highly fibrotic, immune desert-like tumors, BR2 by T-cell-infiltrated hot tumors, and BR1 by predominantly tumor-cell-dominated cold tumors²⁷.

Further analysis revealed significant heterogeneity in the distribution of different microenvironmental cells across the branches, with CAFs notably enriched in the BR3 branch (Fig. 2C). To delve deeper into intratumoral heterogeneity, we dissected the heterogeneity of epithelial, fibroblast, monocyte, and CD8 cells. We identified six epithelial cell subpopulations (E_HSPA1B, E_JUND, E_MALAT1, E_S100A4, E_SPP1, and E_TOP2A), with E_SPP1 and E_TOP2A predominantly found in the BR3 branch (Fig. 2D). Secreted phosphoprotein 1(SPP1) is highly expressed in multiple tumor types and interacts with fibroblasts or T cells via pathways such as SPP1-CD44, correlating with tumor malignancy²⁸. Type IIA topoisomerase(TOP2A), a cell cycle gene, indicates high proliferative states when overexpressed. Cell velocity and pseudo-time analysis revealed that epithelial cells evolve toward the BR3 subpopulation, with terminal pathways enriched in the cell cycle, epithelial differentiation, and lipid metabolism pathways (Fig. 2D).

We identified five fibroblast subpopulations (F_ACTA2, F_COL1A1, F_JUNB, F_LGALS3, and F_MMP11), most of which were concentrated in the BR3 branch (Fig. 2D). Cell velocity analysis showed enrichment of cell adhesion-related pathways at the onset and metabolic, chemokine, and MAPK pathways at the end (Fig. 2D). Eight macrophage subpopulations were identified (M_A2M, M_CCL3, M_CXCL10, M_HLA-DPB1, M_HMGB2, M_HSPA1B, M_MMP9, and M_SPP1), with M_CCL3, M_CXCL10, M_HLA-DPB1, and M_SPP1 enriched in the BR3 branch (Fig. 2D). Temporal analysis revealed the enrichment of lipid metabolism, granulocyte differentiation, and c-Jun N-terminal kinase (JNK) pathway pathways in the terminal stage (Fig. 2D). Four CD8 cell subpopulations were identified, with CD8_GZMA and CD8_TIGIT enriched in BR1 and BR2, and BR3 dominated by CD8_CXCR4. Cell velocity analysis showed enrichment of T-cell activation, lymphocyte differentiation, and autophagy regulation pathways in the terminal stage of cell differentiation (Fig. 2D).

These results underscore significant inter- and intratumoral cellular heterogeneity in HGSOC, with epithelial cells in BR3 samples exhibiting high proliferation rates and infiltration by diverse fibroblast subtypes.

Tumor functional clonality is associated with patient prognosis and TME composition

In our previous analysis, we divided all single-cell samples into three branches and identified significant heterogeneity between these branches, along with 150 characteristic genes for each branch (Supplementary Table S2). We further examined the prognostic clustering effect of these genes. Prognostic analysis across five datasets (TCGA, GSE14764, GSE26193, and GSE26712; Supplementary Table S1) revealed strong prognostic clustering in these datasets, with the poorest prognosis samples exhibiting high levels of fibroblast infiltration (G1 in TCGA and G2 in GSE14764, GSE26193, and GSE26712; Fig. 3A,B). Correlation analysis also indicated a significant high correlation between fibroblast content and fibroblast-associated marker genes (COL3A1, COL1A2, COL1A1, FN1, and DCN; Supplementary Figure S2). Subsequently, we used the obtained intertumoral heterogeneity branch genes as a reference gene set to perform a hypergeometric test in different bulk samples to calculate the distribution characteristics of samples in the single-cell BR branches. Bulk samples with the poorest prognosis and high fibroblast infiltration were significantly enriched in BR3 scores (Fig. 3C), further validating the reliability of using single-cell samples for studying intertumoral heterogeneity in tumors. Differential gene analysis in bulk samples also revealed the enrichment of numerous pathways related to epithelial cell migration, cell adhesion, and chemokines in samples with poor prognosis (Fig. 3D). These findings suggest that the fibroblast content in HGSOC samples may be highly correlated with tumor malignancy.

Tumor branching evolution reveals that intertumoral heterogeneity and the proportion of fibroblasts promote the poor prognosis of high-grade serous ovarian cancer. (A) Overall survival curves showing the prognosis of the three subtypes (G1, G2, and G3) obtained from non-negative matrix factorization (NMF) clustering using the 150 tumor evolution features in TCGA and GEO cohorts. (B) Boxplots showing the immune cell infiltrates ratio in the three distinct malignant subtypes in the significantly enriched patients via xCell (ns, not significant, *P < 0.05, **P < 0.01, and ***P < 0.001). Pairwise comparisons were conducted using the Wilcoxon rank-sum test in the RNA cohort. For the boxplot, the centerline represents the median, and the box limits represent the upper and lower quartiles. (C) Boxplot showing the GSVA enrichment scores in the poorest prognosis using the branch features of tumor evolution analysis in scRNA datasets. Boxplots showing the mean expression level of BR3 genes in the three subtypes of bulk RNA datasets. ns, not significant, *P < 0.05, **P < 0.01, and ***P < 0.001 by Wilcoxon rank-sum test. (D) GO enrichment analysis of upregulated genes of the poorest prognosis group (G1 in TCGA; G2 in the other cohorts). GO:0030198: Extracellular matrix organization, GO:0043062:extracellular structure organization, GO:0071559: Response to transforming growth factor beta, GO:0031589: Cell-substrate adhesion, GO:0060562: Epithelial tube morphogenesis, GO:0010631: Epithelial cell migration, GO:0045123: Cellular extravasation, GO:0061448: Connective tissue development, GO:0048660: Regulation of smooth muscle cell proliferation, GO:0050673: Epithelial cell proliferation, GO:0007043: Cell-cell junction assembly, GO:0032642: Regulation of chemokine production, GO:0007160: Cell-matrix adhesion, GO:0048639: Positive regulation of developmental growth, GO:0016049: Cell growth, GO:0001558: Regulation of cell growth, GO:0042692: Muscle cell differentiation, GO:0030336: Negative regulation of cell migration, GO:0045229: External encapsulating structure organization, GO:0090130: Tissue migration, GO:0032602: Chemokine production.

Full size image

Identification of a subpopulation of CXCL12-positive fibroblasts associated with chemotherapy resistance and poor prognosis

Previous research has demonstrated a strong correlation between the heterogeneity of TME composition, especially the infiltration of fibroblasts, and prognosis². Given that chemotherapy is commonly used in the clinical treatment of ovarian cancer, we explored the relationship between intratumoral heterogeneity and chemotherapy resistance. We analyzed the differences in the cellular composition of the TME in tumors with varying responses to chemotherapy. Using the 150 characteristic genes obtained (Supplementary Table S2), we conducted a tumor evolution analysis on five samples within the single-cell dataset (GSE154600). The chemotherapy-sensitive samples GSM4675276 and GSM4675276 belonged to the BR1 subgroup, while the chemotherapy-resistant samples GSM4675273 and GSM4675273 were classified under the BR3 subgroup (Fig. 4A; Supplementary Table S3). Combining these five samples, a total of 50,571 cells were obtained, including B cells, CAF, monocytes, T cells, EPI_1 (epithelial cells with high levels of Matrix metallopeptidase 7 (MMP7) and E74 like ETS transcription factor 3 (ELF3)), and EPI_2 (epithelial cells with a high level of hes family bHLH transcription factor 1 (HES1)) (Fig. 4B and Supplementary Figure S3). Among these, B cells were significantly enriched in samples derived from BR1 and BR2, while CAFs were significantly enriched in BR3 samples, indicating that the degree of fibroblast infiltration not only affects patient prognosis but also correlates with the treatment response (Fig. 4C). To elucidate the relationship between CAFs and ovarian cancer treatment response, we further sub-classified CAFs and identified seven distinct subtypes (Fig. 4D, Supplementary Figure S4). F_CCL21 and F_RAMP3 were significantly enriched in BR1, and F_PGF, F_IGFBP4, F_IGFBP7, and F_CXCL12 were significantly enriched in BR3, particularly the F_CXCL12 subtype, which was almost exclusively found in BR3 (Fig. 4E). We then used gene co-expression analysis to identify key characteristic genes in different CAF subpopulations, resulting in four gene modules, with the blue module showing the highest correlation with the BR3-specific F_CXCL12 subgroup (R = 0.88, P < 0.001, Fig. 4F) and the yellow module highly correlated with the BR1-specific F_RAMP3 module (R = 0.97, P < 0.001; Fig. 4F). Analysis of the hub genes in the blue module revealed that these genes were primarily enriched in pathways related to extracellular matrix remodeling and tube morphogenesis (Fig. 4G and Supplementary Figure S5). Given the significant enrichment of CAFs in BR3, we explored the communication between these fibroblasts and tumor cells and found high-intensity interactions between the fibroblast subpopulations F_CXCL12 and F_IGFBP7 and tumor cells in EPI_1 (Fig. 4H). Specifically, F_CXCL12 interacts with tumor epithelial cells via the ligand-receptor pair CXCL12/CXCR4 and through the secretion of placental growth factor (PGF) (Fig. 4I). These findings suggest that CXCL12-positive CAF cells may be associated with chemotherapy resistance.

Intra- and intertumoral heterogeneity of cancer-associated fibroblasts. (A) Tumor phylogenetic tree constructed by hierarchical clustering using the 150 branch genes. (B) UMAP plot showing the major cell types in dataset GSE154600. (C) Bar plot showing the origins of cell types in three subtypes of branching evolution. (D) UMAP plot showing the subtypes of cancer-associated fibroblasts (CAFs). (E) Bar plot showing the origins of CAFs in the three evolutionary subtypes. (F) WGCNA results showing the gene modules in distinct CAF subtypes. Columns represent cell types. The colors from blue to red indicate low to high correlation between the gene module and cell subtypes (Pearson correlation test). (G) GO enrichment analysis of hub genes of the BR3 enrichment subtype (F_CXCL12). (H) Number of significant ligand-receptor pairs between CAF and epithelial subtypes. The edge width is proportional to the indicated number of ligand-receptor pairs. EPI_1, epithelial subtype with high expression of MMP7 and ELF3; EPI_3, epithelial subtype with high expression of HES1 and CD24. (I) Dot plot showing the ligand-receptor pairs between CAFs and epithelial cells. Rows represent the ligand receptor (L-R) pairs, and columns represent cell subset–cell subset pairs. The color gradient from black/blue to red indicates the mean values of the L–R pairs from low to high, and the circle size indicates the significance of the pairs. P-values were calculated via a permutation test using CellChat.

Full size image

To validate the relationship between fibroblasts and chemotherapy, we examined the relationship between fibroblasts and treatment response in the other dataset (GSE165897) that included chemotherapy response information (Supplementary Table S4). Using the previously identified 150 characteristic genes, we performed a hierarchical evolutionary analysis between samples and conducted permutation tests. The platinum-free interval (PFI) values of samples of EOC3, EOC349, EOC540, and EOC87 in the BR3 subgroup were significantly lower than those in the groups BR1 and BR2, indicating significant chemotherapy resistance, while the samples EOC136 and EOC153 in BR1 showed chemotherapy sensitivity, again confirming the high correlation between intertumoral heterogeneity and treatment response (Fig. 5A; Supplementary Table S4). From this dataset of nine samples, we obtained 17,898 cells, comprising 11 subpopulations, including two fibroblast subpopulations, F_CXCL12 and F_COL1A1, with significant F_CXCL12 infiltration in BR3 and significant F_COL1A1 enrichment in BR1 (Fig. 5B,C and Supplementary Figure S6). We then re-clustered the fibroblasts and identified five fibroblast subpopulations, including the F_CXCL12 subgroup (Fig. 5D and Supplementary Figure S6). Using gene co-expression analysis, we also obtained a gene module that was highly correlated with the F_CXCL12 subgroup (R = 0.85, P < 0.001; Fig. 5E), with characteristic genes primarily enriched in TGF-beta, MAPK, and cytokine interaction signaling pathways (Supplementary Figure S7). We also verified the correlation between CXCL12-positive fibroblasts and drug sensitivity in the GM dataset (Fig. 5F and Supplementary Figure S8).

Heterogeneity of cancer-associated fibroblasts is associated with chemotherapy treatment outcomes. (A) Platinum-free interval values in the three tumor evolution branches in dataset GSE165897. (B) UMAP (i) and bar plot (ii) showing the major cell types and their origin. (C) Volcano plot showing the differential genes for cancer-associated fibroblast (CAF) subtypes. Upregulated genes are indicated in red, while downregulated ones are indicated in blue. The X-axis represents the log2 fold change, and the Y-axis represents the single-cell clustering information. (D) UMAP(i) and bar plot (ii) showing the subtypes and origins of CAFs. Scatter plot showing the expression level of marker genes (iii). (E) WGCNA results showing the gene modules of distinct CAF subtypes in GSE165897. (F) Heatmap showing the CAF subtype correlation between datasets GSE154600 and GSE165897.

Full size image

In these datasets, we identified a fibroblast subgroup with high CXCL12 expression, and a similarity correlation analysis between subgroups confirmed the high similarity of CXCL12-high fibroblast subgroups in both datasets (Fig. 5E).

Identification of a gene set affecting prognosis and drug resistance in HGSOC

Above, we identified a significant correlation between the highly invasive CXCL12-positive CAF subpopulation and both treatment response and prognosis in ovarian cancer. We next conducted a comparative analysis of the WGCNA hub gene results to consolidate the characteristic genes of the CXCL12-positive CAF subpopulation and identified 24 shared genes (Supplementary Table S5). Subsequent Cox regression analysis using the GSE26193 dataset revealed that DCN, CXCL12 and TNFAIP6 were significantly positively associated with poor prognosis in ovarian cancer, with DCN specifically identified as a fibroblast-characteristic gene (Fig. 6A). Using these significantly differential genes and Cox regression coefficients, we performed a risk analysis across four additional datasets (GSE14764, GSE26712, TCGA and GSE9891, Supplementary Table S1). The findings indicated that the risk coefficients for the high gene expression groups were notably higher than those for the low expression groups, and both CXCL12 and DCN showed significant positive correlations with poor overall patient prognosis (Fig. 6A–C). Through Friends analysis, we observed a high correlation of CXCL12 with other genes (Fig. 6D). Additionally, using these 24 genes, we applied enrichment analysis in three treatment datasets (GSE33482, GSE114206²⁹, and GSE189843¹⁵) to validate the association between CAFs and chemotherapy response (Fig. 6E; Supplementary Table S6). Based on gene expression data from these samples, we categorized them according to their enrichment levels in the 24 genes. These 24 genes could effectively stratify samples into high and low groups, with the high CAF subpopulation exhibiting strong drug resistance. Moreover, the 24 genes showed excellent predictive efficacy for drug resistance, with AUC values of 1, 0.89, 0.83 in datasets GSE33482, GSE114206 and GSE189843 respectively (Fig. 6E). These findings underscore the critical role of fibroblasts in OC treatment and prognosis, and through further feature extraction, we have identified a reduced set of genes that are indicative of treatment response.

Correlation between 24 genes of CXCL12-positive fibroblasts and prognosis and drug resistance in high-grade serous ovarian cancer. (A) Forest plot showing the risk prognosis results from 24 shared CXCL12-positive fibroblasts obtained from two single-cell samples via COX regression in GSE14764. (B) Overall survival curves showing the prognosis results with different cancer-associated fibroblast (CAF) risk scores in the four high-grade serous ovarian cancer cohorts. In this analysis, for Fig. 6B, we first identified 24 signature genes associated with CXCL12-positive CAF cells. We then applied the COX regression coefficients obtained from Fig. 6A to compute a risk score for each patient. Using the median risk score as the cutoff, we divided the patients into high-risk and low-risk groups. To assess the prognostic significance, we employed the log-rank test to calculate the p-value, thereby verifying that patients with high expression of CXCL12-positive CAF genes have poorer prognoses. (C) We performed a log-rank test to calculate the p-values and assess the correlation between genes identified in the dataset GSE14764 in Fig. 6A (DCN, CXCL12 and TNFAIP6) that are most strongly associated with poor prognosis in relation to CXCL12-positive CAFs. We conducted the prognostic analysis using 1,232 serous OV patients’ data available on the https://kmplot.com/analysis, by selecting the median expression level of each gene as the cutoff. (D) Box plot showing the Friends analysis results. (E) AUC and Sankey diagram showing the prediction of chemotherapy resistance using the 24 CXCL12 positive CAFs genes. R, chemotherapy resistant. S, chemotherapy sensitive.

Full size image

Spatial distribution characteristics of CXCL12-positive fibroblasts affect the clinical prognostic results of HGSOC

Given the strong correlation between CXCL12-positive fibroblasts and the treatment outcomes and prognoses of HGSOC. Through cell communication analysis, we identified that CXCL12-positive fibroblasts interact with tumor cells via CXCR4. Therefore, we elucidated the interaction mechanisms of CXCL12-positive fibroblasts in vitro and in spatial transcriptomic data. In vitro, we silenced the CXCL12 receptor gene CXCR4 in the ovarian cancer cell line SKOV3 (Fig. 7A). As shown by qPCR, all three siRNA knockdown fragments significantly suppressed CXCR4 mRNA expression levels in SKOV3 cells (P < 0.05; Fig. 7A). Subsequent addition of exogenous CXCL12 protein to both the control group and the CXCR4-si group demonstrated that CXCR4 knockdown significantly inhibited cell viability compared with the siNC control group (P < 0.05; Fig. 7B,C and Supplementary Figure S9). Considering that fibroblasts and tumor cells interact through cellular communication, and the intensity of this communication is highly related to the spatial positioning of genes, we analyzed the spatial distribution of fibroblasts and tumor cells within a treatment response-associated spatial transcriptomic dataset. Using the inferCNV method and cell type identified, we identified tumor regions and the surrounding cellular environments. We analyzed the spatial distribution characteristics of cells in chemotherapy-resistant and chemotherapy-sensitive samples in GSE211956³⁰ dataset respectively (Fig. 7D,E, Supplementary Table S7, Supplementary Figures S10, S11). In the chemotherapy-resistant sample (GSM6506110), we identified 5 clusters (Fig. 7D(i)). Cell type analysis revealed that cluster 4 is enriched with fibroblasts with high expression level both CXCL12 and DCN (Fig. 7D(iii)). Furthermore, in the surrounding area of cluster 4 (Nbs_4), which are enriched in tumor epithelial cells and expressed high level of tumor cell marker KRT19 and CXCR4 (the receptor for CXCL12) (Fig. 7D(iii)). In contrast, in the chemotherapy-sensitive sample (GSM6506114), we identified 4 clusters (Fig. 7E(i)). Among these, clusters 0, 1, and 3 are enriched with fibroblasts (Fig. 7E(ii)). Cluster 1 is specifically enriched with KRT19-positive tumor cells (Fig. 7E). However, in the regions dominated by KRT19-positive tumor cells, there is no high expression of CXCR4, forming a spatial interaction ___domain of CXCL12-positive CAFs. Finally, in our treatment response-associated clinical samples, we also validated the local spatial adjacency of CXCL12/DCN/KRT19 in chemotherapy-resistant samples (Fig. 7F).

Function and spatial distribution characteristics of CXCL12-positive fibroblasts. (A) The bar chart showing the silencing effect of the CXCL12 receptor gene CXCR4 (ns, not significant, *P < 0.05, **P < 0.01, and ***P < 0.001 by t-test). (B) Western blot showing the silencing effect of CXCR4 and the expression levels of CXCR4 protein in the control and CXCR4-silenced groups after the addition of exogenous CXCL12 protein. (C) The CCK8 results show that silencing CXCR4 significantly inhibits the viability of tumor cells (*P < 0.05 by t-test). (i) Clustering and spatial distribution, (ii) cell type composition in each cluster, and (iii) gene profile around the tumor boundary in (D) chemotherapy-resistant samples and (E) chemotherapy-sensitive samples. (F) The multiplex immunofluorescence results show the spatial proximity relationship between fibroblasts and tumor cells in chemotherapy-resistant (i) and chemotherapy-sensitive (ii) and samples. Scale bars = 50 μm.

Full size image

Solid-organ malignancies exhibit heterogeneity that enhances tumor cell survival and drug resistance through various molecular mechanisms³¹. Intratumoral heterogeneity evolves both spatially and temporally during tumor development, reprogramming the tumor microenvironment³². Understanding tumor clonality and its evolutionary path is crucial for comprehending tumor cell behavior and identifying driving factors^19,33. Tumors harbor numerous genetic and epigenetic alterations, many of which complicate the interpretation of their functional effects on tumor evolution²¹. Single-cell and spatial transcriptomics can delineate the roles of different cell types within tumors, providing insights into the functional diversity arising from genetic heterogeneity and adaptation^22,23,24,26. In this study, we integrated single-cell, bulk, and spatial transcriptomic data; elucidated the intratumoral heterogeneity (including tumor and TME cells of HGSOC); explored their evolutionary patterns of intra- and intertumoral heterogeneity; and identified critical factors influencing tumor prognosis and therapeutic responses.

Single-cell data processing faces challenges with large-scale datasets, necessitating efficient computational and statistical methods to handle data sparsity and ensure informative analysis³⁴. Batch effects between datasets require robust correction techniques to guarantee accurate cross-dataset comparisons and integrations³⁴. Determining the optimal number of clusters is essential for the analysis of intratumoral heterogeneity, yet this process is complicated by the lack of clear standards³⁵. In this study, we initially identified intratumoral heterogeneity using results from Seurat’s initial clustering combined with CNV data. We extracted tumor-associated feature genes and identified “representative cells” that reflect the characteristics of different tumor clones from subpopulations using the TF-IDF method of text mining. We then analyzed the intratumoral heterogeneity of these cells. This approach ensures the extraction of the most accurate and essential tumor clone information while minimizing computational resource consumption. Additionally, it mitigates the effect of sparsity in single-cell transcriptome data on data analysis. Furthermore, we used the representative cells to determine the clustering numbers automatically and therefore confirm the intratumoral heterogeneity using SC3 consensus clustering and support vector machine methods. This strategy addresses the subjectivity associated with the manual definition of cluster numbers. Finally, by obtaining precise intratumoral heterogeneity and related feature genes, we constructed a tumor evolution model based on their expression levels and analyzed the characteristics of intra- and intertumoral heterogeneity in HGSOC. Through the integrated analysis of intra- and intertumoral heterogeneity, we can more comprehensively elucidate the roles of tumor and microenvironmental cells in the prediction of prognosis and chemoresistance of HGSOC. This provides a novel strategy for interpreting the correlation between intra- and intertumoral heterogeneity using expression data.

Using the aforementioned strategy, we analyzed and validated across multiple datasets that HGSOC patients can be classified into three evolutionary branches: BR1, BR2, and BR3. Notable heterogeneity was observed within each branch concerning the quantity and composition of tumor and microenvironmental cells. The BR1 branch was predominated by tumor cells and resembled cold tumors, BR2 was infiltrated with high level of B and T cells and resembled hot tumors, and BR3 was heavily infiltrated by fibroblasts, akin to immunologically excluded tumors²⁷. By correlating tumor evolution with prognosis and resistance to therapy, we found that a high infiltration of fibroblasts was associated with a tumor cell subpopulation enriched in cell cycle pathways, low CD8 cell activity, increased cell proliferation, and malignant transformation in cell morphology, resulting in significantly worse prognosis and therapy resistance. Fibroblasts can be broadly classified into myofibroblasts, TGF-β-driven CAFs, inflammatory fibroblasts, and antigen-presenting fibroblasts and are highly associated with poor prognosis in various cancers^20,36,37. In HGSOC, fibroblasts can form a dense physical barrier around tumor cells, impede the infiltration of immune cells, and promote tumor cell proliferation through the secretion of cytokines^38,39. Combining with WGCNA method, we identified a clinical predictive model comprising 24 genes and discovered that fibroblasts with high CXCL12 expression are significantly related to poor prognosis and chemotherapy resistance. Our study demonstrates that CXCL12 exerts its biological functions through its receptor CXCR4, significantly affecting cell viability. The siRNA-mediated knockdown of CXCR4 further confirmed its pivotal role in this process. The reduction in cell viability upon CXCR4 knockdown aligns with previous research findings, reinforcing the central importance of CXCL12 in cell survival via CXCR4. The differences in the spatial distribution of cells, particularly the characteristics and functions of cells near the tumor boundary, play a crucial role in the development of tumor heterogeneity²⁵. To this end, we developed a bioinformatics analytical workflow for tumor boundary identification and found that this type of fibroblast significantly infiltrates the margins of chemotherapy-resistant tumors. Spatial interactions between cell clusters may have a more profound effect on chemo responsiveness than the composition of the clusters alone. In HGSOC, fibroblasts can interact with other cell types through various pathways, ultimately contributing to chemotherapy resistance³⁷. CXCL12-positive CAFs could promote cancer cell migration and invasion and upregulate the expression of PDL1 in bladder and pancreatic cancer cells^40,41. They can also attract monocytes through the CXCL12-CXCR4 pathway and induce their differentiation into M2 macrophages, which leads to enhanced tumor cell proliferation and reduced apoptosis in oral squamous cell carcinoma³². It has also demonstrated that in mycosis fungoides (MF), CAFs promote MF cell survival and migration through the secretion of CXCL12. Targeting the CXCL12/CXCR4 axis may enhance the effectiveness of anticancer therapies by disrupting the protective tumor microenvironment created by CAFs⁴². CXCL12 produced by CAFs and regulates cancer stem cell phenotypes could promote chemoresistance through CXCR4 interaction and activation of Wnt/β-catenin or PI3K/Akt pathways in multiple tumors⁴³. In this study, by analyzing the evolutionary pathways of tumor and microenvironmental cells, we determined that CXCL12-expressing fibroblasts, through spatially proximal interactions with tumor cells, influence patient prognosis and therapeutic outcomes, offering a new perspective for the prevention and treatment of HGSOC.