Decoding the effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer through tumor evolution analysis

Decoding the effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer through tumor evolution analysis | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Decoding the effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer through tumor evolution analysis tingjie wang This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4827560/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Tumor heterogeneity is associated with poor prognosis and drug resistance, leading to therapeutic failure. Here, we used tumor evolution analysis to determine the intra- and intertumoral heterogeneity of high-grade serous ovarian cancer (HGSOC) and analyze the correlation between tumor heterogeneity and prognosis, as well as chemotherapy response, through single-cell and spatial transcriptomic analysis. We collected and curated 28 HGSOC patients’ single-cell transcriptomic data from five datasets. Then, we developed a novel text-mining-based machine-learning approach to deconstruct the evolutionary patterns of tumor cell functions. We then identified key tumor-related genes within different evolutionary branches, characterized the microenvironmental cell compositions that various functional tumor cells depend on, and analyzed the intra- and intertumoral heterogeneity as well as the tumor microenvironments. These analyses were conducted in relation to the prognosis and chemotherapy response in HGSOC patients. We validated our findings in two spatial and seven bulk transcriptomic datasets (total: 1,030 patients). Using transcriptomic clusters as proxies for functional clonality, we identified a significant increase in tumor cell state heterogeneity that was strongly correlated with patient prognosis and treatment response. Furthermore, increased intra- and intertumoral functional clonality was associated with the characteristics of cancer-associated fibroblasts (CAFs). The spatial proximity between CXCL12-positive CAFs and tumor cells, mediated through the CXCL12/CXCR4 interaction, was highly positively correlated with poor prognosis and chemotherapy resistance in HGSOC. In this study, we constructed a panel of 24 genes through statistical modeling that correlate with CXCL12-positive fibroblasts and can predict both prognosis and the response to chemotherapy in HGSOC patients. Tumor cell state Functional clonality Tumor evolution HGSOC Tumor transcriptomic heterogeneity Microenvironment Spatiotemporal transcriptome CXCL12 Cancer-associated fibroblasts Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 1. Introduction Ovarian cancer (OC) ranks first in mortality among gynecologic malignancies [ 1 ]. High-grade serous ovarian cancer (HGSOC) is the most common and deadly subtype of OC [ 2 ]. Over 75% of HGSOC patients are diagnosed at an advanced stage with extensive malignant ascites and omental metastasis [ 3 ]. Complete staging surgery followed by platinum-based chemotherapy is the standard treatment for HGSOC [ 4 ]. 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 [ 5 ]. 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 [ 6 ]. 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) [ 6 ]. 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[ 3 ]. In HGSOC, it could be classified into four molecular subtypes: mesenchymal, immunoreactive, differentiated, and proliferative [ 7 ]. Among these, the immunoreactive subtype demonstrates a better prognosis and treatment response [ 7 ]. 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 [ 8 ]. However, intratumoral heterogeneity causes drug resistance, leading to therapeutic failure [ 6 ]. 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 [ 9 ]. Natural selection drives clonal expansion, leading to the emergence of subclones with varying proliferative, migratory, and invasive capabilities [ 10 ]. 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 – 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 [ 13 ]. 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 [ 15 ]. 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 [ 7 ]. Understanding the evolutionary trajectories of tumor and TME cells is crucial to develope more effective and personalized treatment strategies [ 16 ]. The multifaceted intratumoral heterogeneity and the TME may shape the evolutionary dynamics of OC [ 17 ]. Understanding the tumoral heterogeneity and evolutionary biodiversity could provide conceptual knowledge about the evolving tumor cell community to improve patient outcomes. 2. Results 2.1 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. 1 A). This dataset comprises four single-cell datasets and two spatial transcriptomics datasets from GEO datasets (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. 1 B). We also curated two bulk RNA-sequencing datasets that include chemotherapy treatment responses from GEO datasets, encompassing 24 samples (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 intra-tumoral 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. 1 C). 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%). We used the CNV method and calculated the epithelial cell score to identify 18,273 individual tumor cells (Fig. 1 C). 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. 2.2 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. 2 A). 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 (Figure S1). Having established the clonal evolutionary relationships across different samples, we identified key genes within these branches and conducted functional analyses (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. 2 B). 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 [ 18 ]. Further analysis revealed significant heterogeneity in the distribution of different microenvironmental cells across the branches, with CAFs notably enriched in the BR3 branch (Fig. 2 C). 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. 2 D). 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[ 19 ]. 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. 2 D). 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. 2 D). Cell velocity analysis showed enrichment of cell adhesion-related pathways at the onset and metabolic, chemokine, and MAPK pathways at the end (Fig. 2 D). 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. 2 D). Temporal analysis revealed the enrichment of lipid metabolism, granulocyte differentiation, and c-Jun N-terminal kinase (JNK) pathway pathways in the terminal stage (Fig. 2 D). 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. 2 D). 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. 2.3 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 (Table S2). We further examined the prognostic clustering effect of these genes. Prognostic analysis across five datasets (TCGA, GSE14764, GSE26193, and GSE26712; 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. 3 A, B). Correlation analysis also indicated a significant high correlation between fibroblast content and fibroblast-associated marker genes ( COL3A1 , COL1A2 , COL1A1 , FN1 , and DCN ; 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. 3 C), 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. 3 D). These findings suggest that the fibroblast content in HGSOC samples may be highly correlated with tumor malignancy. 2.4 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[ 2 ]. 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 (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. 4 A; 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 )) (Figs. 4 B and 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. 4 C). To elucidate the relationship between CAFs and ovarian cancer treatment response, we further sub-classified CAFs and identified seven distinct subtypes (Fig. 4 D). 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. 4 E). 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. 4 F) and the yellow module highly correlated with the BR1-specific F_RAMP3 module (R = 0.97, P < 0.001; Fig. 4 F). 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 (Figs. 4 G and 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. 4 H). 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. 4 I). These findings suggest that CXCL12-positive CAF cells may be associated with chemotherapy resistance. 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 (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. 5 A; 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 (Figs. 5 B, C and S6). We then re-clustered the fibroblasts and identified five fibroblast subpopulations, including the F_CXCL12 subgroup (Figs. 5 D and 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. 5 E), with characteristic genes primarily enriched in TGF-beta, MAPK, and cytokine interaction signaling pathways (Figure S7). In both 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. 5 E). 2.5 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 (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. 6 A). Using these significantly differential genes and Cox regression coefficients, we performed a risk analysis across four additional datasets(GSE14764, GSE26712, TCGA and GSE9891, 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. 6 A, B). Through Friends analysis, we observed a high correlation of CXCL12 with other genes (Fig. 6 C). Additionally, using these 24 genes, we applied enrichment analysis in three treatment datasets (GSE33482, GSE114206 [ 20 ], and GSE189843 [ 15 ]) to validate the association between CAFs and chemotherapy response (Fig. 6 D; 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. 6 D). 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. 2.6 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, as well as their interaction with tumor cells via CXCR4 , 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. 7 A). As shown by qPCR, all three siRNA knockdown fragments significantly suppressed CXCR4 mRNA expression levels in SKOV3 cells ( P < 0.05; Fig. 7 A). 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. 7 B, C). 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 (GSM6506110) and chemotherapy-sensitive (GSM6506114) samples in GSE211956 dataset respectively (Fig. 7 D, 7 E, Table S7). In the chemotherapy-resistant samples, we identified 5 clusters (Fig. 7 D(i)). Cell type analysis revealed that cluster 4 is enriched with fibroblasts with high expression level both CXCL12 and DCN (Fig. 7 D(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. 7 D(iii)). In contrast, in the chemotherapy-sensitive samples, we identified 4 clusters (Fig. 7 E(i)). Among these, clusters 0, 1, and 3 are enriched with fibroblasts (Fig. 7 E(ii)). Cluster 1 is specifically enriched with KRT19-positive tumor cells (Fig. 7 E). 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. 7 E). 3. Discussion Solid-organ malignancies exhibit heterogeneity that enhances tumor cell survival and drug resistance through various molecular mechanisms [ 30 ]. Intratumoral heterogeneity evolves both spatially and temporally during tumor development, reprogramming the tumor microenvironment [ 31 ]. Understanding tumor clonality and its evolutionary path is crucial for comprehending tumor cell behavior and identifying driving factors [ 32 , 33 ]. Tumors harbor numerous genetic and epigenetic alterations, many of which complicate the interpretation of their functional effects on tumor evolution [ 34 ]. 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 [ 35 – 39 ].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 [ 21 ]. Batch effects between datasets require robust correction techniques to guarantee accurate cross-dataset comparisons and integrations [ 21 ]. 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 [ 22 ]. 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 [ 18 ]. 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 [ 23 – 25 ]. 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 [ 26 , 27 ]. 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. 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 [ 28 ]. 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 [ 25 ]. CXCL12-positive CAFs could promote cancer cell migration and invasion and upregulate the expression of PDL1 in bladder and pancreatic cancer cells [ 29 , 30 ]. 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 [ 31 ]. 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. 4. Materials and Methods 4.1 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 , Table S1). Raw gene expression matrices were processed using the Seurat R package (version 3.2.3) [ 32 ]. 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 [ 33 ]. 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%. 4.2 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. 4.3 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 [ 24 ]. 4.4 Tumor evolution analysis 4.4.1 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. $$\:{S}_{i,j}=\frac{{N}_{i,j}}{\sum\:_{i=1}^{m}\:{N}_{i,j}}*log(1+\frac{n}{\sum\:_{j=1}^{n}\:{N}_{i,j}})$$ 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: $$\:\text{T}=\frac{1}{2}{Q}_{2}+\frac{1}{4}({Q}_{1}+{Q}_{3})$$ 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. 4.4.2 Clustering and refinement Using the high-purity tumor cells obtained in the first step for each tumor subgroup, we employed consensus clustering method (SC3[ 34 ]) 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 [ 34 ], 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. 4.4.3 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. 4.4 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 [ 24 ]. 4.5 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 [ 35 ]. 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. [ 36 ] 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. 4.6 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. 4.7 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). 4.8 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. 4.9 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 [ 37 ], we calculated the functional correlations among genes linked to the prognosis of HGSOC and drug resistance. 4.10 Spatial transcriptome data analysis We collected two HGSCO spatial transcriptome datasets including chemotherapy treatment response (GSE189843 and GSE211956, 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 ) [ 28 ]. 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. 4.11 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 [ 38 ]. The microarray datasets underwent quantile normalization, while the RNA-seq dataset was quantified in terms of fragments per kilobase million(FPKM) 4.12 Cell transfection The human ovarian cancer cells SKOV3 were 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 hours 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. 4.13 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 2 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 minutes and 85°C for 2 minutes. 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 sec. The primer sequences used were as follows: CXCR4-h-F ACTACACCGAGGAAATGGGCT; CXCR4-h-R CCCACAATGCCAGTTAAGAAGA; GAPDH-h-F TGACAACTTTGGTATCGTGGAAGG; GAPDH-h-R AGGCAGGGATGATGTTCTGGAGAG. 4.14 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 PMSF༈Phenylmethanesulfonyl fluoride༉and protease inhibitors; sonicated on ice; and centrifuged at 4°C, 12,000 g for 10 minutes. 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 minutes, 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 hour at room temperature, the membrane was incubated withCXCR4 (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 minutes each, incubate the membrane with Anti-rabbit IgG HRP-linked Antibody(7074S, 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. 4.15 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 hours, the medium was replaced with complete medium containing 0, 50, 100, 200, or 300 ng/mL CXCL12 and cultured for 48 hours. 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 hours after plating. Six hours later, CXCL12 (200 ng/mL) was added, and the CCK8 assay was conducted using the same method. 5. Conclusions In this study, we aimed to dissect the complex tumor heterogeneity of high-grade serous ovarian cancer (HGSOC) and its implications on prognosis and chemotherapy response using advanced single-cell and spatial transcriptomics. By analyzing data from 28 HGSOC patients across five datasets, we developed an innovative text mining-based machine learning method to unravel the evolutionary dynamics of tumor cell functions. This approach revealed critical tumor-related genes and the varied microenvironmental compositions on which different tumor cell functions rely. Our findings highlighted a significant connection between increased tumor cell state heterogeneity and worse patient outcomes, including prognosis and treatment resistance. We validated these insights using additional spatial and bulk transcriptomic datasets, encompassing a total of 1,030 patients. We further identified that heightened intra- and inter-tumoral functional clonality is closely linked with the characteristics of cancer-associated fibroblasts (CAFs). Notably, the spatial proximity between CXCL12-positive CAFs and tumor cells, facilitated by the CXCL12/CXCR4 axis, emerged as a strong predictor of poor prognosis and chemotherapy resistance. Moreover, we developed a panel of 24 genes that are highly correlated with CXCL12-positive fibroblasts. This gene panel effectively predicts both prognosis and chemotherapy response in HGSOC patients. Our study underscores the critical role of tumor heterogeneity in therapeutic outcomes and offers new insights into potential biomarkers for personalized treatment strategies. Declarations Supplementary Materials The following supporting information can be downloaded online. Conflicts of Interest The authors have declared that no competing interest exists. Funding This study was supported by the Henan Province and Ministry of Health of Medical Science and Technology Program (SBGJ202302028 for Tingjie Wang and SBGJ202101009 for Yongjun Guo), Dalian Science and Technology Innovation Fund (2022JJ12SN049 for Jun Yang), and the Fundamental Research Funds for the Central Universities. Author contributions The study design and supervision were conducted by Yongjun Guo, Jun Yang, Jun Li, and Tingjie Wang. Data analysis, in vitro experiments, and manuscript writing were performed by Tingjie Wang, Lingxi Tian, and Ruitao Long. Clinical data collection, sorting, analysiss and multicolor immunofluorescence staining analysis were performed by Bing Wei, Cuiyun Zhang, Bo Wang, and Yougai Zhang. The patients’ clinical data and their therapeutic responses were verified and evaluated by Yougai Zhang and Xiaofei Zhu. The manuscript has been reviewed and approved by all authors. Acknowledgments The authors thank Dr. Mengyun Ke for the sample collection and storage. Data Availability Statement The datasets analyzed in this study have been deposited in the Gene Expression Omnibus (GEO) repository under the accession numbers in supplementary tables. References Coburn SB, Bray F, Sherman ME, Trabert B (2017) International patterns and trends in ovarian cancer incidence, overall and by histologic subtype. Int J Cancer 140(11):2451–2460 Chowdhury S, Kennedy JJ, Ivey RG, Murillo OD, Hosseini N, Song X, Petralia F, Calinawan A, Savage SR, Berry AB, Reva B, Ozbek U, Krek A, Ma W, da Veiga Leprevost F, Ji J, Yoo S, Lin C, Voytovich UJ, Huang Y, Lee SH, Bergan L, Lorentzen TD, Mesri M, Rodriguez H, Hoofnagle AN, Herbert ZT, Nesvizhskii AI, Zhang B, Whiteaker JR, Fenyo D, McKerrow W, Wang J, Schurer SC, Stathias V, Chen XS, Barcellos-Hoff MH, Starr TK, Winterhoff BJ, Nelson AC, Mok SC, Kaufmann SH, Drescher C, Cieslik M, Wang P, Birrer MJ, Paulovich (2023) A. G., Proteogenomic analysis of chemo-refractory high-grade serous ovarian cancer. Cell 186(16):3476–3498e35 Zheng X, Wang X, Cheng X, Liu Z, Yin Y, Li X, Huang Z, Wang Z, Guo W, Ginhoux F, Li Z, Zhang Z, Wang X (2023) Single-cell analyses implicate ascites in remodeling the ecosystems of primary and metastatic tumors in ovarian cancer. Nat Cancer 4(8):1138–1156 Matthews BG, Bowden NA, Wong-Brown MW (2021) Epigenetic Mechanisms and Therapeutic Targets in Chemoresistant High-Grade Serous Ovarian Cancer. Cancers (Basel) 13, (23) Silva R, Glennon K, Metoudi M, Moran B, Salta S, Slattery K, Treacy A, Martin T, Shaw J, Doran P, Lynch L, Jeronimo C, Perry AS, Brennan DJ (2023) Unveiling the epigenomic mechanisms of acquired platinum-resistance in high-grade serous ovarian cancer. Int J Cancer 153(1):120–132 Liu J, Dang H, Wang XW (2018) The significance of intertumor and intratumor heterogeneity in liver cancer. Exp Mol Med 50, (1), e416 Olbrecht S, Busschaert P, Qian J, Vanderstichele A, Loverix L, Van Gorp T, Van Nieuwenhuysen E, Han S, Van den Broeck A, Coosemans A, Van Rompuy AS, Lambrechts D, Vergote I (2021) High-grade serous tubo-ovarian cancer refined with single-cell RNA sequencing: specific cell subtypes influence survival and determine molecular subtype classification. Genome Med 13(1):111 Zhu JW, Charkhchi P, Akbari MR (2022) Potential clinical utility of liquid biopsies in ovarian cancer. Mol Cancer 21(1):114 Zellmer VR (2014) Evolving concepts of tumor heterogeneity. Cell bioscience 4:1–8 Polyak K (2014) Tumor heterogeneity confounds and illuminates: a case for Darwinian tumor evolution. Nat Med 20(4):344–346 Cordani M, Dando I, Ambrosini G, Gonzalez-Menendez P (2024) Signaling, cancer cell plasticity, and intratumor heterogeneity. Cell Commun Signal 22(1):255 Ciriello G, Magnani L, Aitken SJ, Akkari L, Behjati S, Hanahan D, Landau DA, Lopez-Bigas N, Lupianez DG, Marine JC, Martin-Villalba A, Natoli G, Obenauf AC, Oricchio E, Scaffidi P, Sottoriva A, Swarbrick A, Tonon G, Vanharanta S, Zuber J (2024) Cancer Evolution: A Multifaceted Affair. Cancer Discov 14(1):36–48 de Visser KE, Joyce JA (2023) The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 41(3):374–403 Ferri-Borgogno S, Zhu Y, Sheng J, Burks JK, Gomez JA, Wong KK, Wong STC, Mok SC (2023) Spatial Transcriptomics Depict Ligand-Receptor Cross-talk Heterogeneity at the Tumor-Stroma Interface in Long-Term Ovarian Cancer Survivors. Cancer Res 83(9):1503–1516 Stur E, Corvigno S, Xu M, Chen K, Tan Y, Lee S, Liu J, Ricco E, Kraushaar D, Castro P, Zhang J, Sood AK (2022) Spatially resolved transcriptomics of high-grade serous ovarian carcinoma. iScience 25(3):103923 Khatib S, Pomyen Y, Dang H, Wang XW (2020) Understanding the Cause and Consequence of Tumor Heterogeneity. Trends Cancer 6(4):267–271 Iacobuzio-Donahue CA, Litchfield K, Swanton C (2020) Intratumor heterogeneity reflects clinical disease course. Nat Cancer 1(1):3–6 Hegde PS, Chen DS (2020) Top 10 Challenges in Cancer Immunotherapy. Immunity 52(1):17–35 Ma L, Wang L, Khatib SA, Chang CW, Heinrich S, Dominguez DA, Forgues M, Candia J, Hernandez MO, Kelly M, Zhao Y, Tran B, Hernandez JM, Davis JL, Kleiner DE, Wood BJ, Greten TF, Wang XW (2021) Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma. J Hepatol 75(6):1397–1408 Veskimae K, Scaravilli M, Niininen W, Karvonen H, Jaatinen S, Nykter M, Visakorpi T, Maenpaa J, Ungureanu D, Staff S (2018) Expression Analysis of Platinum Sensitive and Resistant Epithelial Ovarian Cancer Patient Samples Reveals New Candidates for Targeted Therapies. Transl Oncol 11(5):1160–1170 Lahnemann D, Koster J, Szczurek E, McCarthy DJ, Hicks SC, Robinson MD, Vallejos CA, Campbell KR, Beerenwinkel N, Mahfouz A, Pinello L, Skums P, Stamatakis A, Attolini CS, Aparicio S, Baaijens J, Balvert M, Barbanson B, Cappuccio A, Corleone G, Dutilh BE, Florescu M, Guryev V, Holmer R, Jahn K, Lobo TJ, Keizer EM, Khatri I, Kielbasa SM, Korbel JO, Kozlov AM, Kuo TH, Lelieveldt BPF, Mandoiu II, Marioni JC, Marschall T, Molder F, Niknejad A, Raczkowska A, Reinders M, Ridder J, Saliba AE, Somarakis A, Stegle O, Theis FJ, Yang H, Zelikovsky A, McHardy AC, Raphael BJ, Shah SP, Schonhuth (2020) A., Eleven grand challenges in single-cell data science. Genome Biol 21, (1), 31 Yu L, Cao Y, Yang JYH, Yang P (2022) Benchmarking clustering algorithms on estimating the number of cell types from single-cell RNA-sequencing data. Genome Biol 23(1):49 Zhang L, Cascio S, Mellors JW, Buckanovich RJ, Osmanbeyoglu HU (2024) Single-cell analysis reveals the stromal dynamics and tumor-specific characteristics in the microenvironment of ovarian cancer. Commun Biol 7(1):20 Wang T, Xu C, Zhang Z, Wu H, Li X, Zhang Y, Deng N, Dang N, Tang G, Yang X, Shi B, Li Z, Li L, Ye K (2022) Cellular heterogeneity and transcriptomic profiles during intrahepatic cholangiocarcinoma initiation and progression. Hepatology 76(5):1302–1317 Hussain A, Voisin V, Poon S, Karamboulas C, Bui NHB, Meens J, Dmytryshyn J, Ho VW, Tang KH, Paterson J, Clarke BA, Bernardini MQ, Bader GD, Neel BG, Ailles LE (2020) Distinct fibroblast functional states drive clinical outcomes in ovarian cancer and are regulated by TCF21. J Exp Med 217, (8). Mao X, Xu J, Wang W, Liang C, Hua J, Liu J, Zhang B, Meng Q, Yu X, Shi S (2021) Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives. Mol Cancer 20(1):131 Cai J, Tang H, Xu L, Wang X, Yang C, Ruan S, Guo J, Hu S, Wang Z (2012) Fibroblasts in omentum activated by tumor cells promote ovarian cancer growth, adhesion and invasiveness. Carcinogenesis 33(1):20–29 Xun Z, Ding X, Zhang Y, Zhang B, Lai S, Zou D, Zheng J, Chen G, Su B, Han L, Ye Y (2023) Reconstruction of the tumor spatial microenvironment along the malignant-boundary-nonmalignant axis. Nat Commun 14(1):933 Feig C, Jones JO, Kraman M, Wells RJ, Deonarine A, Chan DS, Connell CM, Roberts EW, Zhao Q, Caballero OL, Teichmann SA, Janowitz T, Jodrell DI, Tuveson DA, Fearon DT (2013) Targeting CXCL12 from FAP-expressing carcinoma-associated fibroblasts synergizes with anti-PD-L1 immunotherapy in pancreatic cancer. Proc Natl Acad Sci U S A 110(50):20212–20217 Zhang Z, Yu Y, Zhang Z, Li D, Liang Z, Wang L, Chen Y, Liang Y, Niu H (2023) Cancer-associated fibroblasts-derived CXCL12 enhances immune escape of bladder cancer through inhibiting P62-mediated autophagic degradation of PDL1. J Exp Clin Cancer Res 42(1):316 Li X, Bu W, Meng L, Liu X, Wang S, Jiang L, Ren M, Fan Y, Sun H (2019) CXCL12/CXCR4 pathway orchestrates CSC-like properties by CAF recruited tumor associated macrophage in OSCC. Exp Cell Res 378(2):131–138 Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd;, Hao Y, Stoeckius M, Smibert P, Satija R Comprehensive Integration of Single-Cell Data. Cell 2019, 177, (7), 1888–1902 e21. Lun AT, McCarthy DJ, Marioni JC (2016) A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor. F1000Res 5:2122 Kiselev VY, Kirschner K, Schaub MT, Andrews T, Yiu A, Chandra T, Natarajan KN, Reik W, Barahona M, Green AR, Hemberg M (2017) SC3: consensus clustering of single-cell RNA-seq data. Nat Methods 14(5):483–486 Tosches MA, Naumann YTM (2018) Evolution of pallium, hippocampus, and cortical cell types revealed by single-cell transcriptomics in reptiles. Science 360(6391):881–888 Zhang K, Erkan EP, Jamalzadeh S, Dai J, Andersson N, Kaipio K, Lamminen T, Mansuri N, Huhtinen K, Carpen O, Hietanen S, Oikkonen J, Hynninen J, Virtanen A, Hakkinen A, Hautaniemi S, Vaharautio A (2022) Longitudinal single-cell RNA-seq analysis reveals stress-promoted chemoresistance in metastatic ovarian cancer. Sci Adv 8(8):eabm1831 Yu G, Li F, Qin Y, Bo X, Wu Y, Wang S (2010) GOSemSim: an R package for measuring semantic similarity among GO terms and gene products. Bioinformatics 26(7):976–978 Aran D, Hu Z, Butte AJ (2017) xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 18(1):220 Disclaimer/Publisher’s Note The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content Additional Declarations The authors declare no competing interests. Supplementary Files SupplementaryMaterial.docx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4827560","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":333659242,"identity":"f8b43dbc-1b9c-4d03-950e-6eccc0b96c4a","order_by":0,"name":"tingjie wang","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA8klEQVRIie2RsWrDMBCGTwjOi+KsNh7yCgKD6WCcV7ExKEuGjh0FAU95gHrpM3TMKCKiSQ+QwdCWgqcMztYhlKpjFttjofqmO/g/jrsD8Hj+Jpuhun0XCFqpgefFtIAg4CpVHRJTfTw/inqOgqSVirxQm6ZsOBI5ZfA32X8uDh1FNCLJuaIQ6NPrqGJgky5sHyLTJtnyLgQmxHlMyQyUCUNKMXJTtrynELFslkKa1SVLHrgmcoYi4rbRpAG3PsxR1qas+WBFjeCOvOeumNol3tnqvXzKi5V0r/y65cUy0GZUAWDlfY/j8V8CNZ3xeDye/80PJMhRMJJ9CewAAAAASUVORK5CYII=","orcid":"","institution":"","correspondingAuthor":true,"prefix":"","firstName":"tingjie","middleName":"","lastName":"wang","suffix":""}],"badges":[],"createdAt":"2024-07-30 09:36:16","currentVersionCode":1,"declarations":{"humanSubjects":false,"vertebrateSubjects":true,"conflictsOfInterestStatement":false,"humanSubjectEthicalGuidelines":false,"humanSubjectConsent":false,"humanSubjectClinicalTrial":false,"humanSubjectCaseReport":false,"vertebrateSubjectEthicalGuidelines":true},"doi":"10.21203/rs.3.rs-4827560/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4827560/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":61811185,"identity":"6b75b636-5b7d-4da0-8193-efc7fe264c15","added_by":"auto","created_at":"2024-08-05 20:23:52","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":939927,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eHigh-grade serous ovarian cancer transcriptome atlas.\u003c/strong\u003e (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.\u003c/p\u003e","description":"","filename":"image1.png","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/a077e5663a71e46b310d954b.png"},{"id":61811186,"identity":"97813775-c339-4945-810f-69fc2f127cc7","added_by":"auto","created_at":"2024-08-05 20:23:52","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1370420,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eTumor branching evolution reveals intratumor heterogeneity in high-grade serous ovarian cancer. \u003c/strong\u003e(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. (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.\u003c/p\u003e","description":"","filename":"image2.png","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/d78aa52ec661c77e57740993.png"},{"id":61811187,"identity":"9a64a355-f3cc-4ba0-974a-76bdc42e1929","added_by":"auto","created_at":"2024-08-05 20:23:52","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":390058,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eTumor branching evolution reveals that intertumoral heterogeneity and the proportion of fibroblasts promote the poor prognosis of high-grade serous ovarian cancer. \u003c/strong\u003e(A)\u003cstrong\u003e \u003c/strong\u003eOverall 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, *\u003cem\u003eP\u003c/em\u003e \u0026lt; 0.05, **\u003cem\u003eP\u003c/em\u003e \u0026lt; 0.01, and ***\u003cem\u003eP\u003c/em\u003e \u0026lt; 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, *\u003cem\u003eP\u003c/em\u003e \u0026lt; 0.05, **\u003cem\u003eP\u003c/em\u003e\u0026lt; 0.01, and ***\u003cem\u003eP\u003c/em\u003e \u0026lt; 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).\u003c/p\u003e","description":"","filename":"image3.png","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/1c7da1a1253ad43a7cb2fb85.png"},{"id":61811192,"identity":"7eed7f72-fcc4-4646-bae4-a276c60e4535","added_by":"auto","created_at":"2024-08-05 20:23:53","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":551459,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIntra- and intertumoral heterogeneity of cancer-associated fibroblasts. \u003c/strong\u003e(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 \u003cstrong\u003e(\u003c/strong\u003eCAFs). (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)\u003cstrong\u003e \u003c/strong\u003eGO 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 \u003cem\u003eMMP7\u003c/em\u003e and \u003cem\u003eELF3\u003c/em\u003e; EPI_3, epithelial subtype with high expression of \u003cem\u003eHES1\u003c/em\u003e and \u003cem\u003eCD24\u003c/em\u003e. (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.\u003c/p\u003e","description":"","filename":"image4.png","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/6bbb245d264db6b1981d851f.png"},{"id":61811189,"identity":"d89fd138-c8c0-4c6c-837b-12d0838a3d7a","added_by":"auto","created_at":"2024-08-05 20:23:53","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":623680,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eHeterogeneity of cancer-associated fibroblasts is associated with chemotherapy treatment outcomes. \u003c/strong\u003e(A) Platinum-free interval values in the three tumor evolution branches in dataset GSE165897. (B) UMAP and bar plot showing the major cell types and their origin. (C) Volcano plot showing the differential genes for cancer-associated fibroblast \u003cstrong\u003e(\u003c/strong\u003eCAF) subtypes. Upregulated genes are indicated in red, while downregulated ones are indicated in blue. (D) UMAP and bar plot showing the subtypes and origins of CAFs. (E) WGCNA results showing the gene modules of distinct CAF subtypes in GSE165897. (E) Heatmap showing the CAF subtype correlation between datasets GSE154600 and GSE165897.\u003c/p\u003e","description":"","filename":"image5.png","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/aaea8409d9fdc3b3a0f9142a.png"},{"id":61811188,"identity":"81cfc8c2-5911-48a4-a466-a904ad493806","added_by":"auto","created_at":"2024-08-05 20:23:53","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":684768,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCorrelation between 24 genes of CXCL12-positive fibroblasts and prognosis and drug resistance in high-grade serous ovarian cancer. \u003c/strong\u003e(A) Forest plot showing the risk prognosis results from 24 shared CXCL12-positive fibroblasts obtained from two single-cell samples via COX regression. (B) Overall survival curves showing the prognosis results with different cancer-associated fibroblast \u003cstrong\u003e(\u003c/strong\u003eCAF) risk scores in the four high-grade serous ovarian cancer cohorts using the weighted sum with cox coefficients regression of 24 genes. Statistical significance was calculated using the log-rank test. (C) Overall survival curves showing the prognosis result of the high riskcores genes. (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.\u003c/p\u003e","description":"","filename":"image6.png","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/8d5666d9e4da3723b1aa49d4.png"},{"id":61811800,"identity":"0cad79b6-c621-4114-8419-e68793f4badd","added_by":"auto","created_at":"2024-08-05 20:31:53","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":1899587,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFunction and spatial distribution characteristics of CXCL12-positive fibroblasts. \u003c/strong\u003e(A) The bar chart showing the silencing effect of the CXCL12 receptor gene \u003cem\u003eCXCR4\u003c/em\u003e (ns, not significant, *\u003cem\u003eP\u003c/em\u003e \u0026lt; 0.05, **\u003cem\u003eP\u003c/em\u003e \u0026lt; 0.01, and ***\u003cem\u003eP\u003c/em\u003e\u0026lt; 0.001 by t-test). (B) Western blot showing the silencing effect of \u003cem\u003eCXCR4\u003c/em\u003e 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 (*\u003cem\u003eP\u003c/em\u003e\u0026lt; 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-sensitive and chemotherapy-resistant samples. Scale bars = 50 µm.\u003c/p\u003e","description":"","filename":"image7.png","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/8fdec6370d9203576ee8ad30.png"},{"id":61812029,"identity":"8eb586f6-84be-426d-bdab-573ac0c0746b","added_by":"auto","created_at":"2024-08-05 20:39:56","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":6767298,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/49638624-4c00-42bd-8dab-03971f0db2c7.pdf"},{"id":61811190,"identity":"ac90af3f-2d88-44c8-82e5-22c9eeeb4194","added_by":"auto","created_at":"2024-08-05 20:23:53","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":4761719,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryMaterial.docx","url":"https://assets-eu.researchsquare.com/files/rs-4827560/v1/259d78192527e6c9eb2b0b92.docx"}],"financialInterests":"The authors declare no competing interests.","formattedTitle":"\u003cp\u003eDecoding the effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer through tumor evolution analysis\u003c/p\u003e","fulltext":[{"header":"1. Introduction","content":"\u003cp\u003eOvarian cancer (OC) ranks first in mortality among gynecologic malignancies [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e]. High-grade serous ovarian cancer (HGSOC) is the most common and deadly subtype of OC [\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. Over 75% of HGSOC patients are diagnosed at an advanced stage with extensive malignant ascites and omental metastasis [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]. Complete staging surgery followed by platinum-based chemotherapy is the standard treatment for HGSOC [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e]. However, over 25% of patients develop chemotherapy resistance within 6 months of initial treatment, and 70% experience recurrence within 2\u0026ndash;3 years, ultimately succumbing to acquired drug resistance [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. 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.\u003c/p\u003e \u003cp\u003eInter- and intratumoral heterogeneity is a crucial factor influencing patient prognosis and treatment outcomes [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]. 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) [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e].\u003c/p\u003e \u003cp\u003ePrevious 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[\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]. In HGSOC, it could be classified into four molecular subtypes: mesenchymal, immunoreactive, differentiated, and proliferative [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. Among these, the immunoreactive subtype demonstrates a better prognosis and treatment response [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. 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 [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e]. However, intratumoral heterogeneity causes drug resistance, leading to therapeutic failure [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]. 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.\u003c/p\u003e \u003cp\u003eTumor 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 [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e]. Natural selection drives clonal expansion, leading to the emergence of subclones with varying proliferative, migratory, and invasive capabilities [\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e]. 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 [\u003cspan additionalcitationids=\"CR12\" citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e]. 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.\u003c/p\u003e \u003cp\u003eThe TME, which includes stromal cells, immune cells, extracellular matrix, and a variety of soluble factors, can influence tumor progression and the response to therapy [\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e]. Single-cell and spatial transcriptomics technologies have significantly advanced the investigation of heterogeneity in HGSOC [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]. 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 [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]. 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 [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eUnderstanding the evolutionary trajectories of tumor and TME cells is crucial to develope more effective and personalized treatment strategies [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. The multifaceted intratumoral heterogeneity and the TME may shape the evolutionary dynamics of OC [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. Understanding the tumoral heterogeneity and evolutionary biodiversity could provide conceptual knowledge about the evolving tumor cell community to improve patient outcomes.\u003c/p\u003e"},{"header":"2. Results","content":"\u003cp\u003e \u003cem\u003e2.1 Construction of a multi-omics atlas of HGSOC including chemotherapy response characteristics using single-cell and spatial transcriptomics\u003c/em\u003e \u003c/p\u003e \u003cp\u003eTo 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.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA). This dataset comprises four single-cell datasets and two spatial transcriptomics datasets from GEO datasets (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.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eB). We also curated two bulk RNA-sequencing datasets that include chemotherapy treatment responses from GEO datasets, encompassing 24 samples (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 intra-tumoral 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.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eC). 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%). We used the CNV method and calculated the epithelial cell score to identify 18,273 individual tumor cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eC). 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\u0026rsquo;s tumor.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003e2.2 Inter- and intratumoral heterogeneity in high-grade serous ovarian cancer\u003c/h2\u003e \u003cp\u003ePatient 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.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eA). 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 (Figure S1).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eHaving established the clonal evolutionary relationships across different samples, we identified key genes within these branches and conducted functional analyses (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.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eB). 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 [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eFurther analysis revealed significant heterogeneity in the distribution of different microenvironmental cells across the branches, with CAFs notably enriched in the BR3 branch (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eC). 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.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD). Secreted phosphoprotein 1(\u003cem\u003eSPP1\u003c/em\u003e) is highly expressed in multiple tumor types and interacts with fibroblasts or T cells via pathways such as SPP1-CD44, correlating with tumor malignancy[\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e]. Type IIA topoisomerase(\u003cem\u003eTOP2A\u003c/em\u003e), 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.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD).\u003c/p\u003e \u003cp\u003eWe 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.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD). Cell velocity analysis showed enrichment of cell adhesion-related pathways at the onset and metabolic, chemokine, and MAPK pathways at the end (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD). 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.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD). Temporal analysis revealed the enrichment of lipid metabolism, granulocyte differentiation, and c-Jun N-terminal kinase (JNK) pathway pathways in the terminal stage (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD). 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.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD).\u003c/p\u003e \u003cp\u003eThese 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.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003e2.3 Tumor functional clonality is associated with patient prognosis and TME composition\u003c/h2\u003e \u003cp\u003eIn 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 (Table S2). We further examined the prognostic clustering effect of these genes. Prognostic analysis across five datasets (TCGA, GSE14764, GSE26193, and GSE26712; 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.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA, B). Correlation analysis also indicated a significant high correlation between fibroblast content and fibroblast-associated marker genes (\u003cem\u003eCOL3A1\u003c/em\u003e, \u003cem\u003eCOL1A2\u003c/em\u003e, \u003cem\u003eCOL1A1\u003c/em\u003e, \u003cem\u003eFN1\u003c/em\u003e, and \u003cem\u003eDCN\u003c/em\u003e; 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.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC), 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.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD). These findings suggest that the fibroblast content in HGSOC samples may be highly correlated with tumor malignancy.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003e2.4 Identification of a subpopulation of CXCL12-positive fibroblasts associated with chemotherapy resistance and poor prognosis\u003c/h2\u003e \u003cp\u003ePrevious research has demonstrated a strong correlation between the heterogeneity of TME composition, especially the infiltration of fibroblasts, and prognosis[\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. 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 (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.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA; 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 (\u003cem\u003eMMP7\u003c/em\u003e) and E74 like ETS transcription factor 3 (\u003cem\u003eELF3\u003c/em\u003e)), and EPI_2 (epithelial cells with a high level of hes family bHLH transcription factor 1 (\u003cem\u003eHES1\u003c/em\u003e)) (Figs.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB and 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.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eC). To elucidate the relationship between CAFs and ovarian cancer treatment response, we further sub-classified CAFs and identified seven distinct subtypes (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eD). 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.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eE). 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\u0026thinsp;=\u0026thinsp;0.88, P\u0026thinsp;\u0026lt;\u0026thinsp;0.001, Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eF) and the yellow module highly correlated with the BR1-specific F_RAMP3 module (R\u0026thinsp;=\u0026thinsp;0.97, P\u0026thinsp;\u0026lt;\u0026thinsp;0.001; Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eF). 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 (Figs.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eG and 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.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eH). Specifically, F_CXCL12 interacts with tumor epithelial cells via the ligand-receptor pair CXCL12/CXCR4 and through the secretion of placental growth factor (\u003cem\u003ePGF\u003c/em\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eI). These findings suggest that CXCL12-positive CAF cells may be associated with chemotherapy resistance.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo 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 (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.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA; 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 (Figs.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB, C and S6). We then re-clustered the fibroblasts and identified five fibroblast subpopulations, including the F_CXCL12 subgroup (Figs.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eD and S6). Using gene co-expression analysis, we also obtained a gene module that was highly correlated with the F_CXCL12 subgroup (R\u0026thinsp;=\u0026thinsp;0.85, P\u0026thinsp;\u0026lt;\u0026thinsp;0.001; Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eE), with characteristic genes primarily enriched in TGF-beta, MAPK, and cytokine interaction signaling pathways (Figure S7).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eIn both 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.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eE).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003e2.5 Identification of a gene set affecting prognosis and drug resistance in HGSOC\u003c/h2\u003e \u003cp\u003eAbove, 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 (Table S5). Subsequent Cox regression analysis using the GSE26193 dataset revealed that \u003cem\u003eDCN\u003c/em\u003e, \u003cem\u003eCXCL12\u003c/em\u003e, and \u003cem\u003eTNFAIP6\u003c/em\u003e were significantly positively associated with poor prognosis in ovarian cancer, with \u003cem\u003eDCN\u003c/em\u003e specifically identified as a fibroblast-characteristic gene (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA). Using these significantly differential genes and Cox regression coefficients, we performed a risk analysis across four additional datasets(GSE14764, GSE26712, TCGA and GSE9891, 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 \u003cem\u003eCXCL12\u003c/em\u003e and \u003cem\u003eDCN\u003c/em\u003e showed significant positive correlations with poor overall patient prognosis (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA, B). Through Friends analysis, we observed a high correlation of \u003cem\u003eCXCL12\u003c/em\u003e with other genes (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eC). Additionally, using these 24 genes, we applied enrichment analysis in three treatment datasets (GSE33482, GSE114206 [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e], and GSE189843 [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]) to validate the association between CAFs and chemotherapy response (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eD; 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.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eD). 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.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003e2.6 Spatial distribution characteristics of CXCL12-positive fibroblasts affect the clinical prognostic results of HGSOC\u003c/h2\u003e \u003cp\u003eGiven the strong correlation between CXCL12-positive fibroblasts and the treatment outcomes and prognoses of HGSOC, as well as their interaction with tumor cells via \u003cem\u003eCXCR4\u003c/em\u003e, we elucidated the interaction mechanisms of CXCL12-positive fibroblasts in vitro and in spatial transcriptomic data. In vitro, we silenced the CXCL12 receptor gene \u003cem\u003eCXCR4\u003c/em\u003e in the ovarian cancer cell line SKOV3 (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eA). As shown by qPCR, all three siRNA knockdown fragments significantly suppressed \u003cem\u003eCXCR4\u003c/em\u003e mRNA expression levels in SKOV3 cells (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05; Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eA). 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 (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05; Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eB, C). 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 (GSM6506110) and chemotherapy-sensitive (GSM6506114) samples in GSE211956 dataset respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eD, \u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eE, Table S7). In the chemotherapy-resistant samples, we identified 5 clusters (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eD(i)). Cell type analysis revealed that cluster 4 is enriched with fibroblasts with high expression level both \u003cem\u003eCXCL12\u003c/em\u003e and \u003cem\u003eDCN\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eD(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 \u003cem\u003eKRT19\u003c/em\u003e and \u003cem\u003eCXCR4\u003c/em\u003e (the receptor for \u003cem\u003eCXCL12\u003c/em\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eD(iii)). In contrast, in the chemotherapy-sensitive samples, we identified 4 clusters (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eE(i)). Among these, clusters 0, 1, and 3 are enriched with fibroblasts (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eE(ii)). Cluster 1 is specifically enriched with KRT19-positive tumor cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eE). However, in the regions dominated by KRT19-positive tumor cells, there is no high expression of \u003cem\u003eCXCR4\u003c/em\u003e, 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 \u003cem\u003eCXCL12/DCN/KRT19\u003c/em\u003e in chemotherapy-resistant samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eE).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e"},{"header":"3. Discussion","content":"\u003cp\u003eSolid-organ malignancies exhibit heterogeneity that enhances tumor cell survival and drug resistance through various molecular mechanisms [\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]. Intratumoral heterogeneity evolves both spatially and temporally during tumor development, reprogramming the tumor microenvironment [\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]. Understanding tumor clonality and its evolutionary path is crucial for comprehending tumor cell behavior and identifying driving factors [\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e, \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e]. Tumors harbor numerous genetic and epigenetic alterations, many of which complicate the interpretation of their functional effects on tumor evolution [\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e]. 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 [\u003cspan additionalcitationids=\"CR36 CR37 CR38\" citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e].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.\u003c/p\u003e \u003cp\u003eSingle-cell data processing faces challenges with large-scale datasets, necessitating efficient computational and statistical methods to handle data sparsity and ensure informative analysis [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e]. Batch effects between datasets require robust correction techniques to guarantee accurate cross-dataset comparisons and integrations [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e]. 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 [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]. In this study, we initially identified intratumoral heterogeneity using results from Seurat\u0026rsquo;s initial clustering combined with CNV data. We extracted tumor-associated feature genes and identified \u0026ldquo;representative cells\u0026rdquo; 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.\u003c/p\u003e \u003cp\u003eUsing 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 [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e]. 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 [\u003cspan additionalcitationids=\"CR24\" citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]. 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 [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e, \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. Combining with WGCNA method, we identified a clinical predictive model comprising 24 genes and discovered that fibroblasts with high \u003cem\u003eCXCL12\u003c/em\u003e expression are significantly related to poor prognosis and chemotherapy resistance. 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 [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e]. 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 [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]. CXCL12-positive CAFs could promote cancer cell migration and invasion and upregulate the expression of PDL1 in bladder and pancreatic cancer cells [\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e, \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]. 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 [\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]. 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.\u003c/p\u003e"},{"header":"4. Materials and Methods","content":"\u003cdiv id=\"Sec10\" class=\"Section2\"\u003e \u003ch2\u003e4.1 Preprocessing of single-cell RNA-seq data\u003c/h2\u003e \u003cp\u003eWe collected single-cell RNAseq matrix from five HGSOC datasets (GSE140819, GSE154600, GSE165897, GSE184880 and data in \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://blueprint.lambrechtslab.org\u003c/span\u003e\u003cspan address=\"http://blueprint.lambrechtslab.org\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e, Table S1). Raw gene expression matrices were processed using the Seurat R package (version 3.2.3) [\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e]. 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 [\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e]. 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%.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003e4.2 Identification of tumor cells\u003c/h2\u003e \u003cp\u003eCopy 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. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/broadinstitute/inferCNV\u003c/span\u003e\u003cspan address=\"https://github.com/broadinstitute/inferCNV\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). 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 \u0026minus;\u0026thinsp;1 to 1. The CNV for each cell was calculated as the sum of the squares of the CNV values across the respective regions.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003e4.3 Calculation of epithelial cell scores\u003c/h2\u003e \u003cp\u003eWe used the expression levels of epithelial scores based on the mean expression of five well-established markers: \u003cem\u003eKRT19\u003c/em\u003e, \u003cem\u003eKRT7\u003c/em\u003e, \u003cem\u003eKRT18\u003c/em\u003e, \u003cem\u003eKRT8\u003c/em\u003e, and \u003cem\u003eEPCAM\u003c/em\u003e, as reported in previous studies [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e].\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003e4.4 Tumor evolution analysis\u003c/h2\u003e \u003cdiv id=\"Sec14\" class=\"Section3\"\u003e \u003ch2\u003e4.4.1 Selection of representative tumor cells\u003c/h2\u003e \u003cp\u003eWe 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\u0026thinsp;\u0026lt;\u0026thinsp;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.\u003c/p\u003e \u003cp\u003eWe 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.\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:{S}_{i,j}=\\frac{{N}_{i,j}}{\\sum\\:_{i=1}^{m}\\:{N}_{i,j}}*log(1+\\frac{n}{\\sum\\:_{j=1}^{n}\\:{N}_{i,j}})$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eWhere m and n represent the number of genes and cells in the single-cell expression matrix, respectively. \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{N}_{i,j}\\)\u003c/span\u003e\u003c/span\u003e is the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:m\\)\u003c/span\u003e\u003c/span\u003e by \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:n\\)\u003c/span\u003e\u003c/span\u003e 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 \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}_{i,j}\\)\u003c/span\u003e\u003c/span\u003e, 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 \u0026times; A. Any gene weight below this threshold in \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}_{i,j}\\)\u003c/span\u003e\u003c/span\u003e was set to 0, resulting in a corrected \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}_{i,j}\\)\u003c/span\u003e\u003c/span\u003e matrix. Subsequently, for each cell group, we summed the weights of the characteristic genes to derive a total weight and set the cutoff value \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\text{T}\\)\u003c/span\u003e\u003c/span\u003e to establish the cell weight threshold utilizing the formula:\u003cdiv id=\"Equb\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e\n$$\\:\\text{T}=\\frac{1}{2}{Q}_{2}+\\frac{1}{4}({Q}_{1}+{Q}_{3})$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003ewhere 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.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section3\"\u003e \u003ch2\u003e4.4.2 Clustering and refinement\u003c/h2\u003e \u003cp\u003eUsing the high-purity tumor cells obtained in the first step for each tumor subgroup, we employed consensus clustering method (SC3[\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e]) 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\u0026rsquo;s subclassification obtained through the consensus clustering method SC3 [\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e], 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.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section3\"\u003e \u003ch2\u003e4.4.3 Hierarchical evolutionary analysis\u003c/h2\u003e \u003cp\u003eWe 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\u0026rsquo;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\u0026thinsp;=\u0026thinsp;1000).\u003c/p\u003e \u003cp\u003eFollowing the aforementioned step, we obtained highly purified subpopulations within tumors, facilitating the analysis of intratumoral heterogeneity.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003e4.4 Cell-type enrichment analysis\u003c/h2\u003e \u003cp\u003eEach 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 [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e].\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003e4.5 Weighted correlation network analysis in cell subclusters and acquisition of fibroblast-related Genes\u003c/h2\u003e \u003cp\u003eThe 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 \u0026ldquo;pseudo cells,\u0026rdquo; which represent the average gene expression of 10 randomly selected cells within each distinct cell type [\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e]. Network construction was achieved using the \u0026ldquo;blockwiseModules\u0026rdquo; 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\u0026rsquo;s correlation test. Subsequently, hub genes within significant modules were identified based on their modular connectivity, which refers to the absolute value of Pearson\u0026rsquo;s correlation between genes (module membership) and their relationship with clinical traits, defined as the absolute value of Pearson\u0026rsquo;s correlation between individual gene expression and cell type. We performed WGCNA in the fibroblast subpopulations of datasets GSE154600 and GSE165897 from GEO database. [\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e] 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.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003e4.6 Similarity analysis among cell subpopulations\u003c/h2\u003e \u003cp\u003eWe 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 (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/hannet91/ggcor\u003c/span\u003e\u003cspan address=\"https://github.com/hannet91/ggcor\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) to compute the correlations between cell subpopulations and between cell subpopulations and phenotypes and visualize these correlations.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003e4.7 Clustering of bulk RNA data samples\u003c/h2\u003e \u003cp\u003eBulk RNA data were retrieved from the Cancer Genome Atlas (TCGA; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.cancer.gov/tcga\u003c/span\u003e\u003cspan address=\"https://www.cancer.gov/tcga\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) 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).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec21\" class=\"Section2\"\u003e \u003ch2\u003e4.8 Mapping of single-cell and Bulk RNA subpopulations\u003c/h2\u003e \u003cp\u003eTo 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.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec22\" class=\"Section2\"\u003e \u003ch2\u003e4.9 Friends analysis\u003c/h2\u003e \u003cp\u003eThe 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 [\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e], we calculated the functional correlations among genes linked to the prognosis of HGSOC and drug resistance.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec23\" class=\"Section2\"\u003e \u003ch2\u003e4.10 Spatial transcriptome data analysis\u003c/h2\u003e \u003cp\u003eWe collected two HGSCO spatial transcriptome datasets including chemotherapy treatment response (GSE189843 and GSE211956, 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 (\u003cem\u003ePTPRC\u003c/em\u003e, \u003cem\u003eCD2\u003c/em\u003e, \u003cem\u003eCD3D\u003c/em\u003e, \u003cem\u003eCD3E\u003c/em\u003e, \u003cem\u003eCD3G\u003c/em\u003e, \u003cem\u003eCD5\u003c/em\u003e, \u003cem\u003eCD7\u003c/em\u003e, \u003cem\u003eCD79A\u003c/em\u003e, \u003cem\u003eMS4A1\u003c/em\u003e, and \u003cem\u003eCD19\u003c/em\u003e) [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e]. 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; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://ludvigla.github.io/STUtility_web_site/\u003c/span\u003e\u003cspan address=\"https://ludvigla.github.io/STUtility_web_site/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e), 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.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec24\" class=\"Section2\"\u003e \u003ch2\u003e4.11 Immune infiltration analysis\u003c/h2\u003e \u003cp\u003eWe 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 [\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e]. The microarray datasets underwent quantile normalization, while the RNA-seq dataset was quantified in terms of fragments per kilobase million(FPKM)\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec25\" class=\"Section2\"\u003e \u003ch2\u003e4.12 Cell transfection\u003c/h2\u003e \u003cp\u003eThe human ovarian cancer cells SKOV3 were purchased from the ATCC cell bank.\u003c/p\u003e \u003cp\u003ewere digested with trypsin and resuspended. After cell counting, cells were plated at a density of 4 \u0026times; 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 \u0026micro;L lipofectamine per well. After 48 hours of incubation, cells were collected for qRT-PCR and Western blot. The sequences of siRNA for the knockdown and control groups were as follows:\u003c/p\u003e \u003cp\u003esiCXCR4-1 GGCAAUGGAUUGGUCAUCCUGGUCA;\u003c/p\u003e \u003cp\u003esiCXCR4-2 UGGUUGGCCUUAUCCUGCCUGGUAU;\u003c/p\u003e \u003cp\u003esiCXCR4-3 UGUUUCCACUGAGUCUGAGUCUUCA;\u003c/p\u003e \u003cp\u003esiNC UUCUCCGAACGUGUCACGUTT.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec26\" class=\"Section2\"\u003e \u003ch2\u003e4.13 qRT-PCR\u003c/h2\u003e \u003cp\u003eqRT-PCR was used to detect the transcriptional level of the \u003cem\u003eCXCR4\u003c/em\u003e 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 \u0026micro;L Rnase-free H\u003csub\u003e2\u003c/sub\u003eO. The concentration and A260/280 of RNA samples were measured using a Microplate UV-Vis Spectrophotometer. Reverse transcription was performed according to the Goldenstar\u0026reg; RT6 cDNA Synthesis Kit Ver.2 ༈TSK302S, Tsingke, Beijing, China)instructions, and 2 \u0026micro;g of total RNA was used as the template. The reaction mixture was incubated at 50\u0026deg;C for 5 minutes and 85\u0026deg;C for 2 minutes. 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\u0026deg;C 30 s; cycling stage (40 cycles), step 1 95\u0026deg;C 15 s, step 2 60\u0026deg;C 30 s; and melt curve stage: 95\u0026deg;C 15 s, 60\u0026deg;C 60 s, 95\u0026deg;C 15 sec. The primer sequences used were as follows:\u003c/p\u003e \u003cp\u003eCXCR4-h-F ACTACACCGAGGAAATGGGCT;\u003c/p\u003e \u003cp\u003eCXCR4-h-R CCCACAATGCCAGTTAAGAAGA;\u003c/p\u003e \u003cp\u003eGAPDH-h-F TGACAACTTTGGTATCGTGGAAGG;\u003c/p\u003e \u003cp\u003eGAPDH-h-R AGGCAGGGATGATGTTCTGGAGAG.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec27\" class=\"Section2\"\u003e \u003ch2\u003e4.14 Western blot\u003c/h2\u003e \u003cp\u003eWestern 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 PMSF༈Phenylmethanesulfonyl fluoride༉and protease inhibitors; sonicated on ice; and centrifuged at 4\u0026deg;C, 12,000 \u003cem\u003eg\u003c/em\u003e for 10 minutes. Then, the supernatant was taken for protein quantification using the BCA (Bicinchoninic Acid Assay )method. After mixing with loading buffer and denaturing at 100\u0026deg;C for 5 minutes, 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 hour at room temperature, the membrane was incubated withCXCR4 (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 minutes each, incubate the membrane with Anti-rabbit IgG HRP-linked Antibody(7074S, 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.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec28\" class=\"Section2\"\u003e \u003ch2\u003e4.15 CCK8 assay\u003c/h2\u003e \u003cp\u003eThe CCK8 assay was used to detect the effect of CXCL2/CXCR4 on cell viability. When cells reached 80\u0026ndash;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\u0026deg;C incubator for 24 hours, the medium was replaced with complete medium containing 0, 50, 100, 200, or 300 ng/mL CXCL12 and cultured for 48 hours. CCK8 reagent (10 \u0026micro;L) was added per well, and the absorbance at 450 nm was measured using a microplate reader after a 1-hour incubation at 37\u0026deg;C. In separate experiments, siCXCR4 and siNC transfections were performed 24 hours after plating. Six hours later, CXCL12 (200 ng/mL) was added, and the CCK8 assay was conducted using the same method.\u003c/p\u003e \u003c/div\u003e"},{"header":"5. Conclusions","content":"\u003cp\u003eIn this study, we aimed to dissect the complex tumor heterogeneity of high-grade serous ovarian cancer (HGSOC) and its implications on prognosis and chemotherapy response using advanced single-cell and spatial transcriptomics. By analyzing data from 28 HGSOC patients across five datasets, we developed an innovative text mining-based machine learning method to unravel the evolutionary dynamics of tumor cell functions. This approach revealed critical tumor-related genes and the varied microenvironmental compositions on which different tumor cell functions rely. Our findings highlighted a significant connection between increased tumor cell state heterogeneity and worse patient outcomes, including prognosis and treatment resistance. We validated these insights using additional spatial and bulk transcriptomic datasets, encompassing a total of 1,030 patients. We further identified that heightened intra- and inter-tumoral functional clonality is closely linked with the characteristics of cancer-associated fibroblasts (CAFs). Notably, the spatial proximity between CXCL12-positive CAFs and tumor cells, facilitated by the CXCL12/CXCR4 axis, emerged as a strong predictor of poor prognosis and chemotherapy resistance. Moreover, we developed a panel of 24 genes that are highly correlated with CXCL12-positive fibroblasts. This gene panel effectively predicts both prognosis and chemotherapy response in HGSOC patients. Our study underscores the critical role of tumor heterogeneity in therapeutic outcomes and offers new insights into potential biomarkers for personalized treatment strategies.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e \u003ch2\u003eSupplementary Materials\u003c/h2\u003e \u003cp\u003eThe following supporting information can be downloaded online.\u003c/p\u003e \u003c/p\u003e\u003cp\u003e \u003ch2\u003eConflicts of Interest\u003c/h2\u003e \u003cp\u003eThe authors have declared that no competing interest exists.\u003c/p\u003e \u003c/p\u003e\u003ch2\u003eFunding\u003c/h2\u003e \u003cp\u003eThis study was supported by the Henan Province and Ministry of Health of Medical Science and Technology Program (SBGJ202302028 for Tingjie Wang and SBGJ202101009 for Yongjun Guo), Dalian Science and Technology Innovation Fund (2022JJ12SN049 for Jun Yang), and the Fundamental Research Funds for the Central Universities.\u003c/p\u003e\u003ch2\u003eAuthor contributions\u003c/h2\u003e \u003cp\u003eThe study design and supervision were conducted by Yongjun Guo, Jun Yang, Jun Li, and Tingjie Wang. Data analysis, in vitro experiments, and manuscript writing were performed by Tingjie Wang, Lingxi Tian, and Ruitao Long. Clinical data collection, sorting, analysiss and multicolor immunofluorescence staining analysis were performed by Bing Wei, Cuiyun Zhang, Bo Wang, and Yougai Zhang. The patients\u0026rsquo; clinical data and their therapeutic responses were verified and evaluated by Yougai Zhang and Xiaofei Zhu. The manuscript has been reviewed and approved by all authors.\u003c/p\u003e\u003ch2\u003eAcknowledgments\u003c/h2\u003e \u003cp\u003eThe authors thank Dr. Mengyun Ke for the sample collection and storage.\u003c/p\u003e\u003ch2\u003eData Availability Statement\u003c/h2\u003e \u003cp\u003eThe datasets analyzed in this study have been deposited in the Gene Expression Omnibus (GEO) repository under the accession numbers in supplementary tables.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eCoburn SB, Bray F, Sherman ME, Trabert B (2017) International patterns and trends in ovarian cancer incidence, overall and by histologic subtype. Int J Cancer 140(11):2451\u0026ndash;2460\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChowdhury S, Kennedy JJ, Ivey RG, Murillo OD, Hosseini N, Song X, Petralia F, Calinawan A, Savage SR, Berry AB, Reva B, Ozbek U, Krek A, Ma W, da Veiga Leprevost F, Ji J, Yoo S, Lin C, Voytovich UJ, Huang Y, Lee SH, Bergan L, Lorentzen TD, Mesri M, Rodriguez H, Hoofnagle AN, Herbert ZT, Nesvizhskii AI, Zhang B, Whiteaker JR, Fenyo D, McKerrow W, Wang J, Schurer SC, Stathias V, Chen XS, Barcellos-Hoff MH, Starr TK, Winterhoff BJ, Nelson AC, Mok SC, Kaufmann SH, Drescher C, Cieslik M, Wang P, Birrer MJ, Paulovich (2023) A. G., Proteogenomic analysis of chemo-refractory high-grade serous ovarian cancer. Cell 186(16):3476\u0026ndash;3498e35\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZheng X, Wang X, Cheng X, Liu Z, Yin Y, Li X, Huang Z, Wang Z, Guo W, Ginhoux F, Li Z, Zhang Z, Wang X (2023) Single-cell analyses implicate ascites in remodeling the ecosystems of primary and metastatic tumors in ovarian cancer. Nat Cancer 4(8):1138\u0026ndash;1156\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMatthews BG, Bowden NA, Wong-Brown MW (2021) Epigenetic Mechanisms and Therapeutic Targets in Chemoresistant High-Grade Serous Ovarian Cancer. \u003cem\u003eCancers (Basel)\u003c/em\u003e 13, (23)\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSilva R, Glennon K, Metoudi M, Moran B, Salta S, Slattery K, Treacy A, Martin T, Shaw J, Doran P, Lynch L, Jeronimo C, Perry AS, Brennan DJ (2023) Unveiling the epigenomic mechanisms of acquired platinum-resistance in high-grade serous ovarian cancer. Int J Cancer 153(1):120\u0026ndash;132\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu J, Dang H, Wang XW (2018) The significance of intertumor and intratumor heterogeneity in liver cancer. Exp Mol Med 50, (1), e416\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eOlbrecht S, Busschaert P, Qian J, Vanderstichele A, Loverix L, Van Gorp T, Van Nieuwenhuysen E, Han S, Van den Broeck A, Coosemans A, Van Rompuy AS, Lambrechts D, Vergote I (2021) High-grade serous tubo-ovarian cancer refined with single-cell RNA sequencing: specific cell subtypes influence survival and determine molecular subtype classification. Genome Med 13(1):111\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhu JW, Charkhchi P, Akbari MR (2022) Potential clinical utility of liquid biopsies in ovarian cancer. Mol Cancer 21(1):114\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZellmer VR (2014) Evolving concepts of tumor heterogeneity. Cell bioscience 4:1\u0026ndash;8\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePolyak K (2014) Tumor heterogeneity confounds and illuminates: a case for Darwinian tumor evolution. Nat Med 20(4):344\u0026ndash;346\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCordani M, Dando I, Ambrosini G, Gonzalez-Menendez P (2024) Signaling, cancer cell plasticity, and intratumor heterogeneity. Cell Commun Signal 22(1):255\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCiriello G, Magnani L, Aitken SJ, Akkari L, Behjati S, Hanahan D, Landau DA, Lopez-Bigas N, Lupianez DG, Marine JC, Martin-Villalba A, Natoli G, Obenauf AC, Oricchio E, Scaffidi P, Sottoriva A, Swarbrick A, Tonon G, Vanharanta S, Zuber J (2024) Cancer Evolution: A Multifaceted Affair. Cancer Discov 14(1):36\u0026ndash;48\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ede Visser KE, Joyce JA (2023) The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 41(3):374\u0026ndash;403\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFerri-Borgogno S, Zhu Y, Sheng J, Burks JK, Gomez JA, Wong KK, Wong STC, Mok SC (2023) Spatial Transcriptomics Depict Ligand-Receptor Cross-talk Heterogeneity at the Tumor-Stroma Interface in Long-Term Ovarian Cancer Survivors. Cancer Res 83(9):1503\u0026ndash;1516\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eStur E, Corvigno S, Xu M, Chen K, Tan Y, Lee S, Liu J, Ricco E, Kraushaar D, Castro P, Zhang J, Sood AK (2022) Spatially resolved transcriptomics of high-grade serous ovarian carcinoma. iScience 25(3):103923\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKhatib S, Pomyen Y, Dang H, Wang XW (2020) Understanding the Cause and Consequence of Tumor Heterogeneity. Trends Cancer 6(4):267\u0026ndash;271\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eIacobuzio-Donahue CA, Litchfield K, Swanton C (2020) Intratumor heterogeneity reflects clinical disease course. Nat Cancer 1(1):3\u0026ndash;6\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHegde PS, Chen DS (2020) Top 10 Challenges in Cancer Immunotherapy. Immunity 52(1):17\u0026ndash;35\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMa L, Wang L, Khatib SA, Chang CW, Heinrich S, Dominguez DA, Forgues M, Candia J, Hernandez MO, Kelly M, Zhao Y, Tran B, Hernandez JM, Davis JL, Kleiner DE, Wood BJ, Greten TF, Wang XW (2021) Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma. J Hepatol 75(6):1397\u0026ndash;1408\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVeskimae K, Scaravilli M, Niininen W, Karvonen H, Jaatinen S, Nykter M, Visakorpi T, Maenpaa J, Ungureanu D, Staff S (2018) Expression Analysis of Platinum Sensitive and Resistant Epithelial Ovarian Cancer Patient Samples Reveals New Candidates for Targeted Therapies. Transl Oncol 11(5):1160\u0026ndash;1170\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLahnemann D, Koster J, Szczurek E, McCarthy DJ, Hicks SC, Robinson MD, Vallejos CA, Campbell KR, Beerenwinkel N, Mahfouz A, Pinello L, Skums P, Stamatakis A, Attolini CS, Aparicio S, Baaijens J, Balvert M, Barbanson B, Cappuccio A, Corleone G, Dutilh BE, Florescu M, Guryev V, Holmer R, Jahn K, Lobo TJ, Keizer EM, Khatri I, Kielbasa SM, Korbel JO, Kozlov AM, Kuo TH, Lelieveldt BPF, Mandoiu II, Marioni JC, Marschall T, Molder F, Niknejad A, Raczkowska A, Reinders M, Ridder J, Saliba AE, Somarakis A, Stegle O, Theis FJ, Yang H, Zelikovsky A, McHardy AC, Raphael BJ, Shah SP, Schonhuth (2020) A., Eleven grand challenges in single-cell data science. \u003cem\u003eGenome Biol\u003c/em\u003e 21, (1), 31\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYu L, Cao Y, Yang JYH, Yang P (2022) Benchmarking clustering algorithms on estimating the number of cell types from single-cell RNA-sequencing data. Genome Biol 23(1):49\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang L, Cascio S, Mellors JW, Buckanovich RJ, Osmanbeyoglu HU (2024) Single-cell analysis reveals the stromal dynamics and tumor-specific characteristics in the microenvironment of ovarian cancer. Commun Biol 7(1):20\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang T, Xu C, Zhang Z, Wu H, Li X, Zhang Y, Deng N, Dang N, Tang G, Yang X, Shi B, Li Z, Li L, Ye K (2022) Cellular heterogeneity and transcriptomic profiles during intrahepatic cholangiocarcinoma initiation and progression. Hepatology 76(5):1302\u0026ndash;1317\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHussain A, Voisin V, Poon S, Karamboulas C, Bui NHB, Meens J, Dmytryshyn J, Ho VW, Tang KH, Paterson J, Clarke BA, Bernardini MQ, Bader GD, Neel BG, Ailles LE (2020) Distinct fibroblast functional states drive clinical outcomes in ovarian cancer and are regulated by TCF21. J Exp Med 217, (8).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMao X, Xu J, Wang W, Liang C, Hua J, Liu J, Zhang B, Meng Q, Yu X, Shi S (2021) Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives. Mol Cancer 20(1):131\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCai J, Tang H, Xu L, Wang X, Yang C, Ruan S, Guo J, Hu S, Wang Z (2012) Fibroblasts in omentum activated by tumor cells promote ovarian cancer growth, adhesion and invasiveness. Carcinogenesis 33(1):20\u0026ndash;29\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXun Z, Ding X, Zhang Y, Zhang B, Lai S, Zou D, Zheng J, Chen G, Su B, Han L, Ye Y (2023) Reconstruction of the tumor spatial microenvironment along the malignant-boundary-nonmalignant axis. Nat Commun 14(1):933\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFeig C, Jones JO, Kraman M, Wells RJ, Deonarine A, Chan DS, Connell CM, Roberts EW, Zhao Q, Caballero OL, Teichmann SA, Janowitz T, Jodrell DI, Tuveson DA, Fearon DT (2013) Targeting CXCL12 from FAP-expressing carcinoma-associated fibroblasts synergizes with anti-PD-L1 immunotherapy in pancreatic cancer. Proc Natl Acad Sci U S A 110(50):20212\u0026ndash;20217\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang Z, Yu Y, Zhang Z, Li D, Liang Z, Wang L, Chen Y, Liang Y, Niu H (2023) Cancer-associated fibroblasts-derived CXCL12 enhances immune escape of bladder cancer through inhibiting P62-mediated autophagic degradation of PDL1. J Exp Clin Cancer Res 42(1):316\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi X, Bu W, Meng L, Liu X, Wang S, Jiang L, Ren M, Fan Y, Sun H (2019) CXCL12/CXCR4 pathway orchestrates CSC-like properties by CAF recruited tumor associated macrophage in OSCC. Exp Cell Res 378(2):131\u0026ndash;138\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eStuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd;, Hao Y, Stoeckius M, Smibert P, Satija R Comprehensive Integration of Single-Cell Data. Cell 2019, 177, (7), 1888\u0026ndash;1902 e21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLun AT, McCarthy DJ, Marioni JC (2016) A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor. F1000Res 5:2122\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKiselev VY, Kirschner K, Schaub MT, Andrews T, Yiu A, Chandra T, Natarajan KN, Reik W, Barahona M, Green AR, Hemberg M (2017) SC3: consensus clustering of single-cell RNA-seq data. Nat Methods 14(5):483\u0026ndash;486\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTosches MA, Naumann YTM (2018) Evolution of pallium, hippocampus, and cortical cell types revealed by single-cell transcriptomics in reptiles. Science 360(6391):881\u0026ndash;888\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang K, Erkan EP, Jamalzadeh S, Dai J, Andersson N, Kaipio K, Lamminen T, Mansuri N, Huhtinen K, Carpen O, Hietanen S, Oikkonen J, Hynninen J, Virtanen A, Hakkinen A, Hautaniemi S, Vaharautio A (2022) Longitudinal single-cell RNA-seq analysis reveals stress-promoted chemoresistance in metastatic ovarian cancer. Sci Adv 8(8):eabm1831\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYu G, Li F, Qin Y, Bo X, Wu Y, Wang S (2010) GOSemSim: an R package for measuring semantic similarity among GO terms and gene products. Bioinformatics 26(7):976\u0026ndash;978\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAran D, Hu Z, Butte AJ (2017) xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 18(1):220\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDisclaimer/Publisher\u0026rsquo;s Note The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"The Affiliated Cancer Hospital of Zhengzhou University","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Tumor cell state, Functional clonality, Tumor evolution, HGSOC, Tumor transcriptomic heterogeneity, Microenvironment, Spatiotemporal transcriptome, CXCL12, Cancer-associated fibroblasts","lastPublishedDoi":"10.21203/rs.3.rs-4827560/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4827560/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eTumor heterogeneity is associated with poor prognosis and drug resistance, leading to therapeutic failure. Here, we used tumor evolution analysis to determine the intra- and intertumoral heterogeneity of high-grade serous ovarian cancer (HGSOC) and analyze the correlation between tumor heterogeneity and prognosis, as well as chemotherapy response, through single-cell and spatial transcriptomic analysis. We collected and curated 28 HGSOC patients\u0026rsquo; single-cell transcriptomic data from five datasets. Then, we developed a novel text-mining-based machine-learning approach to deconstruct the evolutionary patterns of tumor cell functions. We then identified key tumor-related genes within different evolutionary branches, characterized the microenvironmental cell compositions that various functional tumor cells depend on, and analyzed the intra- and intertumoral heterogeneity as well as the tumor microenvironments. These analyses were conducted in relation to the prognosis and chemotherapy response in HGSOC patients. We validated our findings in two spatial and seven bulk transcriptomic datasets (total: 1,030 patients). Using transcriptomic clusters as proxies for functional clonality, we identified a significant increase in tumor cell state heterogeneity that was strongly correlated with patient prognosis and treatment response. Furthermore, increased intra- and intertumoral functional clonality was associated with the characteristics of cancer-associated fibroblasts (CAFs). The spatial proximity between CXCL12-positive CAFs and tumor cells, mediated through the CXCL12/CXCR4 interaction, was highly positively correlated with poor prognosis and chemotherapy resistance in HGSOC. In this study, we constructed a panel of 24 genes through statistical modeling that correlate with CXCL12-positive fibroblasts and can predict both prognosis and the response to chemotherapy in HGSOC patients.\u003c/p\u003e","manuscriptTitle":"Decoding the effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer through tumor evolution analysis","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-08-05 20:23:47","doi":"10.21203/rs.3.rs-4827560/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"82bae4c2-2e34-402f-b3b9-3986627eaed5","owner":[],"postedDate":"August 5th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2024-08-16T08:59:38+00:00","versionOfRecord":[],"versionCreatedAt":"2024-08-05 20:23:47","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-4827560","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4827560","identity":"rs-4827560","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}